Today we are going to look at Problem 40 from the book 50 Challenging Problems in Probability by Frederick Mosteller. However, instead of solving it analytically, we are going to use the Monte Carlo method (random sampling) to solve it.

The problem goes like this:

Shuffle an ordinary deck of 52 playing cards containing four aces. Then turn up cards from the top until the first ace appears. On the average, how many cards are required to produce the first ace?

If you are dumbfounded, worry not, because so was I when I first saw this problem. It may seem like there is no expedient mathamatical method that would lead us to the answer. So what do we do? This is where Monte Carlo method comes in. You probably know that if you were to toss a coin, there would be a 50-50 chance of it being heads or tails. Now to verify that fact, what you would do, assuming you have a lot of free time, would be to simply toss it many times, and see what percentage of the total tosses ended up heads or tails. It would be unlikely for the percentages to be exactly 50% tails and 50% heads, but after tossing the coin many times, they would undoubtedly be very close to 50-50. The above is a case of using Monte Carlo methods to evaluate a problem deterministic in nature via stochastic sampling.

How can we apply Monte Carlo methods here? Obviously cannot shuffle and draw cards in any reasonable amount of time in real life to obtain sufficient sampling. Luckily, we have the power of modern computers + programming to help us. In other words, we are going to simulate the situation using computers. It is easy to get sufficient sampling with computers, since computers are quite good at doing repetitive tasks. Here we will use python and numpy (a numerical package for python) to build our simulation.

First, let us make our deck as an object in python:

import numpy as np

class deck():
    def __init__(self,deck_size):
        self.deck_size=deck_size
        self.cards=np.zeros(deck_size,dtype='bool')
        self.cards[:4]=True
    def shuffle_deck(self):
        np.random.shuffle(self.cards)
    def draw_until_ace(self):
        for ace_draw in range(self.deck_size):
            if self.cards[ace_draw]:
                break
        return ace_draw+1

It is a pretty straightforward object with two attributes and two methods. deck_size is just the size of the deck, in our case, 52. Since all we are only concerned with drawing aces, we can consider every single card of the deck as ace or not ace. Therefore, cards is just a vector of 52 binaries, 4 of which are True, because we have 4 aces in our deck. it does not matter which four to begin with, since the deck will be shuffled. shuffle_deck allows us to shuffle the deck using np.random.shuffle. The method draw_until_ace returns the number of cards it takes to get an ace in a given shuffle. Now we have a deck that can be shuffled much faster than any human can in real life, so let’s get on with the simulation part.

monte_carlo_deck=deck(52) #initializes the deck
total=0 #total number of draws 
count=0 #total number of shuffles      
while True:
    monte_carlo_deck.shuffle_deck() #shuffle the deck
    ace_draw=monte_carlo_deck.draw_until_ace() #draw until an ace appears
    total+=ace_draw #update number of draws 
    count+=1 #update number of shuffles
    print('Average draws taken to get ace is {:f}'.format(total/count),
          end='\r',flush=True)

You will see that the average draws taken converges very quickly to around 10.6, with the following terminal output:
Average draws taken to get ace is 10.601732.

The book solves this problem using principle of symmetry, where on average the 4 aces cut the rest of the deck into portions of 9.6 non-ace cards each. Therefore, it takes on average 10.6 cards to get an ace. We see that our answer obtained with Monte Carlo methods coincide with the solution provided by the book. Our little simulation proves successful!