Jak wygenerować niecałkowitą liczbę kolejnych sukcesów Bernoulliego?

18

Dany:

Moneta o nieznanym nastawieniu (Głowa). $p$
Ściśle dodatni realny . $a > 0$

Problem:

Wygeneruj losową zmienną Bernoulliego z odchyleniem . $p^{a}$

Czy ktoś wie jak to zrobić? Na przykład, gdy jest dodatnia, to można odwrócić monetę czasu i zobaczyć, czy wszystkie wyniki były Głowy: jeśli są one wtedy problem „0”, w przeciwnym razie problem „1”. Trudność polega na tym, że niekoniecznie jest liczbą całkowitą. Ponadto, gdybym znał stronniczość , mógłbym po prostu zbudować kolejną monetę z pożądanym uprzedzeniem. $a$ $a$ $a$ $p$

sampling

— Pedro A. Ortega
źródło

2

@gung: Myślę, że potrzebny jest algorytm do generowania wariantu Bernoulliego, biorąc pod uwagę monetę.

— Neil G

1

Myślę, że chodzi o to, gdy tylko zachować średnio 1 na każde głów, które pojawia się i gdy , powielać każdej z głowic średnio czasach.

a > 1

$a>1$

a

$a$

a < 1

$a < 1$

1 / a

$1/a$

— Makro

3

@Macro, czy mógłbyś rozwinąć ten pomysł?

— Pedro A. Ortega

1

Drogi Pedro, (+1) do twojego postu, co jest pytaniem, które sprawia, że CV jest bardzo orzeźwiające i stymulujące, przynajmniej dla mnie. Czy mogę zapytać, jakie jest źródło tego pytania?

— kardynał

@cardinal: Jeszcze raz dziękuję za odpowiedź! Ten problem jest częścią próbnika do rozwiązywania stochastycznych problemów kontrolnych, nad którymi pracuję. Powodem, dla którego jest nieznane, jest to, że wymagałoby znajomości stałej normalizującej (która w tym przypadku jest nieprzyjemną funkcją podziału), ale nadal możemy próbkować z niej przy użyciu próbkowania odrzucającego. Przy okazji fajnie byłoby zacytować cię po imieniu, a nie tylko link do CV ;-).

p

$p$

— Pedro A. Ortega

19

Możemy to rozwiązać za pomocą kilku „sztuczek” i odrobiny matematyki.

Oto podstawowy algorytm:

Wygeneruj losową zmienną geometryczną z prawdopodobieństwem sukcesu . $p$
Wynik tej zmiennej losowej określa stałą znaną wartość . $f_n \in [0,1]$
Wygeneruj losową zmienną używając uczciwych rzutów monetą wygenerowanych z blokowanych sparowanych rzutów naszej monety . $\mathrm{Ber}(f_n)$ $\mathrm{Ber}(p)$
Wynikowy wynik to dla dowolnego , co jest wszystkim, czego potrzebujemy. $\mathrm{Ber}(p^a)$ $a \in (0,1)$

Aby uczynić rzeczy bardziej strawnymi, podzielimy je na kawałki.

Część 1 : bez utraty ogólności założono, że . $0 < a < 1$

Jeśli , możemy napisać dla pewnej dodatniej liczby całkowitej a dla . Ale dla dwóch niezależnych Bernoulli mamy Możemy wygenerować Bernoulliego naszych monet w sposób oczywisty. Dlatego musimy się tylko zajmować generowaniem gdy . $a \geq 1$ $p^a = p^n p^b$ $n$ $0 \leq b < 1$

P (X_{1} = X_{2} = 1) = p_{1} p_{2} .

$\renewcommand{\Pr}{\mathbb P} \Pr(X_1 = X_2 = 1) = p_1 p_2 \>.$

p^{n}

$p^n$

B e r (p^{a})

$\mathrm{Ber}(p^a)$

a \in (0, 1)

$a \in (0,1)$

Część 2 : Dowiedz się, jak wygenerować dowolny z rzutu uczciwą monetą. $\mathrm{Ber}(q)$

