Как я могу моделировать сальто до N успехов?

Мы с тобой решили сыграть в игру, в которой мы по очереди подбрасываем монету. Первый игрок, который перевернет 10 голов, выигрывает игру. Естественно, есть спор о том, кто должен идти первым.

Моделирование этой игры показывает, что игрок, который переворачивает первый, выигрывает на 6% больше, чем игрок, который переворачивает второй (первый игрок выигрывает примерно в 53% времени). Я заинтересован в моделировании этого аналитически.

Это не биноминальная случайная величина, так как не существует фиксированного количества испытаний (переворачивайте, пока кто-то не получит 10 голов). Как я могу смоделировать это? Это отрицательное биномиальное распределение?

Чтобы восстановить мои результаты, вот мой код на 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)

Распределение числа хвостов до достижения головок является отрицательным биномиальным с параметрами 10 и 1 / 2 . Пусть F есть функция вероятности и G функция выживания: для каждого п ≥ 0 , ф ( п ) есть шанс игрока из п хвостов перед тем$10$$10$$1/2$$f$$G$$n\ge 0$$f(n)$$n$ голов и G ( п ) есть шанс игрока из русских или более хвостовпрежде чем 10 голов.$10$$G(n)$$n$$10$

Поскольку игроки бросают независимо, шанс того, что первый игрок выиграет с броском ровно $n$ хвостов, получается умножением этого шанса на шанс того, что второй игрок бросит или более хвостов, равный f ( n ) G ( n ) .$n$$f(n)G(n)$

Суммирование по всем возможным дает шансы на победу первого игрока как$n$

Это примерно на больше, чем в половине случаев.$3\%$

В общем случае, заменяя на любое положительное целое число m , ответ можно дать в терминах гипергеометрической функции: он равен$10$$m$

При использовании смещенной монеты с вероятностью голов, это обобщает$p$

Вот Rсимуляция миллиона таких игр. Он сообщает, что оценка $0.5325$ . Тест биномиальной гипотезы сравнить его с теоретическим результатом имеет Z-балл , что является незначительным различием.$-0.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))

Мы можем смоделировать игру так:

• Игрок A многократно подбрасывает монету, получая результаты $A_1, A_2, \dots$ пока они не получат в общей сложности 10 голов. Пусть временной индекс 10-й головы будет случайной величиной$X$ .
• Игрок Б делает то же самое. Пусть индекс времени 10-й головы будет случайной величиной $Y$ , которая является iid-копией$X$ .
• Если $X \le Y$ , игрок A побеждает; в противном случае игрок B выигрывает. То есть $$\begin{align} \Pr(A\text{ wins})&= \Pr(X \ge Y) = \Pr(X > Y) + \Pr(X = Y)\\ \Pr(B\text{ wins})&= \Pr(Y > X) = \Pr(X > Y). \end{align}$$

Разрыв в показателях выигрыша, таким образом ,

Как вы подозревали, $X$ (и $Y$ ) распределены по существу в соответствии с отрицательным биномиальным распределением. Обозначения для этого различаются, но в параметризации Википедии мы имеем головы как «провал» и хвосты как «успех»; нам нужно $r = 10$ «неудач» (голов) до того, как эксперимент будет остановлен, а вероятность успеха $p = \tfrac12$ . Тогда число «успехов», которое составляет$X - 10$, имеет

Таким образом, выигрыш игрока B составляет $\Pr(Y > X) \approx 46.7\%$ , а игрока A - $\frac{619\,380\,496}{1\,162\,261\,467} \approx 53.3\%$.

Пусть $E_{ij}$ будет событием, когда игрок при броске переворачивает мои головы до того, как другой игрок перевернет j голов, и пусть $X$ будет первыми двумя переворотами, имеющими пробное пространство $\{ hh,ht,th,tt\}$ где h означает головы и т хвосты, и пусть $p_{ij} \equiv Pr(E_{ij})$ .

Тогда $p_{ij}=Pr(E_{i-1j-1}|X=hh)*Pr(X=hh)+Pr(E_{i-1j}|X=ht)*Pr(X=ht)+Pr(E_{ij-1}|X=th)*Pr(X=th)+Pr(E_{ij}|X=tt)*Pr(X=tt)$

Assuming a standard coin $Pr(X=*)=1/4$ means that $p_{ij}=1/4*[p_{i-1j-1}+p_{i-1j}+p_{ij-1}+p_{ij}]$

solving for $p_{ij}$, $= 1/3*[p_{i-1j-1}+p_{i-1j}+p_{ij-1}]$

But $p_{0j}=p_{00}=1$ and $p_{i0}=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(i*j)$ 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 $E_{kl}$ 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, $p_{k,l} = 1-1/2*[p_{l,k+1}+p_{l,0}]$ where $p_{il}=p_{ii}=1, p_{ki}=0$

This is a linear system with $i^2$ unknowns and one unique solution.

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

$p_{k,l,n+1} = 1/(1+\epsilon)*[\epsilon*p_{k,l,n} +1-1/2*(p_{l,k+1,n}+p_{l,0,n})]$

Choose $\epsilon$ and $p_{k,l,0}$ wisely and run the iteration for a few steps and monitor the correction term.