4
\$\begingroup\$

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 2flips - 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 2100 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?

\$\endgroup\$
2
  • \$\begingroup\$ What is the aim, to make it as fast as possible? Is just using math an option? \$\endgroup\$
    – Peilonrayz
    Aug 21, 2020 at 22:25
  • 1
    \$\begingroup\$ @Peilonrayz Yes, a math solution sounds very interesting. I'm not entirely sure what I'm looking for, as I believe it's already fairly good. But I've seen reviews with things that hadn't occurred to me, so maybe someone can surprise me here as well. Maybe style comments, as I've had people complain that my solutions are unreadable :-). It's kind of an experiment. \$\endgroup\$ Aug 21, 2020 at 22:52

1 Answer 1

3
\$\begingroup\$

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.

enter image description here

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:

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
\$\endgroup\$
5
  • \$\begingroup\$ Will have a closer look later, looks intriguing! Question already: Why \$O(2^f)\$? My groups never grows larger than \$2^{streak}\$, so I'd say \$O(2^{streak}f)\$ operations. And \$O(2^{streak}f^2)\$ time, as the operations are adding f-bit numbers. \$\endgroup\$ Aug 22, 2020 at 2:25
  • \$\begingroup\$ @superbrain Yes, on second look your code looks like it may not be \$O(2^f)\$. Could possibly be, as you say, \$O(f2^s)\$ however I'm too tired to verify that. I'll edit my answer to address this. \$\endgroup\$
    – Peilonrayz
    Aug 22, 2020 at 2:31
  • \$\begingroup\$ To clarify: I wouldn't call my if new_ending in needles: an optimization. It's a simplification :-). Or rather a non-complification. With s=6, only two of 64 endings don't continue playing. You're right, that doesn't matter speed-wise. That was never the point. The point was that I didn't want to additionally keep track of whether a streak has been found already. Though I've done that now and it was easier than I had guessed. \$\endgroup\$ Aug 22, 2020 at 13:05
  • \$\begingroup\$ Cool observation about the 'naccis. Proof would be nice, if not too complicated, but I'm convinced enough :-). Using deque like that is neat. Slight nitpick: streak=0 should lead to probability 100%. My naive solution does that. Though I admit it's better to reject it (like you do) than to fail it (like I just realized my fast solution does). \$\endgroup\$ Aug 22, 2020 at 13:40
  • \$\begingroup\$ Could also init the deque with just [1]. \$\endgroup\$ Aug 22, 2020 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.