[Python에 의한 과학·기술 계산] 주어진 확률 밀도 함수를 주는 비 균일 난수의 생성, 몬테카를로 시뮬레이션

$[0,1]$의 균일한 난수 $x$로부터 주어진 확률 밀도 함수 $g(t)$를 갖는 비균일한 난수를 구한다.
몬테카를로 방법을 이용한 계산 통계 역학에서는 특정 분포를 갖는 난수를 생성하는 것이 종종 필요하다. 이 논문은 이를 위한 기초가 된다.

일반론 : 비 균일 난수를 생성하는 함수

개략을 설명한다. 여기에서는 수학적 엄밀성은 추구하지 않고 직관적인 이해를 얻을 수 있는 것을 목적으로 하고 있다. 독자 여러분에게는 그 점에 관해 용서해 주시고, 또, 내용에 미비가 있으면 지적해 주셨으면 한다.

음,
구간 [0,1]의 균일 난수의 확률 밀도 함수

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

따라서 임의의 확률 밀도 함수 $ g (t (x)) $를 갖는 비 균일 난수를 제공하기 위해 $ t $와 $ x $의 관계,

$t=w(x)\tag{2}$

하고 싶습니다.

우선, "확률 보존"의 요청으로부터,
$f(x)dx = g(t)dt\tag{3}$

이다. 식 (1)을 사용하면
$dx = g(t)dt\Leftrightarrow x(t) =\int_{−∞}^{t} g(t')dt'+C\tag{3'}$

을 얻는다. 여기서 $ C $는 적분 상수입니다. 충분한 조건으로서 $x(-∞) = C=0$로 취한다. 이때 $\int_{−∞}^{∞} g(t')dt'=1$이므로 $x(∞) = 1$가 된다.

위에서 x와 t의 관계는

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

에 의해 묶여 있는 것을 알 수 있다. $g(t)$는 주는 것이다. 여기서 t를 x의 관계로 구해 식 (2)의 $w(x)$를 결정하면(역함수를 구하면), 그 t는 비균일 분포 g(t)를 줄 확률 변수가 된다. 이상이 임의의 확률 분포 $ g (t) $에 따른 난수를 발생시키는 절차이다.

내용

(1) 지수 함수적인 분포,
$g(t)$=
$ exp(-t)$ (t>0),
0 (t < 0))

에 따라 비 균일 난수를 찾습니다.

(2) 평균값이 $t_{av}$, 표준 편차가 $\sigma$인 가우스 분포에 따른 비 균일 난수를 구한다. 통상의 방법에서는 역함수를 구하기 어렵기 때문에, 2차원 평면에서의 함수 밀도를 극좌표로 변환하는 박스-뮐러법

(1): 지수 분포

식 (4)를 사용한다.

$x(t)=\int_{−∞}^{t} g(t')dt'=\int_{0}^{t} exp(-t') dt' = 1- exp(-t)$
따라서 여기에서 t를 x의 함수로 결정하면

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

을 얻는다. [0,1]의 균일한 난수를 x에 부여하면 식(5)에 의해 t의 분포가 결정되고, 그들은 지수함수 g(t)=exp(-t)로서 분포한다.
그것을 아래의 프로그램으로 확인하자.

"""
x=[0,1]の一様乱数から指数関数型の非一様分布に従う乱数の発生
"""

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

import matplotlib.pyplot as plt

data_lis=[]
N=10**6  # [0,1]の一様乱数のデータ点数 (N->∞で頻度分布が確率密度関数となる)

for n in range(N):
    x=random()
    t=-log(1-x) # x->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]の範囲でNbins等分した規格化された棒グラフの作成

# 理論確率密度分布 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()

결과(1)

x=[0,1]로 100만점 샘플한 것으로부터 식(5)를 이용하여 t를 생성해 그것의 빈도 분포를 그린 것. 이론 곡선 $g(t)=exp(-t)$를 재현하고 있는 것을 알 수 있다.

(2) 가우스 분포

평균값 $t_{av}$, 표준 편차 $\sigma$의 가우스 분포, $g(t)=(2\pi\sigma^2)^{-1/2}\exp(\frac{-(t -t_{av})^2}{2\sigma^2})$
박스-뮐러법에 의해 생성됩니다.
두 개의 균일 한 난수 변수 x와 y에서이 경우 두 개의 독립 가우스 분포를 갖는 비 균일 난수를 따르는 변수 $ t_1 $와 $ t_2 $를 얻습니다.

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

이것을 확인해 보자.
이하의 코드는 $\sigma=1$,$t_{av}=10$로서 몬테카를로 시뮬레이션을 실시한 것이다.

"""
ガウス分布を生成する非一様乱数
ボックス-ミュラーの方法
"""

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]の一様乱数のデータ点数 (N->∞で頻度分布が確率密度関数となる)

sigma=1.0 #標準偏差
ave=10   # 平均値
for n in range(N):
    x1=random()
    x2=random()
    t1=sigma*((-2*log(x1))**(1/2))*(cos(2*pi*x2))+ave             #x->t への変換:
    t2=sigma*((-2*log(x1))**(1/2))*(sin(2*pi*x2))+ave             #x->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]の範囲でNbins等分した規格化された棒グラフの作成
"""
[0,∞]のガウス分布にして正規化してしまうと，本来の積聞値が0.5なのに1にするため数値が頻度が二倍に勘定されてしまうことに注意! 平均値を大きめにすると良い.
"""
ax.hist(data_lis1, range=(0,20),bins=Nbins,normed=True) # [0,5]の範囲でNbins等分した規格化された棒グラフの作成 
#ax.hist(data_lis2, range=(0,20),bins=Nbins,normed=True) # [0,5]の範囲でNbins等分した規格化された棒グラフの作成 


# 理論確率密度分布 
xx=np.linspace(0,20,100)
yy=(1/sqrt(2*pi*sigma**2))*np.exp((-(xx-ave)**2)/(2*sigma**2)) # ガウス分布
plt.plot(xx,yy, label='Theory')
plt.legend(loc='best')

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

plt.show()

결과(2)

오렌지 선은 이론 곡선 (가우스 함수)입니다. 빈도 분포가 가우스 분포를 재현하고 있음을 알 수 있습니다.

부록

역함수를 찾을 수 없는 경우가 많습니다. 그 경우는 수치적인 평가나 다른 방법을 생각할 필요가 있다. 메트로폴리스법은 계산 물리학에서 자주 사용됩니다.

추가:2017년 8월 11일

메트로폴리스법을 이용한 몬테카를로 시뮬레이션의 기사를 투고했습니다.
[Python에 의한 과학 · 기술 계산] 주어진 확률 밀도 함수를주는 비 균일 난수 생성, 몬테카를로 시뮬레이션