Jak mogę modelować flipy, dopóki N nie odniesie sukcesu?

17

Ty i ja decydujemy się zagrać w grę, w której na zmianę podrzucamy monetę. Pierwszy gracz, który rzuci łącznie 10 głów, wygrywa. Oczywiście istnieje spór o to, kto powinien iść pierwszy.

Symulacje tej gry pokazują, że gracz, który przerzuca pierwszy, wygrywa o 6% więcej niż gracz, który przerzuca drugi (pierwszy gracz wygrywa przez około 53% czasu). Jestem zainteresowany modelowaniem tego w sposób analityczny.

To nie jest losowa zmienna dwumianowa, ponieważ nie ma ustalonej liczby prób (odwróć, aż ktoś zdobędzie 10 głów). Jak mogę to wymodelować? Czy to ujemny rozkład dwumianowy?


Aby móc odtworzyć moje wyniki, oto mój kod python:

import numpy as np
from numba import jit


@jit
def sim(N):

    P1_wins = 0
    P2_wins = 0

    for i in range(N):

        P1_heads = 0
        P2_heads = 0
        while True:

            P1_heads += np.random.randint(0,2)

            if P1_heads == 10:
                P1_wins+=1
                break

            P2_heads+= np.random.randint(0,2)
            if P2_heads==10:
                P2_wins+=1
                break
    return P1_wins/N, P2_wins/N


a,b = sim(1000000)
Demetri Pananos
źródło
3
Kiedy rzucasz monetą, aż r , a następnie patrzysz na rozkład liczby sukcesów, które zdarzają się przed zakończeniem takiego eksperymentu, to z definicji jest to rozkład ujemny dwumianowy .
Tim
2
Nie mogę odtworzyć wartości 2%. Uważam, że pierwszy gracz wygrywa 53.290977425133892% czasu.
whuber
1
@ whuber tak, uważam, że masz rację. Przeprowadziłem symulację mniej razy niż powinienem. Moje wyniki są proporcjonalne do twoich.
Demetri Pananos
1
Jeśli jeden wygrywa 53% czasu, drugi powinien wynosić 47%, więc czy opis nie powinien brzmieć „pierwszy gracz wygrywa 6% więcej niż drugi gracz” czy „3% więcej niż połowa czasu”? Nie (jak obecnie mówi) „3% więcej niż gracz, który przerzuca drugie miejsce”
JesseM
3
Czy otrzymałeś to pytanie od FiveThirtyEight Riddler Express ?
foutandabout

Odpowiedzi:

19

Rozkład liczby ogony przed osiągnięciem głowic ujemna dwumianowego parametrów 10 i 1 / 2 . Niech f będzie funkcją prawdopodobieństwa, a G funkcją przetrwania: dla każdego n 0 , f ( n ) oznacza szansę gracza na n- ogona przed 10 główami, a G ( n ) oznacza szansę gracza na n lub więcej ogonów przed 10 głowami.10101/2fGn0f(n)n10G(n)n10

Ponieważ gracze rzucają niezależnie, szansę, którą pierwszy gracz wygrywa, rzucając dokładnie ogonami, uzyskuje się, mnożąc tę ​​szansę przez szansę, że drugi gracz rzuca n lub więcej ogonami, równą f ( n ) G ( n ) .nnf(n)G(n)

Zsumowanie wszystkich możliwych daje szanse wygranej pierwszego gracza jakon

n=0f(n)G(n)53.290977425133892%.

To około więcej niż połowa czasu.3%

Zasadniczo, zastępując dowolną dodatnią liczbą całkowitą m , odpowiedź można podać w kategoriach funkcji hipergeometrycznej: jest równa10m

1/2+22m12F1(m,m,1,1/4).

Gdy używa się stronniczej monety z szansą głów, uogólnia się nap

12+12(p2m)2F1(m,m,1,(1p)2).

Oto Rsymulacja miliona takich gier. Podaje szacunkową wartość . Dwumianowy test hipotezy, aby porównać go z wynikiem teoretycznym, ma wynik Z wynoszący - 0,843 , co jest nieznaczną różnicą.0.53250.843

n.sim <- 1e6
set.seed(17)
xy <- matrix(rnbinom(2*n.sim, 10, 1/2), nrow=2)
p <- mean(xy[1,] <= xy[2,])
cat("Estimate:", signif(p, 4), 
    "Z-score:", signif((p - 0.532909774) / sqrt(p*(1-p)) * sqrt(n.sim), 3))
