Symuluj zmienną Bernoulliego z prawdopodobieństwem

Czy ktoś może mi powiedzieć, jak symulować , gdzie , za pomocą rzutu monetą (tyle razy, ile potrzebujesz) z ?$\mathrm{Bernoulli}\left({a\over b}\right)$$a,b\in \mathbb{N}$$P(H)=p$

Myślałem o użyciu próbkowania odrzucenia, ale nie mogłem tego dopracować.

Ponieważ istnieje niezliczona ilość rozwiązań, znajdźmy skuteczne .

Idea tego zaczyna się od standardowego sposobu implementacji zmiennej Bernoulliego: porównaj jednolitą zmienną losową $U$ do parametru $a/b$. Kiedy$U \lt a/b$, powrót $1$; w przeciwnym razie wróć$0$.

Możemy użyć $p$-coin jako jednolity generator liczb losowych . Aby wygenerować liczbę$U$ równomiernie w dowolnym przedziale czasowym $[x, y)$, rzuć monetą. Kiedy jest to głowa, rekurencyjnie generuj jednolitą wartość$X$ na początku $p$część przedziału; gdy są to ogony, generuj rekurencyjnie$X$ od ostatniego $1-p$część przedziału. W pewnym momencie przedział docelowy stanie się tak mały, że tak naprawdę nie ma znaczenia, jak wybierzesz z niego liczbę: tak zaczyna się rekurencja. Oczywiste jest, że ta procedura generuje jednolite zmienne (do dowolnej pożądanej precyzji), co łatwo udowodnić indukcyjnie.

Ten pomysł nie jest skuteczny, ale prowadzi do wydajnej metody. Ponieważ na każdym etapie będziesz losować liczbę z określonego przedziału$[x,y)$, dlaczego najpierw nie sprawdzić, czy trzeba go narysować? Jeśli wartość docelowa leży poza tym przedziałem, znasz już wynik porównania między wartością losową a wartością docelową. Dlatego algorytm ten szybko się kończy. (Można to interpretować jako procedurę próbkowania odrzucenia wymaganą w pytaniu.)

Możemy dalej optymalizować ten algorytm. Na każdym etapie faktycznie mamy dwie monety, których możemy użyć: poprzez zmianę etykiety naszej monety możemy przekształcić ją w monetę, która jest szansą$1-p$. Dlatego jako wstępne obliczenie możemy rekurencyjnie wybrać, które ponowne etykietowanie prowadzi do niższej oczekiwanej liczby przerzutów potrzebnych do zakończenia. (To obliczenie może być kosztownym krokiem.)

Na przykład nieefektywne jest używanie monety $p=0.9$ naśladować Bernoulliego$(0.01)$zmienna bezpośrednio: średnio zajmuje prawie dziesięć rzutów. Ale jeśli użyjemy$p=1-0.0=0.1$ monetą, to już po dwóch rzutach na pewno się zakończymy, a oczekiwana liczba rzutów jest po prostu $1.2$.

Oto szczegóły.

Podziel dowolny podany półotwarty przedział $I = [x, y)$ w odstępach

To definiuje dwie transformacje $s(*,H)$ i $s(*,T)$ które działają w okresach półotwartych.

Jeśli chodzi o terminologię, jeśli $I$ to dowolny zestaw liczb rzeczywistych, niech wyrażenie

znaczy że $t$ jest dolną granicą dla $I$: $t \lt x$ dla wszystkich $x \in I$. Podobnie,$t \gt I$ znaczy $t$ jest górną granicą dla $I$.

pisać $a/b = t$. (W rzeczywistości nie ma znaczenia, czy$t$jest rzeczywisty zamiast racjonalnego; wymagamy tylko tego$0 \le t \le 1$.)

Oto algorytm tworzenia wariacji $Z$ z żądanym parametrem Bernoulli:

1. Zestaw $n=0$ i $I_n = I_0 = [0,1)$.

2. Podczas $(t\in I_{n})$ {Rzuć monetą, aby wyprodukować $X_{n+1}$. Zestaw$I_{n+1} = S(I_n, X_{n+1}).$ Przyrost $n$.}

3. Gdyby $t \gt I_{n+1}$ następnie ustaw $Z=1$. W przeciwnym razie ustaw$Z=0$.

Realizacja

Aby to zilustrować, oto Rimplementacja alorithm jako funkcji draw. Jego argumenty są wartością docelową$t$ i interwał $[x,y)$, początkowo $[0,1)$. Wykorzystuje simplementację funkcji pomocniczej$s$. Chociaż nie musi, śledzi także liczbę rzutów monetą. Zwraca losową zmienną, liczbę rzutów i ostatni sprawdzony interwał.

s <- function(x, ab, p) {
  d <- diff(ab) * p
  if (x == 1) c(ab[1], ab[1] + d) else c(ab[1] + d, ab[2])
}
draw <- function(target, p) {
  between <- function(z, ab) prod(z - ab) <= 0
  ab <- c(0,1)
  n <- 0
  while(between(target, ab)) {
    n <- n+1; ab <- s(runif(1) < p, ab, p)
  }
  return(c(target > ab[2], n, ab))
}