Istnieje standardowy sposób, aby to zrobić. Rozwiń w rozwinięciu binarnym, a następnie użyj naszych uczciwych rzutów monet, aby „dopasować” cyfry . Pierwsze dopasowanie określa, czy zadeklarujemy sukces („główki”), czy porażkę („ogony”). Jeśli a nasze monetą są głowami, zadeklaruj głowy, jeśli a nasze monetą stanowią reszki, zadeklaruj reszki. W przeciwnym razie rozważ kolejną cyfrę w stosunku do nowego rzutu monetą. $q = 0.q_1 q_2 q_3 \ldots$ $q$ $q_n = 1$ $q_n = 0$

Część 3 : Dowiedz się, jak wygenerować uczciwą monetę od nieuczciwych z nieznanym uprzedzeniem.

Odbywa się to przy założeniu , poprzez przerzucanie monety parami. Jeśli otrzymamy , zadeklaruj głowy; jeśli otrzymamy , zadeklaruj ogony i w przeciwnym razie powtórz eksperyment, aż wystąpi jeden z dwóch wyżej wymienionych wyników. Są one równie prawdopodobne, więc musi mieć prawdopodobieństwo . $p \in (0,1)$ $HT$ $TH$ $1/2$

Część 4 : Trochę matematyki. (Taylor na ratunek.)

Rozwijając wokół , twierdzenie Taylora stwierdza, że Zauważ, że ponieważ każdy termin po pierwszym jest ujemny, więc mamy gdzie są znane z góry . Stąd gdzie , i $h(p) = p^a$ $p_0 = 1$

p^{a} = 1 - a (1 - p) - \frac{a (1 - a)}{2!} (1 - p)^{2} - \frac{a (1 - a) (2 - a)}{3!} (1 - p)^{3} \dots .

$p^a = 1 - a(1-p) - \frac{a(1-a)}{2!} (1-p)^2 - \frac{a(1-a)(2-a)}{3!} (1-p)^3 \cdots \>.$

0 < a < 1

$0 < a < 1$

p^{a} = 1 - \sum_{n = 1}^{\infty} b_{n} (1 - p)^{n},

$p^a = 1 - \sum_{n=1}^\infty b_n (1-p)^n \>,$

0 \leq b_{n} \leq 1

$0 \leq b_n \leq 1$

1 - p^{a} = \sum_{n = 1}^{\infty} b_{n} (1 - p)^{n} = \sum_{n = 1}^{\infty} b_{n} P (G \geq n) = \sum_{n = 1}^{\infty} f_{n} P (G = n) = E f (G),

$1 - p^a = \sum_{n=1}^{\infty} b_n (1-p)^n = \sum_{n=1}^\infty b_n \Pr(G \geq n) = \sum_{n=1}^\infty f_n \Pr(G = n) = \mathbb E f(G),$

G \sim G e o m (p)

$G \sim \mathrm{Geom}(p)$

f_{0} = 0

$f_0 = 0$

f_{n} = \sum_{k = 1}^{n} b_{k}

$f_n = \sum_{k=1}^n b_k$ dla .

n \geq 1

$n \geq 1$

I już wiemy, jak wykorzystać naszą monetę do wygenerowania losowej zmiennej geometrycznej z prawdopodobieństwem sukcesu . $p$

Część 5 : Sztuczka Monte Carlo.

Niech będzie dyskretną zmienną losową przyjmującą wartości w z . Niech . Następnie $X$ $[0,1]$ $\Pr(X = x_n) = p_n$ $U \mid X \sim \mathrm{Ber}(X)$

P (U = 1) = \sum_{n} x_{n} p_{n} .

$\Pr(U = 1) = \sum_n x_n p_n.$

Ale biorąc i , widzimy teraz, jak wygenerować zmienną losową i jest to równoważne z generowaniem jeden. $p_n = p(1-p)^n$ $x_n = f_n$ $\mathrm{Ber}(1-p^a)$ $\mathrm{Ber}(p^a)$

— kardynał
źródło

