Simulieren Sie eine Bernoulli-Variable mit der Wahrscheinlichkeit

Kann mir jemand sagen, wie man simuliert , wobei mit einem Münzwurf (so oft Sie möchten) mit &$\mathrm{Bernoulli}\left({a\over b}\right)$$a,b\in \mathbb{N}$$P(H)=p$

Ich dachte über die Verwendung von Ablehnungsproben nach, konnte sie aber nicht festnageln.

Da es unzählige Lösungen gibt, sollten wir eine effiziente finden .

Die Idee dahinter beginnt mit einer Standardmethode zur Implementierung einer Bernoulli-Variablen: Vergleichen Sie eine einheitliche Zufallsvariable mit dem Parameter . Wenn , geben Sie ; Andernfalls geben Sie .a / b U < a / b 1 $U$$a/b$$U \lt a/b$$1$$0$

Wir können die Münze als einheitlichen Zufallszahlengenerator verwenden$p$ . Um innerhalb eines beliebigen Intervalls gleichmäßig eine Zahl zu erzeugen , werfen Sie die Münze. Wenn es Köpfe sind, erzeugen Sie rekursiv einen einheitlichen Wert im ersten Teil des Intervalls; Wenn es sich um Schwänze handelt, generieren Sie rekursiv aus dem letzten Teil des Intervalls. Irgendwann wird das Zielintervall so klein, dass es nicht wirklich darauf ankommt, wie Sie eine Zahl daraus auswählen: So beginnt die Rekursion. Es ist offensichtlich, dass dieses Verfahren gleichmäßige Variationen erzeugt (bis zu jeder gewünschten Präzision), wie durch Induktion leicht bewiesen werden kann.[ x , y ) X p X 1 - $U$$[x, y)$$X$$p$$X$$1-p$

Diese Idee ist nicht effizient, führt aber zu einer effizienten Methode. Da Sie in jeder Phase eine Zahl aus einem bestimmten Intervall zeichnen, sollten Sie zunächst prüfen, ob Sie sie überhaupt zeichnen müssen. Wenn Ihr Zielwert außerhalb dieses Intervalls liegt, kennen Sie bereits das Ergebnis des Vergleichs zwischen dem Zufallswert und dem Ziel. Daher neigt dieser Algorithmus dazu, schnell zu enden. (Dies könnte als das in der Frage angeforderte Ablehnungsstichprobenverfahren ausgelegt werden .)$[x,y)$

Wir können diesen Algorithmus weiter optimieren. Zu jedem Zeitpunkt haben wir tatsächlich zwei Münzen, die wir verwenden können: Durch Umbenennen unserer Münze können wir daraus eine machen, die Köpfe mit einer Chance von . Daher können wir als Vorberechnung rekursiv auswählen, welche Neukennzeichnung zu der geringeren erwarteten Anzahl von Flips führt, die für die Beendigung benötigt werden. (Diese Berechnung kann ein teurer Schritt sein.)$1-p$

Zum Beispiel ist es ineffizient, eine Münze mit zu verwenden, um eine Bernoulli -Variable direkt zu emulieren : Es dauert durchschnittlich fast zehn Flips. Wenn wir jedoch eine Münze verwenden, sind wir mit nur zwei Flips sicher fertig und die erwartete Anzahl von Flips beträgt nur .( 0,01 ) p = 1 - 0,0 = 0,1 $p=0.9$$(0.01)$$p=1-0.0=0.1$$1.2$

Hier sind die Details.

Partitionieren Sie jedes gegebene halboffene Intervall in die Intervalle$I = [x, y)$

Dies definiert die zwei Transformationen und , die in halboffenen Intervallen arbeiten.s ( ∗ , T $s(*,H)$$s(*,T)$

Aus terminologischen Gründen, wenn eine Menge von reellen Zahlen lass den Ausdruck$I$

bedeuten , dass eine untere Schranke für ist : für alle . In ähnlicher Weise bedeutet dass eine Obergrenze für .$t$$I$$t \lt x$$x \in I$$t \gt I$$t$$I$

Schreiben Sie . (Tatsächlich macht es keinen Unterschied, ob real statt rational ist; wir benötigen nur )$a/b = t$$t$$0 \le t \le 1$

Hier ist der Algorithmus, um eine Variable mit dem gewünschten Bernoulli-Parameter zu erzeugen :$Z$

1. Setze und .$n=0$$I_n = I_0 = [0,1)$

2. Während {Wirf die Münze, um zu erzeugen . Setze Inkrement .}$(t\in I_{n})$$X_{n+1}$$I_{n+1} = S(I_n, X_{n+1}).$$n$

3. Wenn dann setze . Andernfalls setzen Sie .$t \gt I_{n+1}$$Z=1$$Z=0$

Implementierung

Zur Veranschaulichung ist hier eine RImplementierung des Alorithmus als Funktion draw. Seine Argumente sind der Zielwert und das Intervall , anfänglich . Es verwendet die Hilfsfunktion, die implementiert . Obwohl dies nicht erforderlich ist, wird auch die Anzahl der Münzwürfe verfolgt. Es gibt die Zufallsvariable, die Anzahl der Würfe und das zuletzt untersuchte Intervall zurück.$t$$[x,y)$$[0,1)$s$s$

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))
}

