スポンサーリンク

逆関数法とその導出

プログラミング

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

逆関数法とは?

逆関数法は、(コンピュータに標準実装されている)標準一様分布 \(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\) となる(このことについては

を参照)ので

$$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)\) は1回微分可能である」…①。

次に、「 \(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}\) であることがわかる。(定積分の微分については

を参照)

以上の結果から、

$$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()

コメント