I have been reading a lovely book by Paul Nahin, Digital Dice. The book is subtitled Computational Solutions to Practical Probability Problems. I don’t know how practical the problems are, but they are a lot of fun.

Probability theory is rife with seemingly simple problems which turn out to be less simple than first appears. When possible, I like to approach such problems from two sides: analytically and computationally. If I can come up with two solutions, I can test one against the other. If the solutions don’t agree, at least one of them is wrong. Nahin’s book encourages this approach, and provides some fine practice material. For example, Nahin offers this problem, taken from a 1970’s textbook:

A pipe smoker has two booklets of matches in his pocket, each containing 40 matches initially. Whenever a match is required, he picks one of the booklets at random, removing one match. Write a program using random numbers to simulate the situation 100 times and determine the average number of matches that can be removed until one booklet is completely empty.

This is easily simulated in terms of coin tosses. We repeatedly flip a fair two-sided coin; heads represents one book of matches, tails the other. We’ll tally the outcome of each coin toss, and the trial is over as soon as we see the 40th head or the 40th tail, whichever comes first. The outcome of the trial is the total number of coin tosses. What we want to know is the average number of coin tosses per trial.

We’ll do the simulation in Clojure. First we’ll need to model a coin toss:

(defn coin-toss []
  (if (zero? (rand-int 2))
    [1 0]
    [0 1]))

Here’s how we use the coin to simulate one trial.

(def *N* 40)

(defn trial
  ([] (trial 0 0))
  ([heads, tails]
     (if (or (= *N* tails) (= *N* heads))
       (+ heads tails)
       (let [[h t] (coin-toss)]
         (recur (+ heads h) (+ tails t))))))

We need some code to conduct repeated trials, gather the results, and compute the average.

(defn average [coll]
  (float (/ (reduce + coll) (count coll))))

(defn trials [n]
  (average (take n (repeatedly trial))))

With a simulation in hand, let's think about solving the problem mathematically. For fixed N (40, in our case), we will repeatedly flip a coin until we see the Nth head or the Nth tail. Clearly a trial will require at least N coin tosses. The most tosses required will be 2N - 1. That happens when we manage to toss N-1 heads and N-1 tails, for a total of 2N - 2 tosses. The next toss will be decisive, no matter what the outcome, ending the trial at 2N - 1 tosses.

So every trial ends in N + K tosses, where 0 <= K < N. Write P(N, K) for the probability that the trial ends after exactly N+K tosses. The expected value we are trying to compute is the sum of P(N,K) * (N+K), taken over 0 <= K < N. The task boils down to computing P(N,K). How can we calculate P(N,K)? Suppose the trial ends with N heads and K tails. The final toss had to have been heads. Then we must have had K tails out of first N+K-1 tosses. How many ways can that happen? The answer is well known: (N+K-1)!/[( N-1)! K!]. What is the probability of of any one of those outcomes? Given a fair coin it is 1/2^(N+K).

Let's implement the calculation so we can compare it to the results of the earlier simulation.

(defn factorial [n]
  (reduce * (range 1 (inc n))))

(defn choose [n k]
  (/ (factorial n)
     (factorial (- n k))
     (factorial k)))

(defn prob [n k]
  (/ (choose (dec (+ n k)) k)
     (clojure.contrib.math/expt 2 (+ n k))))

(defn expected-outcome []
  (let [n *N*
        probs (for [k (range n)] (prob n k))
        tosses (range n (+ n n))]
    (float (reduce + (map * probs tosses)))))

Now that we have both the simulation code and the calculation code, the only thing needed is a bit of scaffolding to compare the two approaches.

(def *num-trials* 10000)

(defn comparison [n]
  (binding [*N* n]
       [*N* (avg-outcome *num-trials*) (expected-outcome)]))

The comparison function is parameterized so we can vary the problem a bit. First let's try the original N=40.

user> (comparison 40)
[40 72.7887 36.442886]

The simulation and the calculation differ by a factor of about two. Why? Look back to the line that reads "Suppose the trial ends with N heads and K tails". I left out half the cases, the ones that end with N tails and K heads. The thing to do is correct the prob function:

(defn prob [n k]
  (/ (choose (dec (+ n k)) k)
     (clojure.contrib.math/expt 2 (dec (+ n k)))))

After that correction the comparison function shows that the simulation and the calculation agree nicely across a range of values of N:

user> (map comparison (range 20 60 5 ))
([20 34.9941 34.985172] [25 44.3824 44.386242] [30 53.7716 53.84531] [35 63.3641 63.348217] [40 72.8602 72.88577] [45 82.5092 82.4516] [50 92.0544 92.04108] [55 101.6011 101.65071])

In my experience, when the simulation and the calculation disagree, the simulation is usually more trustworthy. Simulations tend to be straightforward, while analytic approaches can be tricky to get right.