(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。
モンテカルロ法の分散低減手法
「モンテカルロ法」で解説したように、モンテカルロの弱点として収束の遅さ($ 1 \;/\; \sqrt{N} $ に比例)がありました。ここでは基本的な分散低減手法を紹介します。
重点的サンプリング (Importance Sampling)
モンテカルロ法 - 式 $ (2) $ ではモンテカルロ法のサンプルに任意の確率密度関数(PDF)が使用できると書きました。モンテカルロ法ではこのPDFの選び方によって推定量の分散が変わるので、PDFの選び方は非常に重要な要素です。それでは一体どのようなPDFが良いのかというと、「被積分関数にできるだけ比例した形状のPDF」が良いPDFとなります。
被積分関数に完全に比例した次のようなPDFを考えます。
\begin{equation*}
p(x) = k f(x)
\end{equation*}
これを分散の計算式に代入してみると、以下の形になります。
\begin{eqnarray*}
\mathbb{\sigma}^2 \Brk{\langle I \rangle} &=& \frac{1}{N} \int_D \Prtc{\frac{f(x)}{p(x)} - I}^2 p(x) \d x
= \frac{1}{N} \int_D \Prtc{\frac{1}{k} - I}^2 k f(x) \d x \\
&=& \frac{\Prtb{\frac{1}{k} - I}^2}{N} \cdot k \int_D f(x) \d x
= \frac{\Prtb{\frac{1}{k} - I}^2}{N} \cdot k I \\
\end{eqnarray*}
ここで、$ p(x) = k f(x) $ の両辺を積分すると $ k = 1 \; / \; I $ が導かれるので
\begin{eqnarray*}
\mathbb{\sigma}^2 \Brk{\langle I \rangle} &=& 0
\end{eqnarray*}
となり、なんと分散が0になることがわかります。$ x_i $ がどの点をサンプルしようとも $ f(x_i) \;/\; p(x_i) $ が定数になるためです。ただ、現実的にモンテカルロ法を使うような場面では、このようなPDFはわかりません。なぜなら理想的なPDFは $ p(x) = f(x) \;/ \int_D f(x) \d x $ の形をとることからわかるように、そもそもPDFの定義に求めたい計算結果が含まれているためです。そこで如何に「できるだけ」比例したPDFからサンプルできるかが重要になってきます。下図のような3種類のPDFを考えると、最も分散が小さい = 結果の収束が速いのは右のPDFとなるでしょう。以上のように推定量にできるだけ比例したPDFに従ってサンプルすることを「重点的サンプリング (importance sampling)」と呼びます。
層化サンプリング (Stratified Sampling)
図2では区間中で一様分布に従って8つのサンプルを発生させていますが、明らかに区間の一部にサンプルが集まってしまっています。このようにモンテカルロ法では、何らかのPDFに従って $ N $ サンプルの分布を発生させたとしても、あくまでその分布が達成されるのは $ N \to \infty $ の場合であって、好ましくない分布になることがあります。サンプルの集中は低精度な推定結果を生み出してしまいます。この問題を解決するための手法に「層化サンプリング」があります。層化サンプリングでは積分区間をいくつかの小区間(strata)に分割し、各区間毎の積分の合計値を推定値とします(図3)。区間 $ \Brk{a, b} $ の1次元積分、分割数 $ m $ を例として式に表すと、次のようになります。 \begin{equation*} \int_a^b f(x) \d x = \int_a^{\alpha_1} f(x) \d x + \int_{\alpha_1}^{\alpha_2} f(x) \d x + \cdots + \int_{\alpha_{m - 2}}^{\alpha_{m - 1}} f(x) \d x + \int_{\alpha_{m - 1}}^{b} f(x) \d x \end{equation*}
各小区間毎で1つ以上のサンプルを発生させることで、サンプルが集中することがなくなり分散を低減することができます。各小区間毎に一様分布から $ n_i $ サンプルとった場合の推定量の分散は次のように計算されます。 \begin{eqnarray*} \mathbb{\sigma}^2 \Brk{\Abk{I}} &=& \sum_{i = 1}^{m} \frac{\alpha_i - \alpha_{i - 1}}{n_i} \int_{\alpha_{i - 1}}^{\alpha_i} f^2(x) \d x - \sum_{i = 1}^{m} \frac{1}{n_i} \Prtc{\int_{\alpha_{i - 1}}^{\alpha_i} f(x) \d x}^2 \end{eqnarray*} 各小区間が全て同じ幅( $ a_i - a_{i - 1} = (b - a) / m $ )、小区間毎のサンプル数は1つ( $ n_i = 1, N = m $ )とすると上の式は以下のように簡単化できます。 \begin{eqnarray*} \mathbb{\sigma}^2 \Brk{\Abk{I}} &=& \sum_{i = 1}^{m} \frac{b - a}{N} \int_{\alpha_{i - 1}}^{\alpha_i} f^2(x) \d x - \sum_{i = 1}^{m} \Prtc{\int_{\alpha_{i - 1}}^{\alpha_i} f(x) \d x}^2 \\ &=& \frac{b - a}{N} \int_{a}^{b} f^2(x) \d x - \sum_{i = 1}^{m} \Prtc{\int_{\alpha_{i - 1}}^{\alpha_i} f(x) \d x}^2 \end{eqnarray*} 純粋なモンテカルロ法の分散の式に、同じ条件を当てはめ式を整理してみるとわかりますが、第一項目は同じに、第二項目が $ (1 / N) \cdot I^2 $ となります。そしてこの項よりも上式の第二項目のほうが大きな値となるようです。したがって層化サンプリングを使えば、必ず分散が小さくなることが言え、また各小区間に2つ以上のサンプルを使うメリットが無いこともわかります。
層化サンプリングは使える箇所では積極的に使って行くべきではありますが、この手法は次元毎に $ N $ サンプル使う場合には $ d $ 次元積分で $ N^d $ サンプル必要となり、次元の呪いを受けてしまいます。層化サンプリングのように、サンプルの集中を抑えつつ、次元の呪いの解決も行った手法にラテン超方格法や準モンテカルロ法があります。
その他の手法
ここで紹介した重点的サンプリング、層化サンプリングの他にも多くの分散低減、またはそれに関連する手法が存在します。
- メトロポリスサンプリング(Metropolis Sampling)
- 制御変量法(Control Variates method)
- N-Rooks、ラテン超方格サンプリング(Latin Hypercube Sampling, LHS)
- 準モンテカルロ法(Quasi Monte Carlo method, QMC)
- 多重重点的サンプリング(Multiple Importance Sampling)
- Resampled Importance Sampling
参考文献
- [Dutré2006] Philip Dutré, Kavita Bala, Philippe Bekaert - "Advanced Global Illumination", 2006, A K Peters