機械学習に詳しくなりたいブログ

機械学習や数学について勉強した内容を中心に書きます。100%趣味です。記事は数学的に厳密でなかったり誤りを含んでいるかもしれません。ご指摘頂ければ幸いです。

逆変換法

一様乱数から任意の確率分布に従う乱数を得る手法の1つ、逆変換法を見ていきます。まだ完全に理解できていない部分があるのですが、、、。

積分の変数変換の式(6)より、


\displaystyle \int f(x)dx = \int f(x)\frac{dx}{dt}dt \tag{1}

でした。

zは区間[0,1]で一様分布g(z)に従い、p(y)に従うyを求めたいとします。式(1)より、


\displaystyle \int g(z)dz = \int g(z)\frac{dz}{dy}dy \tag{2}

です。g(z)p(y)は確率分布なので、全区間で積分すれば1になるはずですから、式(2)の左辺は1になり、 \displaystyle \int p(y)dy とも等しくなります。よって、


\displaystyle \int p(y)dy = \int g(z)\frac{dz}{dy}dy \tag{3}

ですから、


\displaystyle p(y) =  g(z)\left| \frac{dz}{dy} \right| \tag{4}

のように確率分布の変換式が求まります。ちゃんと理解できていないポイント1つ目が式(4)の絶対値。確率分布の変換において積分の方向は考えないということでしょうか?

さて、2変数以上への拡張は、ヤコビアンを用いて


\displaystyle p(y_1,\cdots,y_M) =  g(z_1,\cdots,z_M)\left| \frac{\partial(z_1,\cdots,z_M)}{\partial(y_1,\cdots,y_M)} \right| \tag{5}

と書けるようです。ここで\displaystyle \frac{\partial(z_1,\cdots,z_M)}{\partial(y_1,\cdots,y_M)}はヤコビ行列ではなくヤコビアンです。\displaystyle p(y_1,\cdots,y_M) =  g(z_1,\cdots,z_M)\left| \det J \right| と書いているのと等しいです。積分の変数変換で2変数以上になるとヤコビアンが出てきたのと同じ理屈と思います。

次にちゃんと理解できていないポイント2つ目。式(4)を積分することで


\displaystyle z = h(y) \equiv \int_{-\infty}^y p(\hat{y})d\hat{y} \tag{6}

が得られる、と手元に参考書に書いてあるのですが、一体どう変形したらこうなるのか?誰か教えてください。>< 

ということでいまいち理解できていないのですが、式(6)を使って確率分布を求めることは簡単です。今回仮にp(y)=\exp(-y)の指数分布とすると、式(6)より


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

これをyについて解けば


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

です。一様分布からzをサンプリングし、式(9)で変換すれば指数分布が求められるということです。試してみたところ、以下のグラフのようになり、確かに指数分布からサンプリングできていることが確認できました。逆変換法自体は、累積分布関数を図示すれば直感的に理解ができてそんなに難しくないのですが、式(4)を積分して式(6)になるっていう参考書の説明がとうとう理解できませんでした。うーん。

f:id:opabinia2:20180317215757p:plain

今回のコードです。

# 逆変換法
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()