Jak mogę zacytować ciebie (lub twoje rozwiązanie)?

— Pedro A. Ortega,

2

@Pedro: Przypuszczam, że możesz kliknąć link „udostępnij” na dole tej odpowiedzi. Powinien to być stabilny link. Math.SE ma mechanizm cytowania , który nie wydaje się być włączony na tej stronie, ale możesz go dostosować.

— kardynał

1

Teraz, to jest genialny odpowiedź!

— Zen,

1

Napisałem to na forum dyskusji ogólnej klasy Coursera na temat kombinatoryki analitycznej, ponieważ było to miłe zastosowanie szeregów mocy związanych z niektórymi materiałami tam zawartymi. class.coursera.org/introACpartI-001/forum/thread?thread_id=108

— Douglas Zare

@Douglas: Dzięki! Czy istnieje publicznie dostępna wersja tego wątku, czy też muszę zapisać się na kurs, aby go zobaczyć? Pedro i ja omawialiśmy (za pośrednictwem poczty elektronicznej) możliwe sposoby włączenia tego podejścia do niektórych jego badań.

— kardynał

6

Czy poniższa odpowiedź jest głupia?

Jeśli są niezależne a ma rozkład , wtedy będzie w przybliżeniu dystrybuowane jako , kiedy . $X_1,\dots,X_n$ $\mathrm{Ber}(p)$ $Y_n$ $\mathrm{Ber}\left(\left(\sum_{i=1}^n X_i/n \right)^a\right)$ $Y_n$ $\mathrm{Ber}(p^a)$ $n\to\infty$

Dlatego jeśli nie znasz , ale możesz rzucić tę monetą wiele razy, możesz spróbować (w przybliżeniu) z losowej zmiennej . $p$ $\mathrm{Ber}(p^a)$

Przykładowy Rkod:

n <- 1000000
p <- 1/3 # works for any 0 <= p <= 1
a <- 4
x <- rbinom(n, 1, p)
y <- rbinom(n, 1, mean(x)^a)
cat("p^a =", p^a, "\n")
cat("est =", mean(y))

Wyniki:

p^a = 0.01234568 
est = 0.012291

— Zen
źródło

2

Podoba mi się ta odpowiedź, ale podejrzewam, że pomija sens pytania, które zinterpretowałem jako prośbę o algorytm, który generuje z żądanego rozkładu bez znajomości (lub informacji empirycznych o ). Ale problem zakłada, że możesz wygenerować zmienne losowe , więc jest to całkowicie rozsądna odpowiedź i wcale nie jest głupia! +1

p

$p$

p

$p$

B e r n o u l l i (p)

${\rm Bernoulli}(p)$

— Makro

1

+1: Podoba mi się. Zakładam, że masz na myśli, że jest rozpowszechniany…?

Y_{n}

$Y_n$

— Neil G

Dużo lepiej! Tks, @Neil G!

— Zen,

1

To urocze (+1), ale możemy to zrobić dokładnie w prawie na pewno skończonej liczbie przewrotów (i średnio ta liczba będzie względnie mała).

— kardynał

5

Poniższą ekspozycję tego pytania i odpowiedź kardynała opublikowałem na forum dyskusyjnym bieżącej klasy Analitycznej Kombinatorii na Coursera: „Zastosowanie szeregów mocy do konstruowania zmiennej losowej”. Publikuję tutaj kopię jako wiki społeczności, aby uczynić to publicznie i bardziej trwałym dostępnym.

Na stat.stackexchange.com pojawiło się interesujące pytanie i odpowiedź dotyczące serii potęg: „Jak wygenerować liczbę całkowitą kolejnych sukcesów Bernoulliego?” Sparafrazuję pytanie i odpowiedź kardynałem.

Załóżmy, że mamy potencjalnie niesprawiedliwą monetę, która jest głowami z prawdopodobieństwem , i dodatnią liczbą rzeczywistą . Jak możemy skonstruować zdarzenie, którego prawdopodobieństwem jest ? $p$ $\alpha$ $p^\alpha$

