Probabilità esatta per le serie di coin flip

Question

In un'altra domanda l'utente stava determinando la probabilità di avere una serie di 6 teste o 6 croci in 100 lanci di monete. Per trovare la probabilità, genererebbero 100 lanci di monete casuali e determinerebbero se ci fosse una serie. Avrebbero testato 10.000 di tali sequenze di 100 lanci per scoprire che c'era circa l'80% di possibilità che ci fosse una serie di 100 lanci di monete.

Ho deciso di calcolare la probabilità esatta. Per 100 lanci ci sono \$2^{100}\$possibili esiti. Per determinare la percentuale calcolo quanti di loro hanno una striscia e poi divido per \$2^{100}\$.

La mia soluzione ingenua mi fa ottenere il numero per 20 lanci in pochi secondi:

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))

Risultato:

>>> naive(20, 6)
248384

La mia soluzione veloce mi fa ottenere immediatamente il numero per 100 lanci:

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

L'idea è di avere un pool di giochi ancora in corso, ma raggruppati in base agli ultimi sei lanci, e conta la frequenza con cui quel gruppo è apparso. Quindi esegui i 100 lanci uno alla volta, aggiornando i gruppi e i loro conteggi. Qualsiasi gruppo che a un certo punto finisce con una striscia non continua a giocare, invece lo aggiungo al risultato totale. Il gruppo si è verificato countvolte, ci sono flips - 1 - ilanci rimasti e possono essere qualsiasi cosa, quindi moltiplica countcon 2 ^{lanci - 1 - i} .

Risultati (si noti che il risultato per 20 lanci è lo stesso della soluzione ingenua):

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

E dividendo per 2 ¹⁰⁰ mi dà la percentuale simile a quella degli esperimenti collegati:

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

Eventuali commenti, suggerimenti per il miglioramento?

Answer 1

Il tuo codice sembra buono. È un po' difficile da leggere, ma visto il contesto va bene! Possiamo anche vedere che se new_endingnon è mai presente needles, il tuo codice sembra che verrà eseguito in \$O(f2^s)\$ora, dove \$f\$è flipse \$s\$è streak.

Mentre posso vedere il codice if new_ending in needles:ridurrà il tempo necessario per l'esecuzione del codice. Ad esempio, quando streak=2 consentirà al tuo codice di essere eseguito in tempo lineare, non sarà di grande aiuto su numeri più grandi: il codice tenderà comunque a \$O(f2^s)\$.

Possiamo vedere come stai eseguendo questa ottimizzazione di seguito. Dal momento che non stai cercando il discendente di HH, TT, HTT, THH, ecc., riduce la dimensione dell'albero.

Possiamo vedere chiaramente che croce è solo l'inverso di testa. Se ci concentriamo sulle teste e dividiamo la "base" e la "coda" (i risultati ripetuti) otteniamo quanto segue:

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

Sebbene sia bello, funziona in tempo lineare, non è poi così interessante. E così quando la serie = 2 la possibilità totale per \$f\$ribalta è:

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

Tuttavia, quando osserviamo la sequenza = 3, possiamo vedere l'inizio di uno schema distintivo.

     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

Se prendiamo il conteggio di ogni dimensione, otteniamo:

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

Questo è bello perché è l'inizio dei numeri di Fibonacci . Ho verificato che i primi 30 valori sono gli stessi. E quindi ora possiamo supporre di avere un'equazione per streak=3:

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

Facendo la stessa cosa per streak=4,5,6,10 si ottengono le seguenti sequenze:

• 4 - Tribonaccio
• 5 - Tetranacci
• 6 - Pentanaccio
• 10 - Fibonacci 9 passi

In tutto questo c'è uno schema piuttosto avvincente. E così possiamo scrivere un algoritmo da eseguire in \$O(fs)\$ora in cui \$f\$è capovolge e \$s\$è striature.

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