[Scientific / technical calculation by Python] Generation of non-uniform random numbers giving a given probability density function, Monte Carlo simulation

Find a non-uniform random number with a probability density function $ g (t) $ given by a uniform random number $ x $ of $ [0,1] $. ** Computational statistical mechanics using the Monte Carlo method often requires the generation of random numbers with a specific distribution. ** This article is the basis for that.

General: Functions for generating non-uniform random numbers

Outline it. The purpose here is to obtain an intuitive understanding without pursuing mathematical rigor. I would like readers to forgive me on that point and point out any inadequacies in the content ...

Now, Probability density function of uniform random numbers in interval [0,1]

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

From the relationship between $ t $ and $ x $ to give a non-uniform random number with an arbitrary probability density function $ g (t (x)) $,

t=w(x) \tag{2}

I want to ask.

First, from the request for "preservation of probability" f(x)dx = g(t)dt \tag{3}

Is. Using equation (1) dx = g(t)dt \Leftrightarrow x(t) = \int_{−∞}^{t} g(t')dt'+C \tag{3'}

To get. Here, $ C $ is an integral constant. As a sufficient condition, $ x (-∞) = C = 0 $. At this time, $ \ int _ {−∞} ^ {∞} g (t') dt'= 1 $, so $ x (∞) = 1 $.

From the above, the relationship between ** x and t is **

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

You can see that they are connected by **. $ g (t) $ is what you give. From this, if t is calculated as the relation of x and $ w (x) $ in equation (2) is determined (if the inverse function is calculated), that t becomes a random variable giving a non-uniform distribution g (t). The above is the procedure for generating random numbers according to an arbitrary probability distribution $ g (t) $. ** **

Contents

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

Find a non-uniform random number according to.

(2) Find a non-uniform random number according to a Gaussian distribution with a mean value of $ t_ {av} $ and a standard deviation of $ \ sigma $. ** Since it is difficult to find the inverse function by the usual method, convert the function density in the two-dimensional plane to polar coordinates Box-Muller method % 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) is used. ** **

(1): Exponential distribution

Equation (4) is used.

x(t)=\int_{−∞}^{t} g(t')dt'=\int_{0}^{t} exp(-t') dt' = 1- exp(-t) Therefore, if t is calculated as a function of x from here,

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

To get. ** Given a uniform random number of [0,1] to x, equation (5) determines the distribution of t, which are distributed as an exponential function g (t) = exp (-t) **. Let's check it with the following program.

"""
x=[0,1]Random numbers that follow an exponential non-uniform distribution from uniform random numbers 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]Uniform random number data score of(N->At ∞, the frequency distribution becomes the probability density function)

for n in range(N):
    x=random()
    t=-log(1-x) # x->Conversion to 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]Creating a standardized bar graph that is equally divided into Nbins within the range of

#Theoretical probability density distribution 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()

Result (1)

Using equation (5), t is generated from a sample of 1 million points with x = [0,1] and the frequency distribution of it is drawn. It can be seen that the theoretical curve $ g (t) = exp (-t) $ is reproduced.

(2) Gaussian distribution

Gaussian distribution with mean $ t_ {av} $, standard deviation $ \ sigma $, $ g (t) = (2 \ pi \ sigma ^ 2) ^ {-1 / 2} \ exp (\ frac {-(t) -t_ {av}) ^ 2} {2 \ sigma ^ 2}) $ [Box-Muller Transform](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). From the two uniform random numbers x and y, in this case ** we get variables $ t_1 $ and $ t_2 $ that follow a non-uniform random number with two independent Gaussian distributions **.

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

Let's check this. The following code is a Monte Carlo simulation with $ \ sigma = 1 $, $ t_ {av} = 10 $.

"""
Non-uniform random numbers that generate a Gaussian distribution
box-Muller's method
"""

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]Uniform random number data score of(N->At ∞, the frequency distribution becomes the probability density function)

sigma=1.0 #standard deviation
ave=10   #Average value
for n in range(N):
    x1=random()
    x2=random()
    t1=sigma*((-2*log(x1))**(1/2))*(cos(2*pi*x2))+ave             #x->Conversion to t:
    t2=sigma*((-2*log(x1))**(1/2))*(sin(2*pi*x2))+ave             #x->Conversion to 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]Creating a standardized bar graph that is equally divided into Nbins within the range of
"""
[0,∞]If the Gaussian distribution of is normalized, the original product value is 0..Note that the number is doubled because it is 5 but 1!It is good to increase the average value.
"""
ax.hist(data_lis1, range=(0,20),bins=Nbins,normed=True) # [0,5]Creating a standardized bar graph that is equally divided into Nbins within the range of
#ax.hist(data_lis2, range=(0,20),bins=Nbins,normed=True) # [0,5]Creating a standardized bar graph that is equally divided into Nbins within the range of


#Theoretical probability density distribution
xx=np.linspace(0,20,100)
yy=(1/sqrt(2*pi*sigma**2))*np.exp((-(xx-ave)**2)/(2*sigma**2)) #Gaussian distribution
plt.plot(xx,yy, label='Theory')
plt.legend(loc='best')

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

plt.show()

Result (2)

The orange line is the theoretical curve (Gaussian function). It can be seen that the frequency distribution reproduces the Gaussian distribution.

Addendum

It often happens that the inverse function is not found. In that case, it is necessary to consider numerical evaluation and other methods. [Metropolis Method](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) is often used in computational physics.

Postscript: August 11, 2017

I posted an article about Monte Carlo simulation using the Metropolis method. [Scientific / technical calculation by Python] Generation of non-uniform random numbers given a given probability density function, Monte Carlo simulation