Gdyby była dodatnią liczbą całkowitą, moglibyśmy po prostu odwrócić monetę i pozwolić, aby zdarzenie polegało na tym, że wszystkie rzuty były główkami. Jeśli jednak nie jest liczbą całkowitą, powiedzmy , to nie ma to sensu, ale możemy użyć tego pomysłu, aby zredukować do przypadku, że . Jeśli chcemy skonstruować zdarzenie, którego prawdopodobieństwo wynosi , przyjmujemy przecięcie niezależnych zdarzeń, których prawdopodobieństwem jest i . $\alpha$ $\alpha$ $\alpha$ $1/2$ $0 \lt \alpha \lt 1$ $p^{3.5}$ $p^3$ $p^{0.5}$

Jedyne, co możemy zrobić, to skonstruować zdarzenie z dowolnym znanym prawdopodobieństwem . Aby to zrobić, możemy skonstruować strumień uczciwych bitów, wielokrotnie przewracając monetę, odczytując jako i jako oraz ignorując i . Porównujemy ten strumień z binarnym rozszerzeniem . Zdarzenie, w którym pierwsze nieporozumienie występuje, gdy ma prawdopodobieństwo . Nie znamy , więc nie możemy użyć tego bezpośrednio, ale będzie to przydatne narzędzie. $p' \in [0,1]$ $HT$ $1$ $TH$ $0$ $HH$ $TT$ $p' = 0.a_1a_2a_3..._2$ $a_i=1$ $p'$ $p^\alpha$

Główną ideą jest to, że chcielibyśmy użyć szeregu mocy dla gdzie . Możemy konstruować zdarzenia o prawdopodobieństwie , odwracając monetę razy i sprawdzając, czy wszystkie są ogonami, i możemy wygenerować zdarzenie z prawdopodobieństwem , porównując binarne cyfry ze sprawiedliwym strumieniem bitów jak wyżej i sprawdzenie, czy rzutów są ogony. $p^\alpha = (1-q)^\alpha = 1 - \alpha q - \frac{\alpha(1-\alpha)}{2} q^2 - \frac{\alpha (1-\alpha)(2-\alpha)}{3!}q^3 -...$ $p=1-q$ $q^n$ $n$ $p' q^n$ $p'$ $n$

Skonstruuj geometryczną zmienną losową z parametrem . Jest to liczba ogonów przed pierwszą głową w nieskończonej sekwencji rzutów monetą. . (Niektóre osoby używają definicji, która różni się o ). $G$ $p$ $P(G=n) = (1-p)^np = q^n p$ $1$

Biorąc pod uwagę kolejność , możemy produkować : Przerzuć monetę do pierwszego głowie, a jeśli istnieją ogony przed pierwszym głowy, podejmują element sekwencji indeksu . Jeżeli każdy , możemy porównać z jednolitą zmienną losową w (skonstruowanym jak powyżej), aby otrzymać zdarzenie z prawdopodobieństwem . $t_0, t_1, t_2, ...$ $t_G$ $G$ $G$ $t_n \in [0,1]$ $t_G$ $[0,1]$ $E[t_G] = \sum_n t_n P(G=n) = \sum_n t_n q^n p$

To jest prawie to, czego potrzebujemy. Chcielibyśmy wyeliminować to aby użyć szeregu mocy dla w . $p$ $p^\alpha$ $q$

1 = p + q p + q^{2} p + q^{3} p + . . .

$1 = p + qp + q^2p + q^3p + ...$

q^{n} = q^{n} p + q^{n + 1} p + q^{n + 2} p + . . .

$q^n = q^np + q^{n+1}p + q^{n+2}p + ...$

\begin{array}{rcl} \sum_{n} s_{n} q^{n} & = & \sum_{n} s_{n} (q^{n} p + q^{n + 1} p + q^{n + 2} p + . . .) \\ = & \sum_{n} (s_{0} + s_{1} + . . . + s_{n}) q^{n} p \end{array}

