[Calcul scientifique / technique par Python] Génération de nombres aléatoires non uniformes donnant une fonction de densité de probabilité donnée, simulation Monte Carlo

Trouvez un nombre aléatoire non uniforme avec une fonction de densité de probabilité $ g (t) $ donnée par un nombre aléatoire uniforme $ x $ de $ [0,1] $. ** La mécanique statistique computationnelle utilisant la méthode de Monte Carlo nécessite souvent la génération de nombres aléatoires avec une distribution particulière. ** Cet article est la base pour cela.

Général: fonctions de génération de nombres aléatoires non uniformes

Décrivez-le. Le but ici est d'obtenir une compréhension intuitive sans poursuivre la rigueur mathématique. J'aimerais que les lecteurs me pardonnent sur ce point et signalent toute insuffisance dans le contenu ...

Maintenant, Fonction de densité de probabilité de nombres aléatoires uniformes dans l'intervalle [0,1]

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

À partir de la relation entre $ t $ et $ x $ pour donner un nombre aléatoire non uniforme avec une fonction de densité de probabilité arbitraire $ g (t (x)) $,

t=w(x) \tag{2}

Je veux demander.

Premièrement, de la demande de "préservation de la probabilité" f(x)dx = g(t)dt \tag{3}

Est. Utilisation de l'équation (1) dx = g(t)dt \Leftrightarrow x(t) = \int_{−∞}^{t} g(t')dt'+C \tag{3'}

Obtenir. Ici, $ C $ est une constante d'intégration. Comme condition suffisante, soit $ x (-∞) = C = 0 $. A ce moment, $ \ int_ {−∞} ^ {∞} g (t ') dt' = 1 $, donc $ x (∞) = 1 $.

D'après ce qui précède, la relation entre ** x et t est **

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

Vous pouvez voir qu'ils sont connectés par **. $ g (t) $ est ce que vous donnez. À partir de là, si t est calculé comme la relation de x et $ w (x) $ dans l'équation (2) est déterminée (si la fonction inverse est calculée), cette t devient une variable stochastique donnant une distribution non uniforme g (t). Ce qui précède est la procédure pour générer des nombres aléatoires qui suivent une distribution de probabilité arbitraire $ g (t) $. ** **

Contenu

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

Trouvez des nombres aléatoires non uniformes selon.

(2) Trouvez un nombre aléatoire non uniforme selon une distribution gaussienne avec une valeur moyenne de $ t_ {av} $ et un écart type de $ \ sigma $. ** Comme il est difficile de trouver la fonction inverse par la méthode habituelle, convertissez la densité de fonction dans le plan bidimensionnel en coordonnées polaires [méthode Box-Muller](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) est utilisé. ** **

(1): Distribution exponentielle

L'équation (4) est utilisée.

x(t)=\int_{−∞}^{t} g(t')dt'=\int_{0}^{t} exp(-t') dt' = 1- exp(-t) Par conséquent, si t est calculé en fonction de x à partir d'ici,

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

Obtenir. ** Étant donné un nombre aléatoire uniforme de [0,1] à x, l'équation (5) détermine la distribution de t, qui sont distribuées comme une fonction exponentielle g (t) = exp (-t) **. Vérifions-le avec le programme suivant.

"""
x=[0,1]Génération de nombres aléatoires qui suivent une distribution exponentielle non uniforme
"""

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

import matplotlib.pyplot as plt

data_lis=[]
N=10**6  # [0,1]Nombre aléatoire uniforme de points de données(N->En ∞, la distribution de fréquence devient une fonction de densité de probabilité)

for n in range(N):
    x=random()
    t=-log(1-x) # x->Conversion en 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]Création d'un graphique à barres standardisé également divisé en Nbins dans la plage de

#Distribution de densité de probabilité théorique 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()

Résultat (1)

A partir d'un échantillon de 1 million de points avec x = [0,1], t est généré en utilisant l'équation (5) et sa distribution de fréquence est dessinée. On voit que la courbe théorique $ g (t) = exp (-t) $ est reproduite.

(2) Distribution gaussienne

Distribution gaussienne avec moyenne $ t_ {av} $, écart-type $ \ sigma $, $ g (t) = (2 \ pi \ sigma ^ 2) ^ {-1 / 2} \ exp (\ frac {- (t) -t_ {av}) ^ 2} {2 \ sigma ^ 2}) $ [Méthode Box-Muller](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). À partir des deux variables de nombres aléatoires uniformes x et y, dans ce cas ** nous obtenons les variables $ t_1 $ et $ t_2 $ qui suivent des nombres aléatoires non uniformes avec deux distributions gaussiennes indépendantes **.

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

Vérifions ça. Le code suivant est une simulation de Monte Carlo avec $ \ sigma = 1 $, $ t_ {av} = 10 $.

"""
Nombres aléatoires non uniformes qui génèrent une distribution gaussienne
boîte-Méthode Muller
"""

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]Nombre aléatoire uniforme de points de données(N->En ∞, la distribution de fréquence devient une fonction de densité de probabilité)

sigma=1.0 #écart-type
ave=10   #Valeur moyenne
for n in range(N):
    x1=random()
    x2=random()
    t1=sigma*((-2*log(x1))**(1/2))*(cos(2*pi*x2))+ave             #x->Conversion en t:
    t2=sigma*((-2*log(x1))**(1/2))*(sin(2*pi*x2))+ave             #x->Conversion en 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]Création d'un graphique à barres standardisé également divisé en Nbins dans la plage de
"""
[0,∞]Si la distribution gaussienne de est normalisée, la valeur du produit d'origine est 0..Notez que le nombre est doublé car il est de 5 mais 1!Il est bon d'augmenter la valeur moyenne.
"""
ax.hist(data_lis1, range=(0,20),bins=Nbins,normed=True) # [0,5]Création d'un graphique à barres standardisé également divisé en Nbins dans la plage de
#ax.hist(data_lis2, range=(0,20),bins=Nbins,normed=True) # [0,5]Création d'un graphique à barres standardisé également divisé en Nbins dans la plage de


#Distribution de densité de probabilité théorique
xx=np.linspace(0,20,100)
yy=(1/sqrt(2*pi*sigma**2))*np.exp((-(xx-ave)**2)/(2*sigma**2)) #Distribution gaussienne
plt.plot(xx,yy, label='Theory')
plt.legend(loc='best')

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

plt.show()

Résultat (2)

La ligne orange est la courbe théorique (fonction gaussienne). On constate que la distribution de fréquence reproduit la distribution gaussienne.

Addenda

Il arrive souvent que la fonction inverse ne soit pas trouvée. Dans ce cas, il est nécessaire d'envisager une évaluation numérique et d'autres méthodes. [Méthode Metropolis](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) est souvent utilisé en physique computationnelle.

Postscript: 11 août 2017

J'ai posté un article sur la simulation Monte Carlo en utilisant la méthode métropole. [Calcul scientifique / technique par Python] Génération de nombres aléatoires non uniformes donnant une fonction de densité de probabilité donnée, simulation Monte Carlo