Exact probability for coin flip streaks

Question

In another question the user was determining the probability of having a streak of 6 heads or 6 tails in 100 coin flips. To find the probability they would generate 100 random coin flips and determine if there was a streak. They would test 10,000 such sequences of 100 flips to find that there was about an 80% chance of there being a streak in 100 coin flips.

I decided to compute the exact probability. For 100 flips there are \$2^{100}\$ possible outcomes. To determine the percentage I compute how many of them have a streak, and then divide by \$2^{100}\$.

My naive solution gets me the number for 20 flips in few seconds:

from itertools import product

def naive(flips, streak):
    return sum('h' * streak in ''.join(p) or
               't' * streak in ''.join(p)
               for p in product('ht', repeat=flips))

Result:

>>> naive(20, 6)
248384

My fast solution gets me the number for 100 flips instantly:

from collections import Counter

def fast(flips, streak):
    needles = 'h' * streak, 't' * streak
    groups = {'-' * streak: 1}
    total = 0
    for i in range(flips):
        next_groups = Counter()
        for ending, count in groups.items():
            for coin in 'ht':
                new_ending = ending[1:] + coin
                if new_ending in needles:
                    total += count * 2**(flips - 1 - i)
                else:
                    next_groups[new_ending] += count
        groups = next_groups
    return total

The idea is to have a pool of still ongoing games, but grouped by the last six flips, and counts for how often that group has appeared. Then do the 100 flips one at a time, updating the groups and their counts. Any group that at some point ends with a streak doesn't continue playing, instead I add it to the total result. The group occurred count times, there are flips - 1 - i flips left, and they can be anything, so multiply count with 2^{flips - 1 - i}.

Results (note that the result for 20 flips is the same as with the naive solution):

>>> fast(20, 6)
248384
>>> fast(100, 6)
1022766552856718355261682015984

And dividing by 2¹⁰⁰ gives me the percentage similar to those of the linked-to experiments:

>>> 100 * fast(100, 6) / 2**100
80.68205487163246

Any comments, suggestions for improvement?

Answer 1

Your code looks good. It's a little hard to read, but given the context that is ok! We can also see that if new_ending is never in needles then your code looks like it will run in \$O(f2^s)\$ time, where \$f\$ is flips and \$s\$ is streak.

Whilst I can see the code in if new_ending in needles: will reduce the the time your code takes to run. For example when streak=2 it will allow your code to run in linear time, it's not going to help much on bigger numbers - the code will still tend to \$O(f2^s)\$.

We can see how you're performing this optimization in the following. Since you are not searching the descendent of HH, TT, HTT, THH, etc. it cuts down how big the tree will get.

We can clearly see that tails is just the inverse of heads. If we focus on heads and split the 'base' and 'tail' (the repeating results) we get the following:

     HH 1/2^2
H    TT 1/2^3
HT   HH 1/2^4
HTH  TT 1/2^5
HTHT HH 1/2^6

Whilst it's cool it runs in linear time, it's not really that interesting. And so when streak=2 the total chance for \$f\$ flips is:

$$\Sigma_{n=2}^f \frac{2}{2^n}$$

However when we look at streak=3 we can see the start of a distinguishing pattern.

     HHH 1/2^3
H    TTT 1/2^4
HH   TTT 1/2^5
HT   HHH 1/2^5
HHT  HHH 1/2^6
HTH  TTT 1/2^6
HTT  HHH 1/2^6
HHTH TTT 1/2^7
HHTT HHH 1/2^7
HTHH TTT 1/2^7
HTHT HHH 1/2^7
HTTH TTT 1/2^7

If we take the count of each size then we get:

3: 1
4: 1
5: 2
6: 3
7: 5

This is cool because it's the start of the Fibonacci numbers. I have verified that the first 30 values are the same. And so we now can assume we have an equation for streak=3:

$$\Sigma_{n=3}^f \frac{2F(n-2)}{2^n}$$

Doing the same thing for streak=4,5,6,10 give the following sequences:

• 4 - Tribonacci
• 5 - Tetranacci
• 6 - Pentanacci
• 10 - Fibonacci 9-step

In all this is a pretty compelling pattern. And so we can write an algorithm to run in \$O(fs)\$ time where \$f\$ is flips and \$s\$ is streaks.

import collections
import itertools
from fractions import Fraction


def fibonacci_nth(size):
    store = collections.deque([0] * size, size)
    store.append(1)
    while True:
        yield store[-1]
        store.append(sum(store))


def coin_chance(flips, streak):
    if streak <= 0 or streak % 1:
        raise ValueError("streak must be a positive integer")
    if flips < 0 or flips % 1:
        raise ValueError("flips must be a non-negative integer")
    if streak == 1:
        return Fraction(flips != 0, 1)
    sequence = (
        Fraction(2 * numerator, 2 ** exponent)
        for exponent, numerator in enumerate(fibonacci_nth(streak - 1), streak)
    )
    return sum(itertools.islice(sequence, flips - streak + 1))


# Code to get OEIS sequences
def funky_finder(depth, size):
    desired = (['H'] * size, ['T'] * size)
    stack = [iter("HT")]
    stack_value = []
    while stack:
        try:
            coin = next(stack[-1])
        except StopIteration:
            stack.pop()
            if stack_value:
                stack_value.pop()
            continue
        _stack_value = stack_value + [coin]
        if _stack_value[-size:] in desired:
            yield ''.join(_stack_value)
        elif len(stack) < depth:
            stack_value.append(coin)
            stack.append(iter('HT'))


# I know, I know. But I was using this in a REPL!
size = 3; [i // 2 for i in sorted(collections.Counter(len(i) - size for i in funky_finder(20 + size, size)).values())]

>>> 100 * fast(20, 6) / 2**20
23.687744140625
>>> 100 * float(coin_chance(20, 6))
23.687744140625

>>> 100 * fast(100, 6) / 2**100
80.68205487163246
>>> 100 * float(coin_chance(100, 6))
80.68205487163246