$\begin{eqnarray} \sum_n s_n q^n & = & \sum_n s_n (q^n p + q^{n+1}p + q^{n+2}p + ...) \newline & = & \sum_n (s_0 + s_1 + ... + s_n) q^n p \end{eqnarray}$

Rozważmy . Niech będzie sumą współczynników od do . Następnie . Każdy ponieważ współczynniki są dodatnie i sumują się do , więc możemy skonstruować zdarzenie z prawdopodobieństwem , porównując strumień bitów z ekspansją binarną z . Uzupełnienie ma prawdopodobieństwo zgodnie z wymaganiami. $1-p^\alpha = \alpha q + \frac{\alpha(1-\alpha)}{2} q^2 + ...$ $t_n$ $q$ $q^n$ $1-p^\alpha = \sum_n t_n q^n p$ $t_n\in [0,1]$ $1-0^\alpha = 1$ $1-p^\alpha$ $t_G$ $p^\alpha$

Ponownie argumentem jest kardynał.

— Douglas Zare
źródło

1

(+1) Dziękujemy za włożenie kłopotów z opublikowaniem tego. Różnice w ekspozycjach, choć stosunkowo niewielkie, pomagają wyjaśnić podejście.

— kardynał

4

Bardzo kompletna odpowiedź kardynała i późniejsze wypowiedzi zainspirowały następującą uwagę / wariant.

Niech PZ oznacza „Prawdopodobieństwo zera” i . Jeśli jest identyczną sekwencją Bernoulliego z PZ , to jest rv Bernoulliego z PZ . Teraz czyni losowym, tzn. Zastąpienie go liczbą całkowitą rv prowadzi do Bernoulli rv z Jeśli więc i jeśli weźmiemy z odpowiedzi kardynała , znajdziemy $q:=1-p$ $X_n$ $q$ $M_n := \max(X_1,\,X_2,\,\dots, X_n)$ $q^n$ $n$ $N \geq 1$ $M_N$

P r {M_{N} = 0} = \sum_{n = 1}^{\infty} P r {M_{N} = 0 | N = n} P r {N = n} = \sum_{n = 1}^{\infty} P r {N = n} q^{n} .

$\mathrm{Pr}\{M_N =0\} = \sum_{n=1}^\infty \mathrm{Pr}\{M_N =0 \,\vert\, N =n\} \mathrm{Pr}\{N =n\} = \sum_{n=1}^\infty \mathrm{Pr}\{N =n\} \, q^n.$

0 < a < 1

$0 < a < 1$

P r {N = n} = b_{n}

$\mathrm{Pr}\{N =n\} =b_n$

P r {M_{N} = 0} = 1 - p^{a}

$\mathrm{Pr}\{M_N =0\} = 1- p^a$ a to zgodnie z . Jest to rzeczywiście możliwe, ponieważ współczynniki spełniają i sumują się do .

1 - M_{N}

$1-M_N$

B e r (p^{a})

$\mathrm{Ber}(p^a)$

b_{n}

$b_n$

b_{n} ⩾ 0

$b_n \geqslant 0$

1

$1$

Dyskretny rozkład zależy tylko od z , wycofywanie Ma ciekawe funkcje. Okazuje się, że ma nieskończone oczekiwania i zachowanie ciężkiego ogona przy . $N$ $a$ $0 < a < 1$

P r {N = n} = \frac{a}{n} \prod_{k = 1}^{n - 1} (1 - a / k) (n \geq 1) .

$\mathrm{Pr}\{N =n\} = \frac{a}{n}\,\prod_{k=1}^{n-1}\left(1 - a/k\right) \qquad (n \geq 1).$

n b_{n} \sim c / n^{a}

$n \,b_n \sim c/n^a$

c = - 1 / Γ (- a) > 0

$c = -1/\Gamma(-a) >0$

Chociaż jest maksymalnym rvs, jego określenie wymaga liczby która jest ponieważ wynik jest znany, gdy tylko jeden wynosi . Liczba obliczonych jest rozkładem geometrycznym. $M_N$ $N$ $X_k$ $\leq N$ $X_k$ $1$ $X_k$

