More Monte Carlo Simulations with Python

Following on from my previous post surveying some of the concepts underpinning statistics and probability, I thought it might be interesting to apply some Monte Carlo simulations to Roulette, as a means of consolidating some of the topics covered.

Suppose we have a friend, named Adam, who wants to visit a casino, in either Las Vegas or Monte Carlo this weekend. He intends to gamble but is equipped with nothing other than a burning hole in his wallet and a keen desire ‘to beat the house’. In his great and unmatched wisdom, he has approached you to help him. You accept. With the knowledge acquired in the last post, you rush to your laptop and begin figuring out the odds associated with Roulette.

Part 1

Immediately alarm bells begin to ring in your ears. As was established in the previous article, we know that over a large number of trials, the Law of Large Numbers states that the results will tend towards the expected value. It is reasonable to assume that this is favourable to the casino and is in fact the basis of its profits. But before we crush Adam’s dreams, let’s figure out some of the probabilities associated with roulette and then apply them empirically in Monte Carlo simulations.

To begin with, let’s clarify the distinction between European Roulette and American Roulette. European Roulette wheels have numbers from 0 to 36. The ball landing on 0 is always a house win. Thus there are 37 potential spots for the ball to land. American Roulette wheels are very similar, except that it includes a 00 as well. Thus there are 38 total spots for the ball to land. There you have it, the great continental divide between roulette wheels. Let’s take a look at the implied probabilities on the sample space.

European Roulette:

  • High (>= 19) /Low (<= 18) bet @ odds of 1:1. = (1 x 1837) – (1 x 1937) = -0.02702
  • Red/Black bet @ odds of 1:1 = (1 x 1837) – (1 x 1937) = -0.02702
  • Single Number bet @ 35:1 = (35 x 137) – (1 x 3637) = -0.02702

American Roulette:

  • High (>= 19) /Low (<= 18) bet @ odds of 1:1. = (1 x 1838) – (1 x 2038) = -0.05263
  • Red/Black bet @ odds of 1:1 = (1 x 1838) – (1 x 2038) = -0.05263
  • Single Number bet @ 35:1 = (35 x 138) – (1 x 3738) = -0.05263

Simulations

While this is far from an exhaustive analysis of the bets made at Roulette tables, it does suggest, over the long term at least, that the best destination to play Roulette is Monte Carlo. Best, of course, is relative. The expected return is actually a loss of 2.7p for every £1 bet. Still, we know that the expected value is actually the long term weighted average and that it is still possible to make money at the roulette wheels.

Ever the empiricist, Adam asks you to create some simulations. You construct the most basic model for a single bet of £10 on high numbers(This offers exactly the same odds and probabilities as Red/Black):

num = random.randint(0,36)
if num >= 19:
    print("You earn £20!")
else:
    print("You lose £10!")

You lose £10!

While it seems like you lose money more frequently than you gain money, this is not sufficient. So, you build a simulation of the same bet with a limit of 1000 games:

bankroll = 100

for i in range(1000):
    bankroll -= 10
    num = random.randint(0,36)
    if num >= 19:
        bankroll += 20

print("Bankroll after playing 1000 games: £%d" %bankroll)

Bankroll after playing 1000 games: £460

Great news! You have more than quadrupled your money. However, one has the sneaking suspicion that, aside from the vast amount of time required to spin a roulette wheel 1000 times (The average roulette wheel spins 45 times an hour in US Casinos. 1000 times would likely be just under a day spent sat the table - Uh oh!), there is no guarantee that the bankroll didn’t go into serious negative territory before it came back out. No, Adam demands rigour and visualisations!

bankroll = 100
high_walk_results = []
games = 1


while bankroll >= 0 and bankroll <= 1000 and games <= 1000:
    games += 1
    bankroll -= 10
    num = random.randint(0,36)
    if num >= 19:
        bankroll += 20
        high_walk_results.append(bankroll)
    else:
        high_walk_results.append(bankroll)

print("Bankroll after playing %d games : £%d" %(games, bankroll))

Bankroll after playing 126 games : £-10

This time we have an upper limit of £1000 and a lower limit if the bankroll goes below 0. The game is capped at 1000 spins too. This simulates a disciplined player and ends with a loss after 126 spins:

high walk capped

Although interesting, it might be more insightful to visualise more than one game. This time, without an upper limit on the bankroll and only 100 spins:

multi high walk

While four games out of five resulted in a loss of your bankroll, one was cut off too soon! Adam’s palms get sweaty. How about just the ceiling for the bankroll:

multi high walk bank ceiling

This doesn’t look too promising. Adam insists that you simulate a single number bet, on the basis of a larger potential reward for a successful spin. Undeterred you do:

single single walk

Ahh, that is quite disheartening. But it is only a single example after all. Let’s try many:

multi single walk

That is quite gloomy.