Jako przykład jego zastosowania i sprawdzenia jego dokładności weźmy przykład $t=1/100$ i $p=0.9$. Porysujmy$10,000$ wartości przy użyciu algorytmu, zgłaszają średnią (i jej błąd standardowy) oraz wskazują średnią liczbę zastosowanych przerzutów.

target <- 0.01
p <- 0.9
set.seed(17)
sim <- replicate(1e4, draw(target, p))

(m <- mean(sim[1, ]))                           # The mean
(m - target) / (sd(sim[1, ]) / sqrt(ncol(sim))) # A Z-score to compare to `target`
mean(sim[2, ])                                  # Average number of flips

W tej symulacji $0.0095$klapki były główkami. Chociaż niższy niż cel$0.01$, wynik Z na poziomie $-0.5154$nie jest znaczący: to odchylenie można przypisać przypadkowi. Średnia liczba przerzutów wyniosła$9.886$- trochę mniej niż dziesięć. Gdybyśmy skorzystali z$1-p$ monety, średnia byłaby $0.0094$- nadal nie różni się znacząco od celu, ale tylko $1.177$ przerzuty byłyby potrzebne średnio.

Oto rozwiązanie (trochę niechlujne, ale to moje pierwsze dźgnięcie). Możesz faktycznie zignorować$P(H) = p$ i WLOG zakłada $P(H)=1/2$. Dlaczego? Istnieje sprytny algorytm do generowania obiektywnego rzutu monetą z dwóch stronniczych rzutów monetą. Możemy więc założyć$P(H)=1/2$.

Aby wygenerować $\text{Bernoulli}(\frac{a}{b})$, Mogę wymyślić dwa rozwiązania (pierwsze nie jest moje, ale drugie jest uogólnieniem):

Rozwiązanie 1

Odwróć obiektywną monetę $b$czasy. Gdyby$a$głowy nie są obecne, zacznij od nowa. Gdyby$a$głowice są obecne, zwróć, czy pierwsza moneta jest głowicą, czy nie (ponieważ$P(\text{first coin is heads | $a$ heads in $b$ coins}) = \frac{a}{b}$)

Rozwiązanie 2

Można to rozszerzyć na dowolną wartość $\text{Bernoulli}(p)$. pisać$p$w formie binarnej. Na przykład,$0.1 = 0.0001100110011001100110011... \text{base 2}$

Stworzymy nowy numer binarny za pomocą rzutów monetą. Zacząć od$0.$i dodawaj cyfry w zależności od tego, czy pojawią się głowice (1) czy ogony (0). Przy każdym odwróceniu porównaj swój nowy numer binarny z reprezentacją binarną$p$ do tej samej cyfry . W końcu obie się rozejdą i wrócą, jeśli$bin(p)$ jest większy niż liczba binarna.

W Pythonie:

def simulate(p):
    binary_p = float_to_binary(p)
    binary_string = '0.'
    index = 3
    while True:
        binary_string += '0' if random.random() < 0.5 else '1'
        if binary_string != binary_p[:index]:
            return binary_string < binary_p[:index]
        index += 1

Niektóre dowody:

np.mean([simulate(0.4) for i in range(10000)])

wynosi około 0,4 (jednak nie szybko)

Widzę proste rozwiązanie, ale bez wątpienia istnieje wiele sposobów na zrobienie tego, niektóre prawdopodobnie łatwiejsze niż to. To podejście można podzielić na dwa etapy:

1. Generowanie z dwóch zdarzeń z jednakowym prawdopodobieństwem przy nieuczciwej procedurze losowania monet (połączenie konkretnej monety i metody jej rzucania generuje głowę z prawdopodobieństwem $p$). Możemy nazwać te dwa równie prawdopodobne zdarzenia$H^*$, i $T^*$. [Jest na to proste podejście, które wymaga pary rzutów$H^*=(H,T)$ i $T^*=(T,H)$ aby uzyskać dwa równie prawdopodobne wyniki, przy czym wszystkie inne wyniki prowadzą do wygenerowania nowej pary rzutów, aby spróbować ponownie.]

2. Teraz generujesz losowy spacer z dwoma stanami pochłaniania za pomocą symulowanej uczciwej monety. Wybierając odległość stanów pochłaniania od źródła (jeden powyżej i jeden poniżej niego), możesz ustawić szansę absorpcji, mówiąc, że górny stan pochłaniania jest pożądanym stosunkiem liczb całkowitych. W szczególności, jeśli umieścisz górną barierę pochłaniającą w$a$ i niższy o $-(b-a)$ (i rozpocznij proces od początku) i uruchom losowy spacer aż do absorpcji, prawdopodobieństwo absorpcji przy górnej barierze wynosi $\frac{a}{a+(b-a)} = \frac{a}{b}$.

   (Aby to pokazać, należy wykonać pewne obliczenia, ale można dość łatwo wyliczyć prawdopodobieństwa, pracując z relacjami powtarzalności ... lub możesz to zrobić, sumując nieskończone szeregi ... lub istnieją inne sposoby.)