[Wissenschaftlich-technische Berechnung durch Python] Erzeugung ungleichmäßiger Zufallszahlen mit gegebener Wahrscheinlichkeitsdichtefunktion, Monte-Carlo-Simulation

Finden Sie eine ungleichmäßige Zufallszahl mit einer Wahrscheinlichkeitsdichtefunktion $ g (t) $, die durch eine einheitliche Zufallszahl $ x $ von $ [0,1] $ gegeben ist. ** Die rechnergestützte statistische Mechanik nach der Monte-Carlo-Methode erfordert häufig die Erzeugung von Zufallszahlen mit einer bestimmten Verteilung. ** Dieser Artikel ist die Basis dafür.

Allgemein: Funktionen zum Generieren ungleichmäßiger Zufallszahlen

Skizzieren Sie es. Der Zweck hier ist, ein intuitives Verständnis zu erlangen, ohne mathematische Strenge zu verfolgen. Ich möchte, dass die Leser mir in diesem Punkt vergeben und auf Unzulänglichkeiten im Inhalt hinweisen ...

Jetzt, Wahrscheinlichkeitsdichtefunktion einheitlicher Zufallszahlen im Intervall [0,1]

$ f (x) = 1 \ (0 \ le x \ le1), \ 0 $ (Andere) $ \ tag {1} $

Aus der Beziehung zwischen $ t $ und $ x $, um eine ungleichmäßige Zufallszahl mit einer beliebigen Wahrscheinlichkeitsdichtefunktion $ g (t (x)) $ zu erhalten,

t=w(x) \tag{2}

Ich will fragen.

Erstens aus der Forderung nach "Wahrung der Wahrscheinlichkeit" f(x)dx = g(t)dt \tag{3}

