逆関数法とは？その導出と乱数生成のための実装

この記事では、サンプリング法の１つである逆関数法を導出する。

逆関数法では、乱数を生成したい累積分布関数の逆関数にもとづき、標準一様乱数を任意の乱数に変換する。

参考にした資料はPRML（パターン認識と機械学習―ベイズ理論による統計的予測：C.M.ビショップ著）であり、同書「第11章サンプリング法 11.1.1 標準的な分布」の前半の解説としても有用であることを期待している。

Index

逆関数法とは？
導出
Pythonによる実装

逆関数法とは？

逆関数法は、（コンピュータに標準実装されている）標準一様分布 $U(0,1)$ にしたがう乱数によって、非一様分布に従う乱数を生成する手法であり、その際に、生成したい分布の分位関数（累積分布関数の逆関数）を用いる。

すなわち、非一様分布 $p$ に関して

$$z=h(y)\equiv\int_{-\infty}^{y}p(t)dt$$

と置き、乱数 $z$ を標準一様分布から生成する。

その後

$$y=h^{-1}(z)$$

と変換すると、変換後の $y$ は $p$ に従う乱数となる。

導出

$y,z$ がしたがう確率密度関数をそれぞれ $p_{y}(y),p_{z}(z)$ と書く。

$y,z$ の分布はそれぞれ異なるため、ここでは確率密度関数が同一ではないことを明示した。

これらを用いて $y\leftarrow z$ の変数変換を考える。

微小な区間 $(z, z+\delta z)$ に入る観測値は、変換後は同じく微小な区間 $(y, y+\delta y)$ に入るので

$$p_{y}(y)\delta y\simeq p_{z}(z)\delta z$$

$$p_{y}(y)=p_{z}(z)\left|\frac{dz}{dy}\right|\quad(\because p_{y}(y)\geq 0, p_{z}(z)\geq 0)$$

が成り立つ。

とくに $z$ が標準一様分布 $U(z|0,1)$ にしたがうとき、確率密度関数 $p_{z}(z)=1$ となる（このことについては

一様分布の定義・性質とその証明

定義連続変数 $x$ が、有限区間 $x\in$ で一様分布する場合、その確率密度関数 $p(x)$ は $$p(x)=\frac{1}{b-a}$$ で定義される。より一般に、分布が一様分布であることを明示して $$p(x...

を参照）ので

$$p_{y}(y)=\left|\frac{dz}{dy}\right|$$

である。

ここで $z=h(y)$ と置くと

$$p_{y}(y)=\left|\frac{dh(y)}{dy}\right|$$

$$p_{y}(y)=|h'(y)|\tag{1}$$

となる。

上式が成立するための $h(y)$ の条件について考える。ここで、 $y$ の定義域を $(-\infty,\infty)$ としても一般性を失わない。

まず、「 $h(y)$ は１回微分可能である」…①。

次に、「 $h: y\rightarrow z$ は $(-\infty,\infty)\rightarrow(0,1)$ の写像である」…②。

最後に、この関数 $z=h(y)$ について後から逆関数 $y=h^{-1}(z)$ を考えるので、「 $h: y\rightarrow z$ は全単射である」…③。

以上の①～③を満たす関数として、新たな分布 $p_{t}(t)$ （注：これは後に $p_{y}(y)$ に等しいことが示されるが、現時点ではまだ未定である。なお、 $t$ の定義域も $(-\infty,\infty)$ とする）を用いて

$$z=h(y)\equiv\int_{-\infty}^{y} p_{t}(t)dt$$

と置ける。

この時、 $h(y)$ は $y$ についての単調増加関数である。

これを式 $(1)$ に代入すると

$$p_{y}(y)=\left|\frac{d}{dy}\left(\int_{-\infty}^{y}p_{t}(t)dt\right)\right|$$

$$p_{y}(y)=|p_{t}(y)|=p_{t}(y)$$

より、結果的に $p_{t}=p_{y}$ であることがわかる。（定積分の微分については

定積分を微分する公式とその証明

公式 $f(t)$ を $t$ に関する１変数関数、$p(x), q(x)$ をそれぞれ $x$ に関する１変数関数とすると$$\frac{d}{dx}\int_{p(x)}^{q(x)}f(t)dt=f(q(x))q'(x...

を参照）

以上の結果から、

$$h(y)\equiv\int_{-\infty}^{y} p_{y}(t)dt$$

は $y$ の累積分布関数であり、その逆関数を求めて

$$y=h^{-1}(z)$$

という変換を行うことで、標準一様分布に従う乱数 $z$ から非一様分布 $p_{y}$ に従う乱数 $y$ を生成できることが示された。

Pythonによる実装

import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import expon # 指数分布

z = np.random.random(100000) # 標準一様分布
y = expon.ppf(z) # 分位関数（＝パーセント点関数：percent point function）
plt.hist(y, bins=1000)
plt.title("random exp")
plt.xlabel("y")
plt.savefig("exp.png")
plt.show()

逆関数法とは？

導出

Pythonによる実装

Comments