手帳と試行

学んだことをアウトプットしていきます。 日々、ノートあるのみ。

モンテカルロ積分

数値積分法に乱択アルゴリズムを組み合わせたモンテカルロ積分を導入する。

モンテカルロ積分

ある確率変数 xx の確率密度関数が p(x)p(x) である場合、xx を何らかの関数 gg に入力すると、その出力値 g(x)g(x) も確率的な値になる。

x:stochastic    g(x):stochastic\begin{aligned} &x : && \text{stochastic} \\ \implies &g(x) : && \text{stochastic} \end{aligned}

このとき、出力値の期待値は f(x)f(x)00 まわりの 11 次モーメントとして与えられる。

E[g(x)]g(x)=dxg(x)p(x)\begin{aligned} \mathop \mathrm E[g(x)] \coloneqq \langle g(x) \rangle = \int dx g(x)p(x) \end{aligned}

もしも確率密度関数 p(x)p(x)台 (support) X\mathcal X である、すなわち

p(x)={p(x)xX0otherwise\begin{aligned} p(x) = \left\{\begin{aligned} & p(x) && x \in {\cal X} \\ & 0&& \text{otherwise} \end{aligned}\right. \end{aligned}

であるならば、期待値は次式で計算される。

E[g(x)]=Xdxg(x)p(x)\begin{aligned} \mathop \mathrm E[g(x)] = \int_\mathcal X dx g(x) p(x) \end{aligned}

ここで g(x)=f(x)p(x)g(x) = \dfrac{f(x)}{p(x)} という書き換えを施すと次式が得られる。

E[f(x)p(x)]=Xdxf(x)p(x)p(x)=Xdxf(x)\begin{aligned} \mathop \mathrm E \left[ \frac{f(x)}{p(x)} \right] = \int_\mathcal X dx \frac{f(x)}{p(x)}p(x) = \int_\mathcal X dx f(x) \end{aligned}

標本平均 f(x)p(x)\overline{\dfrac{f(x)}{p(x)}} の期待値 E[f(x)p(x)]\mathop \mathrm E \left[ \overline{\dfrac{f(x)}{p(x)}} \right] は期待値 E[f(x)p(x)]\mathop \mathrm E \left[ \dfrac{f(x)}{p(x)} \right] に一致するから、左辺は次のように近似することができる。

E[f(x)p(x)]f(x)p(x)=1Ni=1Nf(xi)p(xi)\mathop \mathrm E \left[ \frac{f(x)}{p(x)} \right] \simeq \overline{\dfrac{f(x)}{p(x)}} = \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{p(x_i)}

これにより次のような数値積分の計算式が得られる。

Xdxf(x)1Ni=1Nf(xi)p(xi)xip(x)\begin{aligned} \int_\mathcal X dx f(x) \simeq \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{p(x_i)} \end{aligned} \\ x_i \sim p(x)

このような乱択アルゴリズムを組み込んだ数値積分の方法をモンテカルロ積分 (Monte Carlo integration) という。

精度

ここで気になるのは、次の近似がどれほど正しいといえるかである。

E[g(x)]g(x)\begin{aligned} \mathop \mathrm E[g(x)] \simeq \overline{g(x)} \end{aligned}

これを定量的に計るには、標本平均 g(x)\overline{g(x)} の標準誤差 V[g(x)]\sqrt{\mathop{\rm V}[\overline{g(x)}]} を使用すればよい。この量は標準偏差 V[g(x)]\sqrt{\mathop{\rm V}[g(x)]}1N\sqrt\dfrac{1}{N} 倍であるから、オーダーは

V[g(x)]=V[g(x)]N=O(1N)\begin{aligned} \sqrt{\mathop{\rm V}[\overline{g(x)}]} = \sqrt\frac{\mathop{\rm V}[g(x)]}{N} = \mathcal O \left( \sqrt\frac{1}{N} \right) \end{aligned}

となる。よって

E[g(x)]=g(x)+O(1N)\begin{aligned} \mathop{\rm E}[g(x)] = \overline{g(x)} + \mathcal O \left( \sqrt\frac{1}{N} \right) \end{aligned}

と見積もれる。この精度はそのまま積分の精度に伝播する。

Xdxf(x)=1Ni=1Nf(xi)p(xi)+O(1N)\begin{aligned} \int_\mathcal X dx f(x) = \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{p(x_i)} + \mathcal O \left( \sqrt\frac{1}{N} \right) \end{aligned}

これは、精度を nn 倍にしたいならば、サンプルの大きさを n2n^2 倍にする必要があるということを主張している。

連続一様分布からのサンプル

たとえば p(x)p(x) として1次元連続一様分布 U(a,b)\mathcal U(a, b) を仮定する。

p(x)=U(a,b)1ba\begin{aligned} p(x) = \mathcal U(a, b) \coloneqq \frac{1}{b-a} \end{aligned}

この場合、p(x)p(x) の台は X=[a,b]\mathcal X = [a, b] であるから、次のような数値積分公式が得られる。

abdxf(x)=1Ni=1Nf(xi)1ba+O(1N)=baNi=1Nf(xi)+O(1N)xiU(xa,b)\begin{aligned} \int_a^b dx f(x) &= \frac{1}{N} \sum_{i=1}^N \frac{f(x_i)}{\frac{1}{b-a}} + {\cal O} \left( \sqrt\frac{1}{N} \right) \\ &= \frac{b-a}{N} \sum_{i=1}^N f(x_i) + {\cal O} \left( \sqrt\frac{1}{N} \right) \end{aligned} \\ x_i \sim \mathcal U(x | a, b)