The question is: If you flip coins in a row, what is the probability that somewhere in those flips you’ll get a series of (or more) heads in a row?

This is one that comes up periodically on Reddit, and I’ve answered it, though I’m not 100% I’ve always answered it exactly the same way twice. I do know that my answers to this one tend to be a little long and messy. To make a long story short, I always use some version of a recursive approach.

So the purpose of this article is to introduce the idea of recursive probability calculations, for problems which don’t seem amenable to another approach.

It interests me because it seems like there should be a simple formulation, but I always have to work out a fairly complicated approach. I’ve observed over the years that with counting problems there are almost always multiple approaches, and often one that’s quite elegant and simple. I suspect that’s the case here, but if so the alternate methods have so far eluded me.

Anyway, my approach is to use recursion, and if you haven’t seen a recursive probability calculation, then you may find my long-way around useful.

Why is this so hard?

First let’s talk about what makes it hard. Let’s say and . You might be tempted to list the possible outcomes with heads as: HHHxx, xHHHx, and xxHHH, where x means you don’t care what that flip was. Can’t you just count all of the outcomes fitting each of those patterns?

Well, not exactly. For instance HHHxx includes HHHHT, HHHHH, HHHTT and HHHHH. While xHHHx includes HHHHT, THHHT, HHHHH and THHHH. You’ll notice that two of those outcomes, HHHHH and HHHHT, are counted twice, and the pattern xxHHH includes other outcomes that have already been counted. So you wouldn’t just add up the four HHHxx outcomes, the four xHHHX, and the four xxHHH outcomes to say there are twelve ways this can happen. You’ve overcounted, which is always a danger in counting problems.

Now, there are well-established procedures for handling overcounting. First of all, you can try to construct your outcomes so that they don’t overlap. To make your outcomes mutually exclusive in other words. So we could use the pattern HHHTx and count those separately from HHHHx. There is guaranteed to be no overlap between HHHTx and xHHHx or xxHHH. With care, you could construct mutually-exclusive classes to count the outcomes for .

But when I try to do that for larger problems, say , where you could have multiple runs of 3 which are separated or together and a huge number of mutually-exclusive patterns, it gets unwieldy. I’m not saying it can’t be done, just saying it breaks my brain and I haven’t found a way yet.

The other approach to overcounting is inclusion-exclusion. Suppose you want to count the number of things that fall in set A, or in set B or in set C, etc and those sets have overlap as those three patterns HHHxx, xHHHx, and xxHHH do. Then the inclusion-exclusion rule is that to calculate the total number of things in the union of those sets:

• Add up the membership of each individual set, . We know we’re overcounting the overlaps, so next we
• Subtract the membership of each intersecting pair of sets, . Each of those pairs is also counting things that belong to three or more sets, so we’ve over done the correction, which means we have to…
• Add the membership of triplets of sets,
• etc, alternating addition and subtraction (four sets, five sets, etc), until we’ve handled every combination.

This seems to me even worse than trying to construct mutually-exclusive categories. Suppose there were different positions where the HHH could occur. Then you’d have to construct all the combinations of etc outcomes and count each of those. In effect one term for every subset of a collection of items. There are subsets.

Again, there may be a different collection of events you could define and a clever inclusion-exclusion argument that’s tractable. But I haven’t found it yet. So that brings us to…

Recursive counting of coin flips

Recursive definition of anything is a pretty simple idea. You have a rule for defining higher terms (larger ) in terms of lower terms (smaller ). And then you define what the first terms are so the rule has a place to start.

The classic example used in teaching recursive computer programming is the factorial function. can be defined recursively as , and you could define to start things off. So for instance

As another example, the Fibonacci numbers can be defined recursively as , so each term is defined in terms of two previous terms. It thus requires two numbers to be defined to start it out, so we usually define .

For instance, to use the recursive expression to evaluate we would do this: . If we were evaluating we would have to follow the chain from and all the way down to and . It can take a lot of computation to evaluate such an expression for large values of , though a few hundred should not be a problem for any computer (Python limits recursion to a depth of by default.)

Let’s define as the probability of at least one run of heads in flips, and let’s see what probabilities we can easily calculate for low values of and to use as the starting values for the recursions.

First of all, for any . Every flip is a “run” of length 1.

for all because you can’t have a run of length in less than flips. We can include negative values of as well.

is the probability that all flips are heads and is therefore equal to for any .

That handles the cases of and . Now let’s consider a general . There will be a run of length in the flips if either

• There is a run of within the first , probability , or
• The last flips are THH…H, probability , and there is no run of length in the first , probability

This gives us the recursion

Let’s see how this works with, for example, . The outcomes with 3 or more heads in a row are HHHTT, HHHTH, HHHHT, HHHHH, THHHT, THHHH, HTHHH, TTHHH, altogether out of . So .

Recursively:

Here’s an implementation of the recursive calculation in Python.

def P(n, k):
    if n < k:
        return 0
    elif n == k:
        return 1/2**n
    else:
        return P(n-1, k) + (1 - P(n-k-1, k))/2**(k+1)

and here is a function to generate all possible sequences of n coin flips and collect statistics on the longest sequence of heads found (the logical value True is coding for heads):

import itertools as it
import numpy as np
def seq_stats(n):
    
    # All possible flips of length n
    all_n_seqs = it.product({True, False}, repeat=n)
    
    stats = np.zeros((n+1,), dtype=int)
    for flips in all_n_seqs:
        seq_len = 0
        longest_seq = 0
        for flip in flips:
            if flip:
                seq_len += 1
                longest_seq = max(longest_seq, seq_len)
            else:
                seq_len = 0
          
        # print(flips)
        # print(f'Longest seq = {longest_seq}')
        stats[longest_seq] += 1
    
    return stats

This code runs a side-by-side comparison of the entire probability distribution for any n (takes about 10 seconds for n = 20).

# Do a side by side comparison of the two methods
def compare_methods(n):
    t0 = time.perf_counter()
    flip_stats = seq_stats(n)
    # stats are P(>=k) so do the cumulative sum
    empirical_stats = np.cumsum(flip_stats[::-1])[::-1] / 2**n
    t1 = time.perf_counter()
    theoretical_stats = [P(n,k) for k in range(n+1)]
    t2 = time.perf_counter()
    print('  k Empirical P(n,k)  Theoretical P(n,k)')
    for k in range(n+1):
        print(f'{k:3d} {empirical_stats[k]:12.5f} {theoretical_stats[k]:12.5f}')

and finally, here are the results for , validating the recursive formula:

  k Empirical P(n,k)  Theoretical P(n,k)
  0      1.00000      1.00000
  1      0.99902      0.99902
  2      0.85938      0.85938
  3      0.50781      0.50781
  4      0.24512      0.24512
  5      0.10938      0.10938
  6      0.04688      0.04688
  7      0.01953      0.01953
  8      0.00781      0.00781
  9      0.00293      0.00293
 10      0.00098      0.00098