[PYTHON] Erstellen Sie eine Zufallszahl mit einer beliebigen Wahrscheinlichkeitsdichte

TL;DR

Um eine kontinuierliche stochastische Variable $ \ hat {X} $ gemäß der Wahrscheinlichkeitsdichte $ f (x) $ zu erstellen, müssen die kumulative Verteilungsfunktion $ F (x) $ und $ 0 \ leq \ hat {R} <1 $ eine einheitliche Zufallszahl sein Sie können $ \ hat {R} $ verwenden, um $ \ hat {X} = F ^ {-1} (\ hat {R}) $ ([inverse Funktionsmethode](https: //ja.wikipedia.) Zu erstellen. org / wiki /% E9% 80% 86% E9% 96% A2% E6% 95% B0% E6% B3% 95)).

Einführung

Es gibt Zeiten, in denen Sie eine stochastische Variable mit einer beliebigen stochastischen Dichteverteilung erstellen möchten. Ein typisches Beispiel ist die Erzeugung von "Punkten, die gleichmäßig auf einer sphärischen Oberfläche verteilt sind". Dies ist beispielsweise erforderlich, wenn Sie die Richtung im dreidimensionalen Raum zufällig angeben möchten, obwohl der absolute Wert der Geschwindigkeit gleich ist.

Im Allgemeinen sind die von einem Programm erhaltenen (Pseudo-) Zufallszahlen einheitliche Zufallszahlen von 0 bis 1. In diesem Artikel erfahren Sie, wie Sie damit Zufallszahlen mit beliebigen Wahrscheinlichkeitsdichteverteilungen und einigen konkreten Beispielen erstellen.

Wahrscheinlichkeitsdichtefunktion und kumulative Verteilungsfunktion

Angenommen, Sie haben eine Wahrscheinlichkeitsvariable $ \ hat {X} $, die einen kontinuierlichen Wert annimmt. Die Wahrscheinlichkeit, dass diese Variable einen Wert zwischen $ a $ und $ b $ annimmt

P(a \leq \hat{X} < b) = \int_a^b f(x) dx

Wenn wir schreiben können, wird $ f (x) $ als Wahrscheinlichkeitsdichtefunktion bezeichnet. Dies schränkt den Bereich ein

P(x \leq \hat{X} < x+dx) \sim f(x)dx

Sie können lesen "Die Wahrscheinlichkeit, dass $ \ hat {X} $ einen Wert nahe dem Wert $ x $ annimmt, ist proportional zu $ f (x) $". In diesem Sinne ist die Bedeutung der Wahrscheinlichkeitsdichtefunktion leicht zu verstehen, aber theoretisch und numerisch ist die kumulative Verteilungsfunktion, die die Wahrscheinlichkeitsdichtefunktion integriert, häufig einfacher zu verwenden.

Die Definition der kumulativen Verteilungsfunktion lautet wie folgt.

P(\hat{X} < x) = F(x)

Das heißt, $ F (x) $ ist die Wahrscheinlichkeit, dass die Wahrscheinlichkeitsvariable $ \ hat {X} $ einen Wert von weniger als $ x $ annimmt. Beim Schreiben mit der Wahrscheinlichkeitsdichtefunktion

F(x) = \int_{-\infty}^x f(x) dx

Daher ist die kumulative Verteilungsfunktion das Integral der Wahrscheinlichkeitsdichtefunktion.

Zum Beispiel die Wahrscheinlichkeitsdichtefunktion und die kumulative Verteilungsfunktion von einheitlichen Zufallszahlen von 0 bis 1

\begin{align}
f(x) &= 1 \\
F(x) &= x
\end{align}

(Jedoch $ 0 \ leq x <1 $).

Beweis

Wie eingangs erwähnt, wird zum Erstellen einer Wahrscheinlichkeitsvariablen $ \ hat {X} $ mit einer Wahrscheinlichkeitsdichte $ f (x) $ eine einheitliche Zufallszahl $ \ hat {R von $ 0 \ leq \ hat {R} <1 $ erstellt } $ Konvertiert mit der Umkehrfunktion $ F ^ {-1} (x) $ der kumulativen Verteilungsfunktion

\hat{X} = F^{-1}(\hat{R})

Es ist in Ordnung, wenn Sie es tun. Um dies zu beweisen, zeigen Sie, dass die kumulative Verteilungsfunktion der Wahrscheinlichkeitsvariablen $ \ hat {X} $, die aus einheitlichen Zufallszahlen erstellt wurde, $ F (x) $ ist, oder umgekehrt die Wahrscheinlichkeitsdichtefunktion $ f (x). Es ist in Ordnung zu zeigen, dass die Konvertierung der Wahrscheinlichkeitsvariablen $ \ hat {X} $ mit $ in $ \ hat {R} = F (\ hat {X}) $ eine einheitliche Zufallszahl von 0 bis 1 ergibt. Persönlich denke ich, dass die letztere Tatsache wichtiger ist, aber ich werde sie in beide Richtungen zeigen.

Beweis (vorwärts)

Eine einheitliche Zufallszahl von 0 bis 1 $ \ hat {R} $ wird unter Verwendung der Umkehrfunktion der kumulativen Verteilungsfunktion $ F (x) $ konvertiert, und eine neue Wahrscheinlichkeitsvariable $ \ hat {X} $ wird erhalten.

\hat{X} = F^{-1}(\hat{R})

Machen durch. Zu diesem Zeitpunkt ist die Wahrscheinlichkeit, dass $ \ hat {X} $ kleiner als der Wert $ x $ ist, [^ 3]

[^ 3]: Hier nehmen wir an, dass die kumulative Verteilungsfunktion eine Eins-zu-Eins-Zuordnung ist (Beweis erforderlich).

\begin{align}
P(\hat{X} < x) &= P(F^{-1}(\hat{R}) < x) \\
&= P(\hat{R} < F(x))
\end{align}

Die kumulative Verteilungsfunktion $ F_R (r) $ der einheitlichen Zufallszahl $ \ hat {R} $ ist

F_R(r) \equiv P(\hat{R} < r) = r

Damit

P(\hat{R} < F(x)) = F(x)

Von Oben,

P(\hat{X} < x) = F(x)

Dies bedeutet, dass die kumulative Verteilungsfunktion von $ \ hat {X} $ $ F (x) $ ist, dh die Wahrscheinlichkeitsdichtefunktion ist $ f (x) $. Das wollte ich beweisen.

Beweis (umgekehrte Richtung)

Angenommen, Sie erhalten eine Wahrscheinlichkeitsvariable $ \ hat {X} $ mit einer Wahrscheinlichkeitsdichte von $ f_X (x) $. Die kumulative Verteilungsfunktion für diese Variable ist $ F_X (x) $ und ihre Definition ist

F_X(x) = P(\hat{X} < x)

ist. Lassen Sie uns nun diese stochastische Variable mit der Umkehrfunktion $ F_X (x) $ dieser kumulativen Verteilungsfunktion transformieren, um eine neue stochastische Variable $ \ hat {Y} $ zu erstellen.

\hat{Y} = F_X(\hat{X})

Aus der Definition der kumulativen Verteilungsfunktion ergibt sich $ 0 \ leq \ hat {Y} <1 $.

Die Wahrscheinlichkeit, dass diese Wahrscheinlichkeitsvariable $ \ hat {Y} $ kleiner als $ y $ ist, dh die kumulative Verteilungsfunktion $ F_Y (y) $

F_Y(y) = P(\hat{Y} < y)

ist. Nun, weil $ \ hat {Y} = F_X (\ hat {X}) $

F_Y(y) = P(F_X(\hat{X}) < y)

Hier ist die kumulative Verteilungsfunktion $ F_X (x) $ eine monoton ansteigende Funktion, so dass sie als Abbildung reversibel ist [^ 2]. Wenn also $ F_X (x) = y $ ist, dann ist $ x = F_X ^ {-1} (y) $. Damit

[^ 2]: Genau genommen muss ich viel sagen, aber vergib mir, weil es lästig ist.

P(F_X(\hat{X}) < y) = P(\hat{X} < F_X^{-1}(y))

Kann gemacht werden. Dies ist die Definition der kumulativen Verteilungsfunktion selbst, da die rechte Seite die Wahrscheinlichkeit ist, dass $ \ hat {X} $ einen Wert von weniger als $ x = F_X ^ {-1} (y) $ hat. Deshalb

P(\hat{X} < F_X^{-1}(y)) = F_X(F_X^{-1}(y)) = y

Von Oben,

F_Y(y) = y

ist geworden. Da $ 0 \ leq \ hat {Y} <1 $ ist, ist die Wahrscheinlichkeitsvariable $ \ hat {Y} $ eine einheitliche Zufallszahl von 0 bis 1. Mit anderen Worten: ** Wenn Sie eine stochastische Variable mit einer beliebigen Wahrscheinlichkeitsdichte mit einer kumulativen Verteilungsfunktion konvertieren, ist dies eine einheitliche Zufallszahl von 0 bis 1 **.

Wenn Sie dagegen eine einheitliche Zufallszahl mit $ F ^ {-1} (x) $ von 0 nach 1 konvertieren, erhalten Sie eine stochastische Variable mit einer Wahrscheinlichkeitsdichtefunktion $ f (x) $. Das wollte ich beweisen.

Beispiel

Punkte gleichmäßig auf der Oberfläche der Kugel verteilt

Angenommen, Sie möchten Punkte erstellen, die gleichmäßig auf der Oberfläche einer Einheitskugel verteilt sind. Die Koordinaten der Oberfläche der Einheitskugel werden unter Verwendung von zwei Winkeln $ \ theta und \ phi $ berechnet.

\begin{align}
x &= \sin{\theta} \cos{\phi}\\
y &= \sin{\theta} \sin{\phi}\\
z &= \cos{\theta}\\
\end{align}

Kann geschrieben werden. Ein häufiger Fehler besteht darin, $ \ theta $ und $ \ phi $ als einheitliche Zufallszahlen von $ 0 $ bis $ 2 \ pi $ zu verwenden. Machen wir das.

import numpy as np
import matplotlib.pyplot as plt
import random
from math import pi, sin, cos, sqrt
from mpl_toolkits.mplot3d import Axes3D

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  theta = random.uniform(0, pi)
  X.append(sin(theta)*cos(phi))
  Y.append(sin(theta)*sin(phi))
  Z.append(cos(theta))

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(X,Y,Z,marker="o",markersize=1,linestyle='None')
plt.show()

Das Ergebnis ist so.

Sie können sehen, dass es viele Punkte am Nord- und Südpol und wenige Punkte in der Nähe des Äquators gibt.

Wenn die Koordinatenvariablen als $ \ phi $ und $ \ theta $ genommen werden, kann zuerst $ \ phi $ eine einheitliche Zufallszahl sein, aber $ \ theta $ ist anders. Wenn Sie sich für einen Winkel $ \ theta $ entscheiden, ist dies ein Kreis der Kugel mit einem Radius von $ \ sin {\ theta} $. Wenn Sie Punkte auf der Kugel gleichmäßig machen möchten, müssen Sie Punkte proportional zur Länge dieses Umfangs machen.

Da die Länge des Umfangs proportional zu $ \ sin {\ theta} $ ist, ist die Wahrscheinlichkeitsdichtefunktion von $ \ theta $

f(\theta) = \frac{1}{2}\sin{\theta}

ist. Standardisierungsbedingungen

\int_0^\pi f(\theta) d\theta = 1

Sie können den Koeffizienten $ 1/2 $ schreiben, um zu erfüllen.

Verteilungsfunktion

F(\theta) = -\frac{1}{2}\cos{\theta} + \frac{1}{2}

Es wird sein. Um $ \ hat {\ Theta} $ zu generieren, um einer solchen kumulativen Verteilungsfunktion zu folgen, verwenden Sie die einheitliche Zufallszahl $ \ hat {R} $, um $ \ hat {\ Theta} = F ^ {-1} zu generieren Muss (\ hat {R}) $ sein. Wenn umgekehrt die Wahrscheinlichkeitsvariable $ \ hat {\ Theta} $ der Wahrscheinlichkeitsdichtefunktion $ f (\ theta) $ folgt

\hat{R} = -\frac{1}{2}\cos{\hat{\Theta}} + \frac{1}{2}

Die neue Wahrscheinlichkeitsvariable $ \ hat {R} $ wird als einheitliche Zufallszahl von 0 bis 1 erhalten.

Nun, was ich ursprünglich wollte, war nicht der Wert von $ \ theta $ selbst, sondern zum Beispiel $ \ cos {\ theta} $, die die $ z $ -Koordinate ist. Aus der obigen Formel können wir erkennen, dass $ \ cos {\ hat {\ Theta}} $ eine einheitliche Zufallszahl von -1 bis 1 ist. Mit anderen Worten, wenn Sie $ \ hat {Z} $ mit einer einheitlichen Zufallszahl von -1 bis 1 machen, wenn Sie $ \ hat {\ Theta} = \ cos ^ {-1} {\ hat {Z}} $ machen, Sie haben $ \ hat {\ Theta} $ mit der richtigen Wahrscheinlichkeitsdichte generiert. Auch aus $ \ cos ^ 2 {\ theta} + \ sin ^ 2 {\ theta} = 1 $, $ \ sin {\ theta} = \ sqrt {1-z ^ 2} $ und $ \ sin {\ theta Sie können auch den Wert von} $ finden.

Wenn das Obige implementiert ist, wird es so.

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  z = random.uniform(-1,1)
  X.append(sqrt(1-z**2)*cos(phi))
  Y.append(sqrt(1-z**2)*sin(phi))
  Z.append(z)

Das Ausführungsergebnis sieht so aus.

Es sieht so aus, als ob die Punkte gleichmäßig auf der Kugel verteilt sind.

Punkte gleichmäßig innerhalb der Einheitskugel verteilt

Die Punkte, die innerhalb der Einheitskugel gleichmäßig verteilt sind, sind

• Bestimmen Sie den Radius $ r (0 <r <1) $ entsprechend
• Verteilen Sie Punkte gleichmäßig auf einer Kugel mit einem Radius von $ r $

Es ist in Ordnung, wenn Sie die Prozedur befolgen. Um die Punkte auf einer Kugel mit einem Radius von $ r $ gleichmäßig zu verteilen, bestimmen Sie $ \ theta, \ phi $ wie zuvor.

\begin{align}
x &= r \sin{\theta} \cos{\phi}\\
y &= r \sin{\theta} \sin{\phi}\\
z &= r \cos{\theta}\\
\end{align}

Es ist in Ordnung, wenn Sie die Koordinaten des Punktes als festlegen. Das Problem ist die Verteilung von $ r $. Bei einer gleichmäßigen Verteilung erhöht sich die Punktedichte, wenn $ r $ klein ist.

Nun ist die Oberfläche einer Kugel mit einem Radius von $ r $ proportional zu $ r ^ 2 $. Natürlich muss die Anzahl der Punkte auch proportional zu $ r ^ 2 $ generiert werden. Wenn also die Wahrscheinlichkeitsdichtefunktion der Wahrscheinlichkeitsvariablen $ r $ die Standardisierungskonstante enthält

f(r) = 3 r^2

Es wird sein. Verteilungsfunktion

F(r) = r^3

ist. Dies ist intuitiv, da auf einer Kugel mit einem Radius von $ r $ oder weniger Punkte proportional zu $ r ^ 3 $ vorhanden sind.

Umkehrfunktion

F^{-1}(r) = r^{1/3}

ist. Daher ist es in Ordnung, eine einheitliche Zufallszahl $ r $ von 0 bis 1 zu verwenden und $ r ^ {1/3} $ als Radius festzulegen (der Radius r und die Zufallszahl r sind gleich, was verwirrend, aber verwirrend ist. Wird nicht).

Im Code sieht es so aus.

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  z = random.uniform(-1,1)
  r = random.uniform(0,1)
  X.append(r**(1.0/3.0)*sqrt(1-z**2)*cos(phi))
  Y.append(r**(1.0/3.0)*sqrt(1-z**2)*sin(phi))
  Z.append(r**(1.0/3.0)*z)

Das Ausführungsergebnis sieht so aus.

Vielleicht sind die Punkte innerhalb der Einheitskugel gleichmäßig verteilt.

Exponentialverteilung

Angenommen, es gibt ein Phänomen, das $ \ lambda $ mal pro Zeiteinheit auftritt. Die Verteilung der Zeit vom Beginn der Beobachtung bis zum ersten Auftreten des Phänomens ist eine Exponentialverteilung gemäß dem Parameter $ \ lambda $. Ein typisches Beispiel für eine Exponentialverteilung ist die Zerfallswahrscheinlichkeit radioaktiver Materialien. Darüber hinaus ist die Verteilung, wie viel Licht gerade gehen kann, ohne in einem Medium mit gleichmäßig verunreinigten Stoffen gestreut zu werden, ebenfalls eine exponentielle Verteilung.

Angenommen, Sie möchten ein Ereignis simulieren, das einer Exponentialverteilung folgt. Die Zeit vor dem Auftreten eines Ereignisses, das der Exponentialverteilung $ \ lambda $ folgt, $ \ hat {T} $ ist eine stochastische Variable. Da die Wahrscheinlichkeit, dass diese Zeit kleiner als $ t $ ist, die kumulative Verteilungsfunktion ist, ergibt sich aus der Definition der Exponentialverteilung

P(\hat{T}< t) \equiv F(t) = 1 - \exp(-t/\lambda)

Es wird sein. Umkehrfunktion

F^{-1}(t) = -\lambda \log(1-t)

Verwenden Sie also eine einheitliche Zufallszahl von 0 bis 1 $ \ hat {R} $

\hat{T} = -\lambda \log(1-\hat{R})

Dann erhalten Sie die gewünschte Wahrscheinlichkeitsvariable. Praktisch haben $ 1- \ hat {R} $ und $ \ hat {R} $ die gleiche Verteilung

\hat{T} = -\lambda \log(\hat{R})

Es wird oft gesagt.

Zusammenfassung

Ich habe vorgestellt, wie man aus einer einheitlichen Zufallszahl eine stochastische Variable mit einer beliebigen Wahrscheinlichkeitsdichtefunktion erstellt. Es scheint, dass es "inverse Funktionsmethode" genannt wird, weil es die inverse Funktion der kumulativen Verteilungsfunktion verwendet.

Die Methode, $ z = \ cos {\ theta} $ als einheitliche Zufallszahl zu erzeugen, wenn Punkte erzeugt werden, die gleichmäßig auf einer Kugel verteilt sind, ist gut eingeführt, aber es gibt nicht viele Beweise dafür, "warum das in Ordnung ist". Also habe ich es geschrieben. Das Wesentliche ist, dass "die Umwandlung einer beliebigen stochastischen Variablen mit ihrer kumulativen Verteilungsfunktion eine einheitliche Zufallszahl von 0 nach 1 ergibt."

Wenn die stochastischen Variablen eher diskret als kontinuierlich sind, kann die Walker-Alias-Methode einfach verwendet werden. In C ++ wird "std :: discrete_distribution" (wahrscheinlich) von Walkers Alias-Methode implementiert, daher denke ich, dass es einfach zu verwenden ist.