— Yves
źródło

pomysłem byłoby uzależnienie od ekstremalnego indeksu , co oznacza, że ma PZ zamiast . Wzięcie wykonałoby zadanie dla dowolnego . Biorąc pod uwagę sekwencję iid rv.s po standardowym Frechecie, znane są metody generowania zależnej sekwencji ze standardowym marginesem Frecheta i zalecanym wskaźnikiem ekstremalnym . Co się jednak stanie, jeśli zastąpimy Bernoulli ''?

X_{k}

$X_k$

θ

$\theta$

(0 < θ < 1)

$(0 < \theta < 1)$

M_{n}

$M_n$

q^{n θ}

$q^{n\theta}$

q^{n}

$q^n$

n θ = a

$n\theta = a$

a > 0

$a>0$

X_{n}

$X_n$

X_{n}^{⋆}

$X_n^\star$

θ

$\theta$ standard Frechet'' by

— Yves

(+1) Bardzo fajnie, @Yves. Kilka uwag: ( 1 ) Pierwszą część można postrzegać jako uzupełnienie przyjętego przeze mnie podejścia. W rzeczywistości, kiedy po raz pierwszy dostałem serię w , gdy natychmiast zobaczyłem związek z geometrią, najpierw spróbowałem czegoś bardziej bezpośredniego i nie wpadłem na naturalny sposób, aby to zrobić. Twoja odpowiedź rozwiązuje ten problem. ( 2 ) Twoje podejście można również wdrożyć, używając tylko monety . W szczególności można wygenerować, schodząc w dół pełnego drzewa binarnego opartego na rzetelnych rzutach monetą, gdzie lewe węzły są liśćmi, a decyzję podejmuje (...)

b_{n} q^{n}

$b_n q^n$

B e r (p)

$\mathrm{Ber}(p)$

N

$N$

— kardynał

(...) porównując liczbę w skonstruowaną z (częściowej) sekwencji przewracania monet do sum częściowych . Głębokość rozwiązanie daje . ( 3 ) Uważam, że , co zmieni twoje wnioski dotyczące skończoności średniej. ( 4 ) Zarówno w twoim, jak i moim podejściu wydaje się, że nie możemy uniknąć obliczeń , nawet jeśli pozwolimy na jednolite losowe zmienne oprócz naszej monety do celów pobierania próbek. Znalezienie sposobu na uniknięcie tego wydaje się najbardziej oczywistym sposobem na poprawę wydajności.

(0, 1)

$(0,1)$

f_{n} = \sum_{i = 1}^{n} b_{i}

$f_n = \sum_{i=1}^n b_i$

N

$N$

n b_{n} \sim c n^{- (1 + a)}

$n b_n \sim c n^{-(1+a)}$

f_{n} = \sum_{i = 1}^{n} b_{i}

$f_n = \sum_{i=1}^n b_i$

B e r (p)

$\mathrm{Ber}(p)$

— kardynał

1

Dziękuję @cardinal. Zgadzam się ze wszystkimi twoimi komentarzami, z wyjątkiem być może (3). Naprawdę popełniłem błąd, ponieważ wynosi (edytowany), ale wykładnik wydaje się właściwy. Użyłem reprezentacji znalezionej np. Na stronie Wikipedii na nieskończonym produkcie i wziąłem co daje odpowiednik dla produktu . Byłbym bardziej pewny, gdybyś mógł to sprawdzić.

c

$c$

- 1 / Γ (- a)

$-1/\Gamma(-a)$

n

$n$

Γ (z)

$\Gamma(z)$

z := - a

$z:=-a$

\prod_{k = 1}^{n - 1}

$\prod_{k=1}^{n-1}$

— Yves

Drogi @ Yves, (+1) Masz rację co do stałej i około (3). Przepraszam. W jakiś sposób, kiedy zacząłem przepisywać rzeczy na papier, ostatecznie skupiłem się na asymptotykach zamiast . :-)

b_{n}

$b_n$

n b_{n}

$n b_n$

— kardynał