逆関数法 - 機械学習に詳しくなりたいブログ

一様乱数を用いて、ある確率分布に従う乱数を得る手法の１つ、逆関数法を見ていきます。

逆関数法の直感的説明

イメージは図のとおりです。 $p(y)$ からサンプリングしたいが、その手段がないとします。 $h(y)$ は $p(y)$ の累積分布関数です。したがって、取る値の範囲は $[0,1$ ]です。 $z$ は区間 $[0,1$ ]の一様分布とし、いま、 $z^{\prime}$ がサンプリングできたとします。ここで、紫の線で示すように $h(y)$ の逆関数 $h^{-1}(y)$ によって $y^{\prime}$ が得られます。そしてこの $y^{\prime}$ は $p(y)$ に従って分布する、というものです。サンプリングしたい関数 $p(y)$ から累積分布関数を求め、さらにその逆関数が定まればサンプリングができるはずです。数学的な背景はなくとも直感的にもそんな気はします。確率密度が高いところで累積分布関数の傾きも大きくなるので、よりサンプリングされやすく、逆に確率密度の低いところでは累積分布関数の傾きは小さいのでサンプリングされにくくなります。ちなみに図の $p(y)$ は正規分布ですが、 $h(y)$ はフリーハンドで書いているので適当です、、。

証明

比較的簡単です。区間 $[0,1$ ]の一様分布に従う確率変数を $Z$ とします。逆関数によって $Y=h^{-1}(Z)$ が得られるとします。この $Y$ が、 $h(y)$ を累積分布関数に持つ確率分布に従っていてほしいのです。そこで任意の実数 $y$ に対して、 $\mathrm P(Y\le y)$ がどのように表せるのかを確認します。 $\mathrm P(Y\le y)$ はつまり累積分布関数ですから、これが $h(y)$ であってほしいわけです。さて、 $Y=h^{-1}(Z)$ でしたから、

$\mathrm P(Y\le y) = \mathrm P(h^{-1}(Z)\le y) \tag{1}$

です。

そして、上図より、 $h^{-1}(Z) \le Y$ は $Z \le h(y)$ と置き換えられますから、

$\mathrm P(h^{-1}(Z)\le y) = \mathrm P(Z \le h(y)) \tag{2}$

です。（累積分布関数は単調増加関数のため、大小が入れ替わることはありません）

ここで、 $Z$ が $[0,1$ ]の一様分布に従うとき、以下が成り立ちます。

$\mathrm P(Z\le z) = z \tag{3}$

確率変数 $Z$ が $z$ 以下の確率は、一様分布なら $z$ です。例えば $z=0.5$ とすれば、一様分布が0.5以下を取る確率は0.5ですよね。この関係を用いれば、式(2)の右辺は、

$\mathrm P(Z \le h(y)) = h(y) \tag{4}$

と書けます。以上より、

$\mathrm P(Y\le y) = h(y) \tag{5}$

となり、知りたかった累積分布関数 $\mathrm P(Y\le y)$ は $h(y)$ であることがわかりました。累積分布関数と確率分布は1対1ですから、つまり $Y$ は $h(y)$ を累積分布関数に持つ確率分布に従っているといえます。

サンプリング実験

早速試してみます。サンプリングしたい確率分布を $p(y)=\exp(-y)$ の指数分布とすると、その累積分布関数を求め、

$\begin{eqnarray*} \displaystyle z &=& \int_{-\infty}^y \exp(-y)dy \tag{6} \\ &=& -\exp(-y) + e^0\tag{7}\\ \end{eqnarray*}$

これを $y$ について解けば

$y = -\ln(1-z) \tag{8}$

です。一様分布から $z$ をサンプリングし、式(8)で変換すれば指数分布が求められるということです。試してみたところ、以下のグラフのようになり、確かに指数分布からサンプリングできていることが確認できました。

今回のコードです。

# 逆関数法
import numpy as np
import matplotlib.pyplot as plt


# サンプル数
N = 1000000

x = np.random.uniform(0, 1, N)
y = -np.log(1-x)

x2 = np.linspace(0, 5, 100)
y2 = np.exp(-x2)

plt.hist(y, normed=True, bins=50, alpha=0.5)
plt.plot(x2, y2)
plt.xlim(0,5)
plt.show()