Ist. Unter Verwendung von Gleichung (1) dx = g(t)dt \Leftrightarrow x(t) = \int_{−∞}^{t} g(t')dt'+C \tag{3'}

Bekommen. Hier ist $ C $ eine Integrationskonstante. Als ausreichende Bedingung sei $ x (-∞) = C = 0 $. Zu diesem Zeitpunkt ist $ \ int_ {−∞} ^ {∞} g (t ') dt' = 1 $, also $ x (∞) = 1 $.

Aus dem Obigen ergibt sich die Beziehung zwischen ** x und t **

x(t)=\int_{−∞}^{t} g(t')dt'\tag{4}

Sie können sehen, dass sie durch ** verbunden sind. $ g (t) $ geben Sie. Wenn daraus berechnet wird, wie die Beziehung von x und $ w (x) $ in Gleichung (2) bestimmt wird (wenn die Umkehrfunktion erhalten wird), wird t eine stochastische Variable, die eine ungleichmäßige Verteilung g (t) ergibt. Das Obige ist das Verfahren zum Erzeugen von Zufallszahlen, die einer beliebigen Wahrscheinlichkeitsverteilung $ g (t) $ folgen. ** **.

Inhalt

(1) Exponentialverteilung, g(t)= $ exp(-t)$ (t>0), 0 (t< 0))

Finden Sie ungleichmäßige Zufallszahlen nach.

(2) Finden Sie eine ungleichmäßige Zufallszahl gemäß einer Gaußschen Verteilung mit einem Durchschnittswert von $ t_ {av} $ und einer Standardabweichung von $ \ sigma $. ** Da es mit der üblichen Methode schwierig ist, die Umkehrfunktion zu finden, konvertieren Sie die Funktionsdichte in der zweidimensionalen Ebene in Polarkoordinaten [Box-Muller-Methode](https://ja.wikipedia.org/wiki/%E3%83] % 9C% E3% 83% 83% E3% 82% AF% E3% 82% B9% EF% BC% 9D% E3% 83% 9F% E3% 83% A5% E3% 83% A9% E3% 83% BC % E6% B3% 95) wird verwendet. ** **.

(1): Exponentialverteilung

Gleichung (4) wird verwendet.

x(t)=\int_{−∞}^{t} g(t')dt'=\int_{0}^{t} exp(-t') dt' = 1- exp(-t) Wenn also t von hier aus als Funktion von x berechnet wird,

t=w(x)=- ln(1-x) \tag{5}

Bekommen. ** Bei einer einheitlichen Zufallszahl von [0,1] bis x bestimmt Gleichung (5) die Verteilung von t, die als Exponentialfunktion g (t) = exp (-t) verteilt sind **. Lassen Sie es uns mit dem folgenden Programm überprüfen.

"""
x=[0,1]Erzeugung von Zufallszahlen, die einer exponentiellen ungleichmäßigen Verteilung folgen, aus einheitlichen Zufallszahlen in
"""

from random import random
import numpy as np
from math import log

import matplotlib.pyplot as plt

data_lis=[]
N=10**6  # [0,1]Einheitliche Zufallszahl von Datenpunkten(N->Bei ∞ wird die Häufigkeitsverteilung zu einer Wahrscheinlichkeitsdichtefunktion)

for n in range(N):
    x=random()
    t=-log(1-x) # x->Umstellung auf t: t = -log(1-x)
    data_lis.append(t)


#plot
fig=plt.figure()
ax=fig.add_subplot(111)

Nbins=int(N/100)
ax.hist(data_lis, range=(0,5),bins=Nbins,normed=True) # [0,5]Erstellen eines standardisierten Balkendiagramms, das gleichmäßig in Nbins im Bereich von unterteilt ist

#Theoretische Wahrscheinlichkeitsdichteverteilung y= exp(-x)
xx=np.linspace(0,5,100)
yy=np.exp(-xx)
plt.plot(xx,yy, label='Theory: g (t) = exp(-t)')
plt.legend(loc='best')

plt.xlabel('t',size=24)
plt.ylabel('Frequency', size=24)

plt.show()

Ergebnis (1)

Aus einer Stichprobe von 1 Million Punkten mit x = [0,1] wird t unter Verwendung von Gleichung (5) erzeugt und die Häufigkeitsverteilung davon wird gezogen. Es ist ersichtlich, dass die theoretische Kurve $ g (t) = exp (-t) $ reproduziert wird.

(2) Gaußsche Verteilung

Gaußsche Verteilung mit Mittelwert $ t_ {av} $, Standardabweichung $ \ sigma $, $ g (t) = (2 \ pi \ sigma ^ 2) ^ {-1 / 2} \ exp (\ frac {- (t) -t_ {av}) ^ 2} {2 \ sigma ^ 2}) $ [Box-Muller-Methode](https://ja.wikipedia.org/wiki/%E3%83%9C%E3%83%83%E3%82%AF%E3%82%B9%EF%BC%9D % E3% 83% 9F% E3% 83% A5% E3% 83% A9% E3% 83% BC% E6% B3% 95). Aus den beiden einheitlichen Zufallszahlenvariablen x und y erhalten wir in diesem Fall ** Variablen $ t_1 $ und $ t_2 $, die ungleichmäßigen Zufallszahlen mit zwei unabhängigen Gaußschen Verteilungen folgen **.

t_1 = \sigma \ [(-2ln(x))^{1/2} \ cos(2 \pi y)+t_{av} ] t_2 = \sigma \ [(-2ln(x))^{1/2} \ sin(2 \pi y)+t_{av} ]

Lassen Sie uns dies überprüfen. Der folgende Code ist eine Monte-Carlo-Simulation mit $ \ sigma = 1 $, $ t_ {av} = 10 $.

"""
Ungleichmäßige Zufallszahlen, die eine Gaußsche Verteilung erzeugen
Box-Mullers Methode
"""

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

data_lis1=[]
#data_lis2=[]


N=10**5 # [0,1]Einheitliche Zufallszahl von Datenpunkten(N->Bei ∞ wird die Häufigkeitsverteilung zu einer Wahrscheinlichkeitsdichtefunktion)

sigma=1.0 #Standardabweichung
ave=10   #Durchschnittswert
for n in range(N):
    x1=random()
    x2=random()
    t1=sigma*((-2*log(x1))**(1/2))*(cos(2*pi*x2))+ave             #x->Umstellung auf t:
    t2=sigma*((-2*log(x1))**(1/2))*(sin(2*pi*x2))+ave             #x->Umstellung auf t:

    data_lis1.append(t1)
#    data_lis2.append(t2)



#plot
fig=plt.figure()
ax=fig.add_subplot(111)

Nbins=int(N/100)
#ax.hist(data_lis, range=(0,10),bins=Nbins,normed=True) # [0,5]Erstellen eines standardisierten Balkendiagramms, das gleichmäßig in Nbins im Bereich von unterteilt ist
"""
[0,∞]Wenn die Gaußsche Verteilung von normalisiert ist, ist der ursprüngliche Produktwert 0..Beachten Sie, dass die Zahl verdoppelt wird, da sie 5, aber 1 ist!Es ist gut, den Durchschnittswert zu erhöhen.
"""
ax.hist(data_lis1, range=(0,20),bins=Nbins,normed=True) # [0,5]Erstellen eines standardisierten Balkendiagramms, das gleichmäßig in Nbins im Bereich von unterteilt ist
#ax.hist(data_lis2, range=(0,20),bins=Nbins,normed=True) # [0,5]Erstellen eines standardisierten Balkendiagramms, das gleichmäßig in Nbins im Bereich von unterteilt ist


#Theoretische Wahrscheinlichkeitsdichteverteilung
xx=np.linspace(0,20,100)
yy=(1/sqrt(2*pi*sigma**2))*np.exp((-(xx-ave)**2)/(2*sigma**2)) #Gaußsche Verteilung
plt.plot(xx,yy, label='Theory')
plt.legend(loc='best')

plt.xlabel('t',size=24)
plt.ylabel('Frequency', size=24)

plt.show()

Ergebnis (2)

Die orange Linie ist die theoretische Kurve (Gaußsche Funktion). Es ist ersichtlich, dass die Häufigkeitsverteilung die Gaußsche Verteilung reproduziert.

Nachtrag

Es kommt oft vor, dass die Umkehrfunktion nicht gefunden wird. In diesem Fall müssen numerische Auswertungen und andere Methoden berücksichtigt werden. [Metropolis-Methode](https://ja.wikipedia.org/wiki/%E3%83%A1%E3%83%88%E3%83%AD%E3%83%9D%E3%83%AA%E3 % 82% B9% E6% B3% 95) wird häufig in der Computerphysik verwendet.

Nachtrag: 11. August 2017

Ich habe einen Artikel über die Monte-Carlo-Simulation mit der Metropolis-Methode veröffentlicht. Wissenschaftliche / technische Berechnung durch Python Erzeugung ungleichmäßiger Zufallszahlen mit gegebener Wahrscheinlichkeitsdichtefunktion, Monte-Carlo-Simulation http://qiita.com/sci_Haru/items/45ec84e6eb985f5c53d8)