Monte Carlo in Python

Today we look at a very famous method called the Monte Carlo in Python, which can be used to solve any problem having a probabilistic interpretation. Let’s get started with some fun history.

A Bit of Monte Carlo history

It was used to solve complex numerical problems like the Buffon Needle Problem (https://en.wikipedia.org/wiki/Buffon%27s_needle_problem):

Suppose we have a floor made of parallel strips of wood, each the same width, and we drop a needle onto the floor. What is the probability that the needle will lie across a line between two strips?

And has been used since the 1940s for:

• studying neutron diffusion for nuclear weapons projects at the Los Alamos National Laboratory, where the ENIAC was used to perform M-C simulations.
• was also used during the development of the hydrogen bomb
• in fluid mechanics, to solve complex differential equations (nonlinear parabolic PDEs)
• for estimating particle transmission energies
• in advanced signal processing and Bayesian inference
• on genetic type mutation-selection learning machines (an early introduction to today’s field of bioinformatics), etc.

Implementing Monte Carlo in Python

Monte Carlo allows us to simulate seemingly random events, and assess risks (among other results, of course). It has been used to assess the risk of a given trading strategy.

For this tutorial, we will simulate a casino (because we cannot simulate a nuclear weapons test 😀 )

In general, the more complex the problem, the more pseudo-random numerical variables you would need.

Open up your Google Colaboratory, and connect to runtime.

1. Creating the basic roll of a casino wheel

Let’s import our numpy and pandas packages:

import numpy as np
import pandas as pd

• Then we define our “roll” as a number from 1 to 100, and let’s set it at 49-51 odds of winning for the customers.
• This means that for the rolls 1-50 and exactly 100, the house (/casino) wins.
• For rolls 51-99, the player wins.

This seems like a fair chance of winning, and most players will most likely take it.

Let’s simulate this in Python:

def casino_roll():
    roll = np.random.randint(1,100)

    if roll <= 50 or roll == 100:
        print('rolled %d. you lose. House wins! Play again?!'%roll)
    elif 50 < roll < 100:
        print('you win! On a roll !')
        return True

So, now we can just call casino_roll(), and you can see what you get. I tried it 4 times, and I lost thrice.

2. Making a dealer

Next, let’s put some money into the bet. So we’ll create the dealer, who takes in a bet, and:

• if there is a win, he rewards the player.
• If there is a loss, he pockets the money.

def deal_bet(funds,initial_bet,bet_count):
    bet = 0
    funds_total = 0
    while bet < bet_count:
        if casino_roll():
            funds += initial_bet
        else:
            funds -= initial_bet

        bet += 1
        print('Funds : %d'%funds)
        funds_total += funds
    print('average win/loss : %f'%(10000-funds_total/bet_count))

Observe that I found the total funds, and then found the average win/loss. If it’s positive it’s a win. If it’s negative, it’s a loss.

3. For 100 bets of Rs. 100

So running this for a hundred bets, each time betting the same amount of 100 rupees, we get:

deal_bet(10000,100,100)

I had to run the program 5 times before getting this negative value, but observe how many of the values are well over 10000 and yet the player loses money overall.

This is just to show that what we perceive as a profitable deal may not always be so.

4. Making more players

We can make 100 players make this same above bet:

x=0
while x<100:
  deal_bet(10000,100,100)
  x+=1

Comment out the print(funds) statement in deal_bet().

So, now we can see all the profits and losses incurred by each player:

5. Plotting averages with Matplotlib

Let’s plot the data with matplotlib here:

return (10000-funds_total/bet_count)

Add the above line to end of deal_bet() function.

Then we modify:

x=0
avg_list=[]
while x<100:
  avg_list.append(deal_bet(10000,100,100))
  x+=1

And finally, we can plot it:

import matplotlib.pyplot as plt

avg_list = np.array(avg_list).astype(int)
plt.figure(figsize=(20,20))
plt.plot(avg_list)
plt.axhline(np.array(avg_list).mean(), color='k', linestyle='dashed', linewidth=1)

Conclusion

That’s it for the day everyone. I hope this example has helped you understand the Monte Carlo simulation in Python perfectly well.

You could now do this for an aggressive player. Perhaps, every time he wins, he bets double money. Think of a variety of scenarios and modify the betting function yourself.

You’ll find some more scenarios along with the full code if you’ve missed something at my github:

https://github.com/arkaprabha-majumdar/monte-carlo

Until next time !!