Symuluj zmienną Bernoulliego z prawdopodobieństwem

9

Czy ktoś może mi powiedzieć, jak symulować , gdzie , za pomocą rzutu monetą (tyle razy, ile potrzebujesz) z ?Bernoulli(ab)a,bNP(H)=p

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

Abrakadabra
źródło
1
Czy to pytanie pochodzi z kursu lub podręcznika? Jeśli tak, dodaj [self-study]tag i przeczytaj jego wiki . Pamiętaj, że na końcu pytania nie musisz prosić o pomoc - wiemy, że każdy, kto tu zamieszcza, ma nadzieję na pomoc!
Silverfish
1
Jest tu świetny post autorstwa @Glen_b (choć nie pamiętam gdzie) o tym, dlaczego nie ma czegoś takiego jak „stronnicza moneta z prawdopodobieństwem p", ale wiem, że to tylko kwestia peryferyjna twojego pytania!
Silverfish,
2
@dsaxton Pytanie brzmi „tyle, ile potrzebujesz”; będzie skończony z prawdopodobieństwem 1, ale nie będzie ograniczony (możesz przekroczyć dowolną ustaloną liczbę rzutów), ale sprzeciw na tej podstawie byłby jak powiedzenie „podrzucaj uczciwą monetę, dopóki nie zdobędziesz głowy” nie jest opłacalne jako metoda generowania geometrii (12losowe liczby.
Glen_b
1
@AbracaDabra Czy to ćwiczenie dla klasy? Jeśli nie, jak powstaje?
Glen_b
1
@Glen_b: To nie jest ćwiczenie z mojej klasy. Właśnie przyszło mi do głowy w tym łańcuchu myśli: zgodnie z klasycznym prawdopodobieństwem weź uczciwą monetę, gdy zwiększasz liczbę rzutów, stosunek#Heads#tailszbiega się w połowie. Tak też musi być w przypadku stronniczości ... Oznacza to, że moneta zbiega się do określonej liczby, potrzebujeszP(H)być tą liczbą. Teraz pomyślałem, co jeśli chcemy wyprodukować liczbę, ale mamy monetęP(H)inny numer (znany lub nieznany)?
AbracaDabra,

Odpowiedzi:

8

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. KiedyU<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 pczęść przedziału; gdy są to ogony, generuj rekurencyjnieX od ostatniego 1pczęść 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ą1p. 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żyjemyp=10.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

[x,y)=[x,x+(yx)p)[x+(yx)p,y)=s(I,H)s(I,T).

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

t<I

znaczy że t jest dolną granicą dla I: t<x dla wszystkich xI. Podobnie,t>I znaczy t jest górną granicą dla I.

pisać a/b=t. (W rzeczywistości nie ma znaczenia, czytjest rzeczywisty zamiast racjonalnego; wymagamy tylko tego0t1.)

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

  1. Zestaw n=0 i In=I0=[0,1).

  2. Podczas (tIn) {Rzuć monetą, aby wyprodukować Xn+1. ZestawIn+1=S(In,Xn+1). Przyrost n.}

  3. Gdyby t>In+1 następnie ustaw Z=1. W przeciwnym razie ustawZ=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 pomocniczejs. 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. Porysujmy10,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.0095klapki były główkami. Chociaż niższy niż cel0.01, wynik Z na poziomie 0.5154nie jest znaczący: to odchylenie można przypisać przypadkowi. Średnia liczba przerzutów wyniosła9.886- trochę mniej niż dziesięć. Gdybyśmy skorzystali z1p monety, średnia byłaby 0.0094- nadal nie różni się znacząco od celu, ale tylko 1.177 przerzuty byłyby potrzebne średnio.

Whuber
źródło
W mojej odpowiedzi nie mogę nie dostrzec podobieństw między tym rozwiązaniem a rozwiązaniem 2. Podczas gdy zakładam obiektywną monetę (PS naprawdę interesujące rozwiązanie problemu stronniczych monet) i wykonuję wszystkie obliczenia / porównania w bazie 2, wszystkie obliczenia / porównania wykonujesz w bazie 10. Jakie są twoje przemyślenia?
Cam.Davidson.Pilon
1
@cam Myślę, że moje przykłady mogą Cię zwieść: chociaż używają ładnych liczb w bazie 10, konstrukcja nie ma nic wspólnego z żadną konkretną bazą.
whuber
2
(+1) Bardzo ładna rozdzielczość. Optymalizacja obejmuje górną i dolną granicęa/b przez moce takie jak pn(1p)m i / lub (n+mm)pn(1p)m. Byłoby miło znaleźć optymalną dychotomię pod względem liczby symulowanych Bernoullis.
Xi'an
5

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ć Bernoulli(ab), Mogę wymyślić dwa rozwiązania (pierwsze nie jest moje, ale drugie jest uogólnieniem):

Rozwiązanie 1

Odwróć obiektywną monetę bczasy. Gdybyagłowy nie są obecne, zacznij od nowa. Gdybyagłowice obecne, zwróć, czy pierwsza moneta jest głowicą, czy nie (ponieważP(first coin is heads | a heads in b coins)=ab)

Rozwiązanie 2

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

Stworzymy nowy numer binarny za pomocą rzutów monetą. Zacząć od0.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ślibin(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)

Cam.Davidson.Pilon
źródło
Dobra odpowiedź, ale czy możesz wyjaśnić swoją metodą 1, jak zrobić dla irracjonalnego p?
AbracaDabra,
2
@AbracaDabra dlaczego miałaby to racjonalność pmateria?
Glen_b
@AbracaDabra: niezależnie od wartości p, prawdopodobieństwo uzyskania (0,1) i (1,0) są takie same, a mianowicie p(1p), stąd prawdopodobieństwo, że uda się uzyskać jeden przeciw drugiemu 1/2.
Xi'an
4

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 zdarzeniaH, i T. [Jest na to proste podejście, które wymaga pary rzutówH=(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ą wa i niższy o (ba) (i rozpocznij proces od początku) i uruchom losowy spacer aż do absorpcji, prawdopodobieństwo absorpcji przy górnej barierze wynosi aa+(ba)=ab.

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

Glen_b - Przywróć Monikę
źródło