任意の確率分布に従う標本生成：逆関数法

科学ライブラリを用いれば、殆ど全ての既知の確率密度関数に従う標本は生成できるが、ここでは累積分布関数が与えられている任意の確率分布に従うi一次元標本を、連続一様分布$${U(0,1)}$$に従う標本を用いて生成する。
連続型確率密度関数を$${f(x)}$$とし、累積分布関数$${u=\int^{x}f(t)dt=F(x)}$$とする。累積分布関数の逆関数$${x=F^{-1}(u)}$$と与えられる。
この累積分布関数は、ある変数$${a< b}$$に対し、必ず$${F(a)<F(b)}$$であることから、確率変数$${x}$$が、確率変数領域内のある$${v}$$に対し、$${x\leq v}$$となる確率$${Pr(x\leq v)}$$は、
$${Pr(x\leq v)=Pr(F^{-1}(u)\leq v)=Pr(u \leq F(v))=\displaystyle{\int^{F(v)}_0 du=F(v)}}$$
となる。
よって、連続一様分布$${U(0,1)}$$に従う標本を$${u}$$が与えられれば、確率密度関数$${f(x)}$$に従う確率変数$${x}$$は、累積分布関数$$の逆関数を用いて、$${x=F^{-1}(u)}$$で与えられる。
これをラプラス分布について行えば、$${f(x)=\displaystyle{\frac{1}{2}e^{-|x|}}}$$より、累積分布関数$$は、以下の手順で得られる。
$${t>0}$$の時、
$${\displaystyle{\frac{1}{2}\int^t_{-\infty}e^{-x} dx=\frac{1}{2} \{ \int_{-\infty}^0 e^x dx +\int^t_0 e^{-x}dx\} =\frac{1}{2} \{ 1+[-e^{-x}]^t_0\} =\frac{1}{2}(2-e^{-t})}}$$
$${t\leq 0}$$の時
$${\displaystyle{\frac{1}{2}\int^t_{-\infty}e^{-x} dx=\frac{1}{2}\int^t_{-\infty}e^x dx = \frac{1}{2}[e^x]^t_{\infty}=\frac{1}{2}e^t}}$$
これをまとめれば、
$${F(x)=\displaystyle{\frac{1}{2}\Big(1+sign(x)(1-e^{-|x|})\Big)}}$$
この逆関数は、
$${F^{-1}(u)=\displaystyle{-sign(u-\frac{1}{2})\log(1-2\Big|u-\frac{1}{2}\Big|)}}$$
これをpythonで実装し、numpy.randomのlaplaceで発生させた標本と比べるコードは以下のようになる。

import numpy as np
from numpy import random as rand
import matplotlib.pyplot as plt
rand.seed(0)
randnum=1_000_000

u=rand.random(randnum)
scale=1.0
lap1=rand.laplace(loc=0.0, scale=scale, size=randnum)
x=-np.sign(u-0.5)*np.log(1.-2.*np.abs(u-0.5))
plt.hist(x, density=True, bins='auto',alpha=1.0,label=r'$F^{-1}$(u)')
plt.hist(lap1, density=True, histtype="step",bins='auto',alpha=1.0,rwidth=1.3)  
plt.ylabel('Probability')
plt.legend()
plt.xlabel('Data');

Laplace distribution