nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

N次元のモンテカルロ積分

一般的なアイデアとして、積分を 体積 x 密度(割合) に焼き直すところがスタート。

\displaystyle
 I = \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \
 = \left( \int_{ \Omega } d \textbf{ r } \right) \
     \left( \frac{ \int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) }{ \int_{ \Omega } d \textbf{ r } } \right) \
 = V_{\Omega} < f >_{\Omega}

密度(割合)は単位量当たりにどれだけ対象となる事物があるかを表し、いわゆる期待値として焼き直せる。連続量の期待値なので、(確率)密度 or 重みをかけた積分に対応。それを離散量で近似しましょうというノリ。このあたりは本当は測度論とかが顔を出すんだろうなぁと妄想していますが、不勉強なので回避。
 \displaystyle
 < f >_{\Omega} = \int_{ \Omega } d \textbf{ r } \; \rho_{ \Omega }( \textbf{ r } ) f( \textbf{ r } ) \
 \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p_{\Omega}( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } )  f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )

 mはサンプリング数またはメッシュの細かさ。トータルのメッシュ数は M \equiv m^N で定義する。
一次元と三次元に関して先ほどの和を解して書いてみると、

\displaystyle
\sum^{ m }_{ \mu } p_{\Omega}( x_{ \mu } ) f( x_{ \mu } ) = p_{\Omega}( x_1 )  f( x_1 ) + p_{\Omega}( x_2 )  f( x_2 ) + \cdots + p_{\Omega}( x_m )  f( x_m ) \\
\displaystyle
\sum^{ m }_{ \mu_x } \sum^{ m }_{ \mu_y } \sum^{ m }_{ \mu_z } p_{\Omega}( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) f( x_{ \mu_x }, y_{ \mu_y }, z_{ \mu_z } ) \\
\displaystyle
\qquad = p_{\Omega}( x_1, y_1, z_1 )  f( x_1, y_1, z_1 ) + p_{\Omega}( x_2, y_1, z_1 )  f( x_2, y_1, z_1 ) \\
\displaystyle
\qquad + \cdots + p_{\Omega}( x_m, y_m, z_{ m-1 } )  f(  x_m, y_m, z_{ m-1 } ) + p_{\Omega} ( x_m, y_m, z_m )  f(  x_m, y_m, z_m )

一次元において、 0 \leq x \leq a における積分を、等間隔のメッシュを切った長方形で近似する場合で書いてみれば、もう少し期待値の描像が明確になる。
 
\displaystyle
\int^{a}_{0} dx \; f(x) \sim \sum^{ m }_{ i = 1 } \frac{ a }{ m } f( ( i - 1 ) a / m ) \\ 
\displaystyle
\int^{a}_{0} dx \; f(x) = a < f >_a \\
\displaystyle
\therefore \; < f >_a \sim \sum^{ m }_{ i = 1 } \frac{ 1 }{ m } f( ( i - 1 ) a / m ) \
 = \sum^{ m }_{ i = 1 } p f( a ( i - 1 ) / m ) \\
\displaystyle
\qquad \sim \sum^{ m }_{ i = 1 } p f( a R_i )

全空間において一様に離散化・サンプリングすると、それは平均値を取りましょうということになる。
そして、これの近似として、 0 < R_i < 1で規格化された昇順の乱数列 \{ R \}を定義して置き換えたものが、積分領域に制約のない場合のモンテカルロ積分となる。
近似関係が分かりやすいようにあえて昇順と言ったが、結局は和なので順番は関係ない。
とりあえず多重積分したいという場合にはこっちの考え方になるだろう。いちいちメッシュの切り方を自分で考えなくて良いので、手軽で良い。
 pが単純にサンプリング数の逆数なので、モンテカルロ積分の例題で有名な \piを求める問題(円の中に入った点の数を数える)は、この考えでは出来ない。

これまでの話は領域 \Omegaの中で閉じていた。
なので、今までの p_{\Omega}はあまり確率の意味はなく、単純に重みの役割が大きかった。
この領域 \Omegaを扱い易い任意の箱 \Boxに拡張するが、その際に余分な領域での値がゼロになるような関数 \Thetaを定義する。
 \displaystyle
\int_{ \Omega } d \textbf{ r } \; f( \textbf{ r } ) \
 = \int_{ \Box } d \textbf{ r } \; \Theta( \textbf{ r } ) f( \textbf{ r } )
 = V_{ \Box } < \Theta f >_{ \Box } \\
\Theta( \textbf{ r } ) = \left\{ \
\begin{array}{ ll }
 1 & ( \textbf{ r } \in \Omega ) \\
 0 & (otherwise)
\end{array}
\right.

これを先ほどと同様に密度関数で表したのちに離散化すると、

 \displaystyle
< \Theta f >_{\Box} = \int_{ \Box } d \textbf{ r } \; \rho_{ \Box }( \textbf{ r } ) \Theta( \textbf{ r } ) f( \textbf{ r } ) \
 = \int_{\Box} d \textbf{ r } \; \rho( \textbf{ r } ) f( \textbf{ r } ) \\
 \displaystyle
 \sim \sum^{ m }_{ \mu_1 } \cdots \sum^{ m }_{ \mu_N } p( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (1) }_{ \mu_N } )  f( r^{ (1) }_{ \mu_1 }, \cdots, r^{ (N) }_{ \mu_N } )

上と同様に一次元の例で考えれば、
 
\displaystyle
\int^{a}_{0} dx \; f(x) = \int^{b}_{0} dx \; \theta( a - x ) f(x) \sim b \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \\ 
\displaystyle
\int^{a}_{0} dx \; f(x) = b < \theta_a f >_b \\
\displaystyle
\therefore \; < \theta_a f >_b \sim \sum^{ m }_{ i = 1 } \frac{ \sigma( ( i - 1 ) b / m ) }{ m } f( ( i - 1 ) b / m ) \
 = \sum^{ m }_{ i = 1 } p( b ( i - 1 ) / m ) f( b ( i - 1 ) / m ) \\
\displaystyle
\qquad \sim \sum^{ m }_{ i = 1 } p( b R_i ) \; f( b R_i ) \\

ここで、

\displaystyle
\sigma( x ) = \left\{
\begin{array}{ll}
 1 & ( x \in [ 0, \; a ] ) \\
 0 & (otherwise)
\end{array}
\right.
と定義した。

 \sigma及び pで、積分領域の調整を乱数点が領域に入った入らないかで行うことが出来る。積分領域が複雑な時には、この方法が威力を発揮することになる。
乱数自体は [ 0, 1 ]の間でしか値が出ないでの、計算したい系に合わせてスケール因子を掛ける必要がある。
一つの次元を c倍する度に結果を c倍すれば正しい積分結果が得られる。