Whuber
źródło
1
Podobnie jak uwaga, która na pierwszy rzut oka może nie być oczywista, nasze odpowiedzi zgadzają się liczbowo: (.53290977425133892 - .5) * 2 to w zasadzie dokładnie takie prawdopodobieństwo, jakie podałem.
Dougal,
1
@Dougal Dziękujemy za zwrócenie na to uwagi. Spojrzałem na twoją odpowiedź, zobaczyłem i wiedząc, że nie zgadza się z formą odpowiedzi wymaganą w pytaniu, nie rozpoznałem, że poprawnie obliczyłeś. Ogólnie, jeśli to możliwe, dobrze jest ułożyć odpowiedź na dowolne pytanie w żądanym formularzu: ułatwia to rozpoznanie, kiedy jest poprawne i łatwe do porównania odpowiedzi. 6.6%
whuber
1
@ whuber Odpowiedziałem na wyrażenie „Symulacje tej gry pokazują, że gracz, który przerzucił pierwszy, wygrywa 2% (EDYCJA: 3% więcej po symulacji większej liczby gier) więcej niż gracz, który przerzucił drugi”. Zinterpretowałbym „wygrywa 2% więcej” jako ; prawidłowa wartość to rzeczywiście 6,6%. Nie jestem pewien, jak interpretować „wygrywa 2% więcej” oznacza „wygrywa 52% czasu”, choć najwyraźniej tak właśnie było. Pr(A wins)Pr(B wins)=2%
Dougal,
@Dougal Zgadzam się, że opis PO jest mylący, a nawet błędny. Jednak kod i jego wynik jasno wskazywały, że miał na myśli „3% więcej niż połowę czasu” zamiast „3% więcej niż inny gracz”.
whuber
1
@whuber Zgoda. Niestety odpowiedziałem na to pytanie przed opublikowaniem kodu i sam nie przeprowadziłem symulacji. :)
Dougal,
15

Możemy wymodelować grę w następujący sposób:

  • Gracz A wielokrotnie rzuca monetą, uzyskując wyniki A1,A2, aż do uzyskania łącznie 10 głów. Niech indeks czasu z szefów 10th być zmienna losowa X .
  • Gracz B robi to samo. Niech indeks czasu z szefów 10th być zmienna losowa Y , który jest iid kopia X .
  • Jeśli XY , Gracz A wygrywa; w przeciwnym razie Gracz B wygrywa. Oznacza to, że
    Pr(A wins)=Pr(XY)=Pr(X>Y)+Pr(X=Y)Pr(B wins)=Pr(Y>X)=Pr(X>Y).

Różnica w stawkach wygranych wynosi zatem

Pr(X=Y)=kPr(X=k,Y=k)=kPr(X=k)2.

Jak podejrzewasz, X (i Y ) są rozmieszczone zasadniczo zgodnie z ujemnym rozkładem dwumianowym. Oznaczenia tego są różne, ale w parametryzacji Wikipedii mamy głowy jako „porażkę”, a ogony jako „sukces”; potrzebujemy r=10 „awarii” (głów), zanim eksperyment zostanie zatrzymany, a prawdopodobieństwo sukcesu p=12 . Wtedy liczba „sukcesów”, która wynosiX10, ma

Pr(X10=k)=(k+9k)210k,
a prawdopodobieństwo zderzenia wynosi
Pr(X=Y)=k=0(k+9k)222k20,
co Mathematica mówi nam, że ma7649952511622614676.6%.

Tak więc wskaźnik wygranych Gracza B wynosi Pr(Y>X)46.7% , a Gracza A wynosi 619380496116226146753.3%

Dougal
źródło
głowy nie muszą znajdować się w rzędzie, a jedynie 10. Zakładam, że to naprawiasz.
Demetri Pananos
6
(+1) Podobało mi się to podejście bardziej niż to, które opublikowałem, ponieważ jest ono obliczeniowo prostsze: wymaga tylko funkcji prawdopodobieństwa, która ma proste wyrażenie w kategoriach współczynników dwumianowych.
whuber
1
I've submitted an edit replacing the last paragraph questioning the difference from the other answer with an explanation of how their results are actually the same.
Monty Harder
1

Let Eij be the event that the player on roll flips i heads before the other player flips j heads, and let X be the first two flips having sample space {hh,ht,th,tt} where h means heads and t tails, and let pijPr(Eij).

Then pij=Pr(Ei1j1|X=hh)Pr(X=hh)+Pr(Ei1j|X=ht)Pr(X=ht)+Pr(Eij1|X=th)Pr(X=th)+Pr(Eij|X=tt)Pr(X=tt)

Assuming a standard coin Pr(X=)=1/4 means that pij=1/4[pi1j1+pi1j+pij1+pij]

solving for pij, =1/3[pi1j1+pi1j+pij1]

But p0j=p00=1 and pi0=0, implying that the recursion fully terminates. However, a direct naive recursive implementation will yield poor performance because the branches intersect.

An efficient implementation will have complexity O(ij) and memory complexity O(min(i,j)). Here's a simple fold implemented in Haskell:

Prelude> let p i j = last. head. drop j $ iterate ((1:).(f 1)) start where
  start = 1 : replicate i 0;
  f c v = case v of (a:[]) -> [];
                    (a:b:rest) -> sum : f sum (b:rest) where
                     sum = (a+b+c)/3 
Prelude> p 0 0
1.0
Prelude> p 1 0
0.0
Prelude> p 10 10
0.5329097742513388
Prelude> 

UPDATE: Someone in the comments above asked whether one was suppose to roll 10 heads in a row or not. So let Ekl be the event that the player on roll flips i heads in a row before the other player flips i heads in a row, given that they already flipped k and l consecutive heads respectively.

Proceeding as before above, but this time conditioning on the first flip only, pk,l=11/2[pl,k+1+pl,0] where pil=pii=1,pki=0

This is a linear system with i2 unknowns and one unique solution.

To convert it into an iterative scheme, simply add an iterate number n and a sensitivity factor ϵ:

pk,l,n+1=1/(1+ϵ)[ϵpk,l,n+11/2(pl,k+1,n+pl,0,n)]

Choose ϵ and pk,l,0 wisely and run the iteration for a few steps and monitor the correction term.

John Rambo
źródło