Nehmen Sie als Beispiel für seine Verwendung und Prüfung seiner Genauigkeit den Fall und . Zeichnen wir mit dem Algorithmus Werte, berichten über den Mittelwert (und seinen Standardfehler) und geben die durchschnittliche Anzahl der verwendeten Flips an.$t=1/100$$p=0.9$$10,000$

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

In dieser Simulation waren der Flips Köpfe. Obwohl niedriger als das Ziel von , ist der Z-Score von nicht signifikant: Diese Abweichung kann dem Zufall zugeschrieben werden. Die durchschnittliche Anzahl der Flips betrug - etwas weniger als zehn. Wenn wir die Münze verwendet hätten, wäre der Mittelwert immer noch nicht signifikant anders als das Ziel, aber wären nur Flips erforderlich gewesen.$0.0095$$0.01$$-0.5154$$9.886$$1-p$$0.0094$$1.177$

Hier ist eine Lösung (etwas chaotisch, aber es ist mein erster Stich). Sie können tatsächlich ignorieren und WLOG nimmt 1/2 an . Warum? Es gibt einen cleveren Algorithmus , um aus zwei vorgespannten Münzwürfen einen unvoreingenommenen Münzwurf zu generieren. Wir können also 1/2 annehmen .$P(H) = p$$P(H)=1/2$$P(H)=1/2$

Um einen zu generieren , kann ich mir zwei Lösungen vorstellen (die erste ist nicht meine eigene, aber die zweite ist eine Verallgemeinerung):$\text{Bernoulli}(\frac{a}{b})$

Lösung 1

Wirf die unbefangene Münze mal um. Wenn Köpfe vorhanden sind, beginnen Sie von vorne. Wenn Kopf ist vorhanden, Rückkehr , ob die erste Münze ist ein Kopf oder nicht (weil )$b$$a$$a$$P(\text{first coin is heads | $a$ heads in $b$ coins}) = \frac{a}{b}$

Lösung 2

Dies kann auf einen beliebigen Wert von . Schreiben Sie in binärer Form. Beispiel:$\text{Bernoulli}(p)$$p$$0.1 = 0.0001100110011001100110011... \text{base 2}$

Wir werden eine neue Binärzahl mit Münzwürfen erstellen. Beginnen Sie mit Fügen Sie Ziffern hinzu, je nachdem, ob Kopf (1) oder Zahl (0) angezeigt werden. Vergleichen Sie bei jedem Flip Ihre neue Binärzahl mit der Binärdarstellung von bis zur gleichen Ziffer . Schließlich werden die beiden auseinander gehen und zurückkehren, wenn größer als Ihre Binärzahl ist.p b i n ( p $0.$$p$ $bin(p)$

In Python:

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

Einige Beweise:

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

ist ungefähr 0,4 (jedoch nicht schnell)

Ich sehe eine einfache Lösung, aber zweifellos gibt es viele Möglichkeiten, einige davon vermutlich einfacher. Dieser Ansatz kann in zwei Schritte unterteilt werden:

1. Generieren aus zwei Ereignissen mit gleicher Wahrscheinlichkeit bei einem unfairen Münzwurfverfahren (die Kombination der jeweiligen Münze und der Methode, mit der sie geworfen wird, um einen Kopf mit der Wahrscheinlichkeit erzeugen ). Wir können diese beiden gleich wahrscheinlichen Ereignisse und . [Hierfür gibt es einen einfachen Ansatz, bei dem zwei Wurfpaare und , um zwei gleich wahrscheinliche Ergebnisse zu erzielen, wobei alle anderen Ergebnisse zur Erzeugung eines neuen Paares führen von Brötchen, um es erneut zu versuchen.]H * T * H * = ( H , T ) T * = ( T , H $p$$H^*$$T^*$$H^*=(H,T)$$T^*=(T,H)$

2. Jetzt erzeugen Sie mit der simulierten fairen Münze einen zufälligen Spaziergang mit zwei absorbierenden Zuständen. Durch Auswahl des Abstands der absorbierenden Zustände vom Ursprung (einer darüber und einer darunter) können Sie die Absorptionswahrscheinlichkeit einstellen, indem Sie beispielsweise den oberen absorbierenden Zustand als ein gewünschtes Verhältnis von ganzen Zahlen festlegen. Insbesondere wenn Sie die obere Absorptionsbarriere bei und die untere bei (und den Prozess vom Ursprung aus starten) und den Random Walk bis zur Absorption ausführen, beträgt die Wahrscheinlichkeit der Absorption an der oberen Barriere .- ( b - a ) $a$$-(b-a)$$\frac{a}{a+(b-a)} = \frac{a}{b}$

   (Es müssen hier einige Berechnungen durchgeführt werden, um dies zu zeigen, aber Sie können die Wahrscheinlichkeiten ziemlich einfach ermitteln, indem Sie mit Wiederholungsrelationen arbeiten ... oder Sie können dies tun, indem Sie unendliche Reihen summieren ... oder es gibt andere Möglichkeiten.)