Let’s calculate the expected value based on the results. The same formula could be used as is used to calculate the expected value on the sample space, replacing the implied probability of profit/loss with the relative frequency.

However, in the interest of simplicity, loss per unit staked can be used. This is simply the total result divided by the total result. For high numbers:

result = []
stakes = []

for i in range(1000000):
    num = random.randint(0,36)
    stakes.append(10)
    if num >= 19:
        result.append(10)
    else:
        result.append(-10)

sum(result) / sum(stakes)

-0.02682

And for a single number:

result = []
stakes = []

for i in range(1000000):
    num = random.randint(0,36)
    stakes.append(10)
    if num == 8:
        result.append(350)
    else:
        result.append(-10)

sum(result) / sum(stakes)

-0.025408

Lo and behold: the expected value does not deviate greatly from that calculated on the sample space: -0.027.

Strategies

In the meantime, Adam has been reading about various strategies that people online have been using to make a profit on roulette. Casinos hate them! Before you can remind him about the law of large numbers, he insists that you run a number of simulations following the Big-Martingale and D’Alembert strategies. You sigh and turn back to the laptop.

Both strategies are fairly easy to understand and rely on 18 number probability patterns (think High/Low or Red/Black) at odds of 1:1.

D’Alembert Strategy

The D’Alembert Strategy begins with an initial stake of 1 unit (1 chip, £10, $10 etc.). If you lose a round, you increase the stake by 1 unit. If you win a round, you decrease the stake by 1 unit. That’s it.

So you write the following simulation quickly.

def multi_dalembert_roulette(stake=10, bankroll=100, total_result=[], games=0):
    for i in range(10):
        stake=10
        bankroll=100
        games=0
        result = []
        while bankroll >= 0:
            games += 1
            bankroll -= stake
            num = random.randint(0,36)
            if num >= 19 and stake > 10:
                bankroll += stake * 2
                stake = stake - 10
                result.append(bankroll)
            elif num >= 19 and stake == 10:
                bankroll += stake * 2
                stake = 10
                result.append(bankroll)
            else:
                stake = stake + 10
                result.append(bankroll)
        total_result.append(result)
    return total_result

And here visualise the results:

multi_dalembert_walk

How about after 100 spins?

dalembert dist

What about the loss per unit staked?

def expected_value_dalembert_roulette(stake=10, result=[], stakes=[], games=1):
        while games <= 1000000:
            games += 1
            num = random.randint(0,36)
            stakes.append(stake)
            if num >= 19 and stake > 10:
                result.append(stake)
                stake = stake - 10
            elif num >= 19 and stake == 10:
                result.append(stake)
                stake = 10
            else:
                result.append(0-stake)
                stake = stake + 10
        return stakes, result

dalembert_stakes, dalembert_results = expected_value_dalembert_roulette()

sum(dalembert_results) / sum(dalembert_stakes)

-0.0270554915380498

As suggested by the expected value.

Big-Martingale Strategy

The Big-Martingale Strategy also begins with an initial stake of 1 unit. If you lose, you double the stake. Unsurprisingly, consecutive losses snowball very quickly. Below we create a generator to simulate the growing size of ones losses after consecutive unsuccessful spins following the Big-Martingle Strategy.

def big_martingale_sim():

    i = 1

    while True:
        yield i
        i = i * 2

b = big_martingale_sim()

big_martingale_list = [next(b) for i in range(15)]

bigm growth of loss

Still, you promised Adam a simulation.

def multi_big_martingale_roulette(stake=10,bankroll=100, end_result=[], games=1):
    for i in range(10):
        stake=10
        bankroll=100
        games=1
        result = []
        while bankroll >= 0:
            bankroll -= stake
            num = random.randint(0,36)
            games += 1
            if num >= 19:
                bankroll += stake * 2
                stake = 10
                result.append(bankroll)
            else:
                stake = (stake * 2)
                result.append(bankroll)
        end_result.append(result)
    return end_result

Unsurprisingly, the results do not suggest the favourability of either of these strategies.

bigm strat 10 bigm strat dist

Again we use the loss per unit staked as a proxy for the expected value,

def expected_value_bigm(stake=10, result=[], stakes=[], games=1):
        while games <= 1000000:
            games += 1
            num = random.randint(0,36)
            stakes.append(stake)
            if num >= 19:
                result.append(stake)
                stake = 10
            else:
                result.append(0-stake)
                stake = stake * 2
        return (result, stakes)

bigm_results, bigm_stakes = expected_value_bigm()

sum(bigm_results) / sum(bigm_stakes)

-0.027331623484638207

and again the result is -0.027.

Conclusion

This has been a very educational experience for both you and Adam. He has learnt about the Law of Large Numbers and it’s basis as the profitability for casino games. That over the long run, you are practically guaranteed to lose your money at the roulette wheels. You have learnt that he may well have a nascent gambling issue and endeavour to suggest a healthier way to spend his holidays. Perhaps a weekend at the seaside?