月別アーカイブ: 2013年10月

任意の確率密度分布に従う乱数の生成（von Neumannの棄却法）

投稿日時: 投稿者:

任意の分布に関して乱数を作りたかったのでメモ。
分布の逆関数を容易に求められる場合は、逆関数法を使ってしまうのが簡単なのだけど、

  1. どんな（一価の）関数にも適用出来る
  2. アルゴリズムが作りやすい

といった理由でここでは棄却法を使います。

以下のページを参考にしています。
数値計算法 2011/6/29 (pdf)
15 乱数とモンテカルロ法 (pdf)

von Neumannの棄却法

考え方はとてもシンプルで、

  1. 分布関数を含む範囲内で一様乱数のペア (x_{i}, y_{i})を生成

    2. 具体的には分布関数の範囲を[x_{min}, x_{max}], [0, y_{max}]、
    確率変数 x_{rand}, y_{rand} を0 \leq x_{rand}, y_{rand} \leq 1 に従う一様分布と仮定した時に、

    x_{i} = x_{min} + (x_{max} - x_{min})x_{rand}\\ y_{i} = y_{max} y_{rand}\\

    を作る。

  2. 生成した乱数ペア(x_{i}, y_{i}) が分布関数内に入っていれば、x_{i}を分布に従う乱数として採用（下図でいうと赤点が採用、×が不採用）

random_test

アルゴリズム的にはこんな感じ？（コード中にあるテストファイルはこちらから→”random_dist.txt“）

#!/opt/local/bin/python
import numpy as np
from scipy.interpolate import interp1d

def randomgen(func,func_range,func_max):
    while(1):
        x_rand, y_rand = np.random.random(2)
        x = func_range[0] + (func_range[1]-func_range[0]) * x_rand
        if ( func_max * y_rand <= func(x) ):
            return x

random_dist = np.loadtxt('random_dist.txt') # load function
func = interp1d(random_dist[:,0],random_dist[:,1],kind='linear') # interpolate
func_range = [min(random_dist[:,0]),max(random_dist[:,0])] # holizontal range [x_min, x_max]
func_max = max(random_dist[:,1]) # vertical range

# generate random number list
xlist = []
for i in range(1000):
    xlist.append(randomgen(func,func_range,func_max))

# plot result
histresult = hist(xlist, range=[func_range[0],func_range[1]], bins=30,  color='red', alpha=0.5)
plot(random_dist[:,0],random_dist[:,1]*max(histresult[0])/max(random_dist[:,1])),color='blue',linewidth=2.0)

結果のヒストグラム：
random_dist

正規分布（normal distribution）と指数関数（exponential distribution）の畳み込み積分

投稿日時: 投稿者:

手計算でチェックしたかったので過程をメモ。

以下の二つの関数を畳み込み積分した場合、

f(x) = \begin{cases} \lambda \exp(-\lambda x) & (x\ge 0) \\ 0 & (x<0) \end{cases}
\displaystyle g(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}

Exponentially modified Gaussian distribution
にあるようにex-Gaussian distributionと呼ばれる分布が出来る。
\displaystyle f(x)\ast g(x) = \frac{\lambda}{2}e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\mbox{erfc}\left(\frac{\mu+\lambda \sigma^2 -x}{\sqrt{2}\sigma}\right)
ここでerfcはガウスの相補誤差関数（complementary error function）で以下のように定義されている。erfはガウスの誤差関数。
\displaystyle \mbox{erfc}(x) = 1-\mbox{erf}(x) = \frac{2}{\sqrt{\pi}} \int ^{+\infty}_{x} e^{-t^2}dt

Convolution of the normal and exponential probability density functions

以下計算、
\displaystyle f(x)\ast g(x) = \int^{+\infty}_{-\infty} dx'f(x')g(x-x')\\[1.0ex] =\int^{+\infty}_{0} dx' \lambda \exp(-\lambda x') \frac{1}{\sqrt{2\pi}\sigma}\exp\left\{-\frac{(x-x'-\mu)^2}{2\sigma^2}\right\}\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} \int^{+\infty}_{0} dx' \exp\left[-\frac{1}{2\sigma^2} \left\{x'^2+2x'(\mu-x+\lambda \sigma^2)+(x-\mu)^2 \right\} \right]\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} \int^{+\infty}_{0} dx'\exp\left\{-\frac{1}{2\sigma^2}(x'-x+\mu+\lambda \sigma^2)^2 + \frac{1}{2}\lambda^2 \sigma^2 -\lambda x + \lambda \mu \right\}\\[1.0ex] = \frac{\lambda}{\sqrt{2\pi}\sigma} e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\int^{+\infty}_{0} dx'\exp\left\{-\frac{1}{2\sigma^2}(x'-x+\mu+\lambda \sigma^2)^2 \right\}

ここで、変数変換
\displaystyle z = \frac{1}{\sqrt{2}\sigma}(x'-x+\mu+\lambda \sigma^2), \; dz' = \frac{dx'}{\sqrt{2}\sigma}
とガウスの相補誤差関数の定義を使えば、
\displaystyle f(x)\ast g(x) = \frac{\lambda}{\sqrt{\pi}} e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\int^{+\infty}_{\frac{\mu+\lambda \sigma^2-x}{\sqrt{2}\sigma}} dz e^{-z^2}\\[1.0ex] = \frac{\lambda}{2}e^{\lambda \mu + \frac{1}{2}\lambda^2\sigma^2}e^{-\lambda x}\mbox{erfc}\left(\frac{\mu+\lambda \sigma^2 -x}{\sqrt{2}\sigma}\right)\\
が得られる。

Python code

Pythonに用意されてるモジュールについてはnumpy.convolveも参照

結果はこんな感じ。
convolution_test