このページではMathJaxによる数式表示を行っています。対応しているブラウザーであれば少し待てば数式が表示されます。
(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。

モンテカルロ法 (Monte Carlo method)

数値積分の一手法であるモンテカルロ法の基礎の紹介です。
積分等を計算するときに、解析的に計算できればそれが一番良いのですが、それが非常に難しい状況では数値的に計算するしかありません。数値積分の手法には台形公式であったりシンプソン法であったり様々なものが存在しますが、ここで紹介するモンテカルロ法もその一種となります。
CGの計算も基本的には積分を解く問題となっていまして、現在主流のアルゴリズムではモンテカルロ法が使われています。

1次元積分の例

1次元積分
図1. 1次元関数の積分
\begin{equation*} I = \int_{a}^{b} f \Prt{x} dx \end{equation*}

ここでは極々簡単な例として1次元積分の問題を考えます。積分範囲は $ \Brk{a, b} $ となっています。モンテカルロ法の最大の特徴は乱数を使って積分値を求めることで、この積分を次のようにして求めます。

\begin{equation} \Abk{I} = \frac{b - a}{N} \sum_{i = 1}^{N} f \Prt{x_i} \label{eq:MC_uniform} \end{equation}

ここで $ x_i $ は $ \Brk{a, b} $ 中で一様に分布するそれぞれ独立な確率変数、$ N $ はそのサンプル数です。1次元関数の積分は独立変数の座標軸と関数によって囲まれた面積という解釈がありますが、この計算式も、区間中でランダムに発生させたサンプル点の平均値に区間長 $ \Prta{b - a} $ をかけたものと考えれば同じく面積を求めていると言えます。モンテカルロ法はこのように乱数、つまり確率変数を用いて計算を行うという特性上、結果もまた確率変数となります。そのため式\eqref{eq:MC_uniform}の左辺は $ \Abk{I} $ で、あくまで推定値となっています。なお、$ \Abk{I} $ のことをモンテカルロ推定関数と呼びます。この推定関数の期待値が次のように計算されます。

\begin{eqnarray*} E \Brk{\Abk{I}} &=& E \Brk{\frac{b - a}{N} \sum_{i = 1}^N f \Prta{x_i}} = \frac{b - a}{N} \sum_{i = 1}^N E \Brk{f \Prta{x}} \\ &=& \Prta{b - a} \cdot E \Brk{f \Prta{x}} = \Prta{b - a} \cdot \int_a^b f \Prta{x} p \Prta{x} dx \\ &=& \Prta{b - a} \cdot \int_a^b f \Prta{x} \frac{1}{b - a} dx = I \end{eqnarray*}

この期待値が真値に一致することがモンテカルロ法において非常に重要な点となります。

さて、式 $ (1) $ では $ x_i $ を区間中で一様に分布する確率変数と設定していましたが、実は一様分布ではない「任意の」分布(probability density function, PDF)を使用することができます。(個人的に確率の分布を表した関数って確率分布関数と呼びたくなりますが、それだと累積分布関数の意味になってしまうんですよね。ここで言っているのは確率密度関数のことです。今後は、確率密度関数をPDF、累積分布関数をCDFと書くことにします。)その場合のモンテカルロ推定関数は次のようになります。

\begin{equation} \Abk{I} = \frac{1}{N} \sum_{i = 1}^{N} \frac{f \Prt{x_i}}{p \Prt{x_i}} \label{eq:MC} \end{equation}

もちろんこの場合も期待値は真値に一致します。

確率変数である以上、モンテカルロ法による値の推定には必ず分散が付き纏い、式\eqref{eq:MC}の分散は次のようにして計算されます。

\begin{eqnarray*} \sigma^2 \Brk{\langle I \rangle} &=& \sigma^2 \Brk{\frac{1}{N} \sum_{i = 1}^N \frac{f \Prt{x_i}}{p \Prt{x_i}}} = \frac{1}{N^2} \sigma^2 \Brk{\sum_{i = 1}^N \frac{f \Prt{x_i}}{p \Prt{x_i}}} \\ &=& \frac{1}{N^2} \sum_{i = 1}^N \sigma^2 \Brk{\frac{f \Prt{x}}{p \Prt{x}}} = \frac{1}{N} \sigma^2 \Brk{\frac{f \Prt{x}}{p \Prt{x}}} \\ &=& \frac{1}{N} E \Brk{\Prtc{\frac{f \Prt{x}}{p \Prt{x}} - I}^2} \\ &=& \frac{1}{N} \int \Prtc{\frac{f \Prt{x}}{p \Prt{x}} - I}^2 p \Prt{x} dx \end{eqnarray*}

この式から何が言えるかというと、モンテカルロ推定関数の分散がサンプル数 $ N $ に反比例して減少することです。台形公式などは予め設定した区間分割数で精度は決めうちですが、モンテカルロ法はサンプル数を増やせば増やす程、結果が真値に近づき高精度なものとなります。一方で分散が $ N $ に反比例して減少、つまり標準偏差は $ \sqrt{N} $ に反比例して減少するので、結果の誤差はサンプル数を4倍にしても1/2にしかならず、収束がなかなかに遅いという問題があります。

多次元積分への適用

モンテカルロ法は高次元の積分にもそのままの形で適用可能です。以下のような2次元積分を考えます。

\begin{equation*} I = \iint_D f \Prta{x, y} dx dy \end{equation*}

モンテカルロ推定関数は次のように表されます。

\begin{equation*} \Abk{I} = \frac{1}{N} \sum_{i = 1}^{N} \frac{f \Prta{x_i, y_i}}{p \Prta{x_i, y_i}} \end{equation*}

サンプルの各次元の値が独立にサンプルされる場合には次のように表されるでしょう。

\begin{equation*} \Abk{I} = \frac{1}{N} \sum_{i = 1}^{N} \frac{f \Prta{x_i, y_i}} {p \Prta{x_i} \cdot p \Prta{y_i}} \end{equation*}

なぜモンテカルロ法を使うのか

モンテカルロ法の強みとしてそのシンプルさがあります。関数がサンプル可能で、適切なPDFの用意ができれば後はそれらの値を用いて平均値を計算するだけで済みます。そしてもうひとつの強みとして高次元積分に強いことが挙げられます。むしろモンテカルロ法は、今回の例のような低次元の積分には実際は使われないでしょう。台形公式などの決定論的な区分求積法は、各次元に $ N $ サンプル使うとすると、$ d $ 次元積分問題の場合には $ N^d $ 個のサンプルが必要となってしまい、次元が上がるにつれて爆発的にサンプル数が増加してしまいます。これを「次元の呪い(Curse of Dimensionality)」と呼びます。モンテカルロ法は次元数に関係なく任意のサンプル数とすることができ、サンプル数を多くとればとる程正確な結果が得られます。しかしモンテカルロ法には「分散」の概念があり、その収束が遅いという問題があります。幸いなことに分散を抑え収束を高速化する様々な手法が開発されています。モンテカルロ法の分散低減手法」でそれらの一部を紹介します。

参考文献

  • [Dutré2006] Philip Dutré, Kavita Bala, Philippe Bekaert - "Advanced Global Illumination", 2006, A K Peters

Today : 00000071 Total : 00061010
Copyright © 2017 shocker.