モンテカルロ積分

数学において、モンテカルロ積分は乱数を用いた数値積分法である。モンテカルロ法は定積分を数値的に計算する特殊な方法である。他のアルゴリズムは通常、積分対象を規則的なグリッドで評価するが、^[1]モンテカルロは積分対象を評価する点をランダムに選択する。^[2]この方法は特に高次元積分に有効である。^[3]

モンテカルロ積分を実行するには、均一サンプリング、層別サンプリング、重要度サンプリング、逐次モンテカルロ（粒子フィルターとも呼ばれる）、平均場粒子法など、さまざまな方法があります。

数値積分では、台形則などの方法では決定論的なアプローチが使用されます。一方、モンテカルロ積分では非決定論的なアプローチが採用されており、実現ごとに異なる結果が提供されます。モンテカルロでは、最終結果はそれぞれのエラー バーを含む正しい値の近似値であり、正しい値はそれらのエラー バー内に収まる可能性が高くなります。

モンテカルロ積分が扱う問題は 、Ω（のサブセット）の体積 が $${\displaystyle I=\int _{\Omega }f({\overline {\mathbf {x} }})\,d{\overline {\mathbf {x} }}}$$${\displaystyle \mathbb {R} ^{m}}$$${\displaystyle V=\int _{\Omega }d{\overline {\mathbf {x} }}}$$

単純なモンテカルロ法は、Ω上の点を均一にサンプリングするものである。^[4] N個の均一なサンプル が与えられた場合、$${\displaystyle {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in \Omega ,}$$

私は次のように近似できる $${\displaystyle I\approx Q_{N}\equiv V{\frac {1}{N}}\sum _{i=1}^{N}f({\overline {\mathbf {x} }}_{i})=V\langle f\rangle .}$$

これは大数の法則により、 $${\displaystyle \lim _{N\to \infty }Q_{N}=I.}$$

Q _NからIを推定すると、分散の不偏推定値を使用してサンプル分散によってQ _Nのエラーバーを推定できます 。

$${\displaystyle \mathrm {Var} (f)=\mathrm {E} (\sigma _{N}^{2})\equiv {\frac {1}{N-1}}\sum _{i=1 }^{N}\mathrm {E} \left[\left(f({\overline {\mathbf {x} }}_{i})-\langle f\rangle \right)^{2}\right].}$$ それは $${\displaystyle \mathrm {Var} (Q_{N})={\frac {V^{2}}{N^{2}}}\sum _{i=1}^{N}\mathrm {Var} (f)=V^{2}{\frac {\mathrm {Var} (f)}{N}}=V^{2}{\frac {\mathrm {E} (\sigma _{N}^{2})}{N}}.}$$

シーケンスは Var(f) と等価であるために有界であるため、これが有限であると仮定する限り、この分散は 1/ Nとして漸近的にゼロに減少します。したがって、 Q _Nの誤差の推定値は となり、 これは として減少します。これは、平均の標準誤差にを掛けたものです。この結果は積分の次元数に依存しません。これは、次元に指数関数的に依存するほとんどの決定論的方法に対するモンテカルロ積分の約束された利点です。^[5] 決定論的方法とは異なり、誤差の推定は厳密な誤差境界ではないことに注意することが重要です。ランダムサンプリングでは、誤差の過小評価につながる可能性のある積分関数の重要な特徴がすべて明らかにならない可能性があります。 $${\displaystyle \left\{\mathrm {E} (\sigma _{1}^{2}),\mathrm {E} (\sigma _{2}^{2}),\mathrm {E} (\sigma _{3}^{2}),\ldots \right\}}$$$${\displaystyle \delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sqrt {\mathrm {Var} (f)}}{\sqrt {N}}},}$$${\displaystyle {\tfrac {1}{\sqrt {N}}}}$${\displaystyle V}$

単純なモンテカルロ法は単純な例には有効ですが、決定論的アルゴリズムよりも優れた結果を得るには、問題固有の標本分布を使用するアルゴリズムを使用する必要があります。適切な標本分布を使用すると、高次元の積分対象はほとんどすべて非常に局所的であり、積分に大きく寄与するのは小さな部分空間のみであるという事実を利用できます。^[6] モンテカルロに関する文献の大部分は、誤差推定を改善するための戦略の開発に費やされています。特に、層別サンプリング（領域をサブドメインに分割する）と重要度サンプリング（非一様分布からサンプリングする）は、そのような手法の 2 つの例です。

モンテカルロ積分の典型的な例はπの推定である。関数 と集合Ω = [−1,1] × [−1,1]（V = 4） を考える。$${\displaystyle H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}}$$$${\displaystyle I_{\pi }=\int _{\Omega }H(x,y)dxdy=\pi .}$$

したがって、モンテカルロ積分でπの値を計算する大まかな方法​​は、 Ω上のN個の乱数を選び、 $${\displaystyle Q_{N}=4{\frac {1}{N}}\sum _{i=1}^{N}H(x_{i},y_{i})}$$

右の図では、相対誤差がN の関数として測定されており、 が確認されています。 ${\displaystyle {\tfrac {Q_{N}-\pi }{\pi }}}$${\displaystyle {\tfrac {1}{\sqrt {N}}}}$

真の乱数ジェネレータを使用する必要があることに注意してください。

#include <stdio.h> #include <stdlib.h> #include <time.h> int main () { // カウント数を 0 に初期化し、ループ内の合計数を 100000 に設定します。int throws = 99999 , insideCircle = 0 ; double randX , randY , pi ; 
 
 
  
    
          
       

    srand (時間( NULL ) );

    // x と y のランダムなペアが半径 1 の円の内側にあるかどうかを確認します。
for ( int i = 0 ; i < throws ; i ++ ) { randX = rand () / ( double ) RAND_MAX ; randY = rand () / ( double ) RAND_MAX ; if ( randX * randX + randY * randY < 1 ) { insideCircle ++ ; } }             
             
             
                  
            
        
    

    // 円周率を計算して出力します。
pi = 4.0 * insideCircle / throws ; printf ( "%lf \n " , pi ); }          
     

Pythonで作成されました。

 numpyを npとして インポートする

rng  =  np .ランダム. default_rng ( 0 )

投げる =  2000
半径 =  1


# 0,0 を中心としたランダムな X と Y データを選択しますx  =  rng .uniform ( - radius , radius , throws ) y = rng .uniform ( - radius , radius , throws )  
    

# (x, y) が円の内側にある回数を数えます。
これは sqrt(x^2 + y^2) <= radius のときに発生します
。inside_circle  =  np . count_nonzero ( np . hypot ( x ,  y )  <=  radius )

# 面積を計算して出力します。投げる回数が増えるにつれて円周率に近づくはずです。 
area  =  ( 2  *  radius ) ** 2  *  inside_circle  /  throws 
print ( area )

以下のコードは、Mathematicaでモンテカルロ法を使用して 関数を積分するプロセスを示しています 。 $${\displaystyle f(x)={\frac {1}{1+\sinh(2x)\log(x)^{2}}}}$$${\displaystyle 0.8<x<3}$

func [ x_ ] := 1 / ( 1 + Sinh [ 2 * x ] * ( Log [ x ]) ^ 2 );    

(*収束を早めるために切り捨てられた正規分布からサンプリング*) 
Distrib [ x_ , average_ , var_ ] := PDF [ NormalDistribution [ average , var ], 1.1 * x - 0.1 ]; n = 10 ; RV = RandomVariate [ TruncatedDistribution [{ 0.8 , 3 }, NormalDistribution [ 1 , 0.399 ]], n ]; Int = 1 / n Total [ func [ RV ] / Distrib [ RV , 1 , 0.399 ]] * Integrate [ Distrib [ x , 1 , 0.399 ], { x , 0.8 , 3 }]          
  
      
          

NIntegrate [ func [ x ], { x , 0.8 , 3 }] (*実際の答えと比較*)    

再帰的層別サンプリングは、1 次元の適応求積法を多次元積分に一般化したものです。各再帰ステップで、積分と誤差は単純なモンテ カルロ アルゴリズムを使用して推定されます。誤差推定値が必要な精度より大きい場合、積分ボリュームはサブボリュームに分割され、サブボリュームに手順が再帰的に適用されます。

通常の「2 で割る」戦略は、サブボリュームの数が追跡できないほど急速に増加するため、多次元では機能しません。代わりに、どの次元に沿って細分化すれば最大の利益が得られるかを推定し、その次元に沿ってのみボリュームを細分化します。

層別サンプリング アルゴリズムは、図に示すように、関数の分散が最大となる領域にサンプリング ポイントを集中させ、全体の分散を減らしてサンプリングをより効率的にします。

一般的な MISER ルーチンは同様のアルゴリズムを実装します。

MISERアルゴリズムは再帰的層別サンプリングに基づいています。この手法は、積分点を最も分散の大きい領域に集中させることで、全体的な積分誤差を減らすことを目的としています。^[7]

層別サンプリングの考え方は、積分と分散とをモンテカルロ推定した2つの離れた領域aとbについて、結合推定値の 分散Var( f )が 次のように与えられるという 観察から始まります。${\displaystyle E_{a}(f)}$${\displaystyle E_{b}(f)}$${\displaystyle \sigma _{a}^{2}(f)}$${\displaystyle \sigma _{b}^{2}(f)}$$${\displaystyle E(f)={\tfrac {1}{2}}\left(E_{a}(f)+E_{b}(f)\right)}$$$${\displaystyle \mathrm {Var} (f)={\frac {\sigma _{a}^{2}(f)}{4N_{a}}}+{\frac {\sigma _{b}^{2}(f)}{4N_{b}}}}$$

この分散は、点を次のように分布させることで最小化されることがわかる。 $${\displaystyle {\frac {N_{a}}{N_{a}+N_{b}}}={\frac {\sigma _{a}}{\sigma _{a}+\sigma _{b}}}}$$

したがって、各サブ領域内の関数の標準偏差に比例してサンプル ポイントを割り当てることによって、最小の誤差推定値が得られます。

MISER アルゴリズムは、積分領域を 1 つの座標軸に沿って二分し、各ステップで 2 つのサブ領域を作成します。方向は、すべての可能な二分法を調べ、2 つのサブ領域の合計分散が最小になるものを選択することで選択されます。サブ領域の分散は、現在のステップで使用可能なポイントの総数の一部を使用してサンプリングすることで推定されます。次に、最適な二分法から 2 つの半空間のそれぞれに対して同じ手順が再帰的に繰り返されます。残りのサンプル ポイントは、N _aとN _bの式を使用してサブ領域に割り当てられます。この積分ポイントの再帰的な割り当ては、ユーザーが指定した深度まで続けられ、そこで各サブ領域は単純なモンテ カルロ推定を使用して積分されます。これらの個々の値とそれらの誤差推定は、上方に結合されて、全体的な結果とその誤差の推定値が得られます。

重要度サンプリングアルゴリズムには様々なものがあり、例えば

重要度サンプリングはモンテカルロ積分を実行するための非常に重要なツールを提供します。^{[3] [8]}この方法に対する重要度サンプリングの主な結果は、均一サンプリングがより一般的な選択の特殊なケースであり、サンプルは任意の分布から抽出されるということです。その考え方は、測定値Q _Nの分散を減らすように選択できるというものです。 ${\displaystyle {\overline {\mathbf {x} }}}$${\displaystyle p({\overline {\mathbf {x} }})}$${\displaystyle p({\overline {\mathbf {x} }})}$

0 を中心とし、σ = 1 のガウス関数を -1000 から 1000 まで数値積分する次の例を考えてみましょう。当然、サンプルが区間 [-1000, 1000] で均一に抽出された場合、積分に重要な部分はごく一部になります。サンプルが選択された分布とは異なる分布を選択することで、この状況を改善できます。たとえば、0 を中心とし、σ = 1 のガウス分布に従ってサンプリングします。もちろん、「正しい」選択は積分対象に大きく依存します。

正式には、分布から選択されたサンプルの集合が与えられた場合、 I の推定値は^[3]で与えられる。$${\displaystyle p({\overline {\mathbf {x} }}):\qquad {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in V,}$$ $${\displaystyle Q_{N}\equiv {\frac {1}{N}}\sum _{i=1}^{N}{\frac {f({\overline {\mathbf {x} }}_{i})}{p({\overline {\mathbf {x} }}_{i})}}}$$

直感的に言えば、これは、特定のサンプルを他のサンプルの 2 倍選択した場合、そのサンプルを他のサンプルの半分の重み付けにすることを意味します。この推定値は、が定数である場合の均一サンプリングに対して当然有効です。 ${\displaystyle p({\overline {\mathbf {x} }})}$

メトロポリス・ヘイスティングスアルゴリズムは、からを生成するために最もよく使われるアルゴリズムの1つであり、^[3]積分を計算する効率的な方法を提供します。 ${\displaystyle {\overline {\mathbf {x} }}}$${\displaystyle p({\overline {\mathbf {x} }})}$

VEGAS アルゴリズムは、関数fのヒストグラムを作成する積分領域を何回か通過することで、正確な分布を近似します。各ヒストグラムは、次の通過のサンプリング分布を定義するために使用されます。この手順は、漸近的に目的の分布に収束します。^{[9]ヒストグラムのビンの数が}K ^{d の}ように増加するのを避けるために、確率分布は分離可能な関数で近似されます 。 これにより、必要なビンの数はKdのみになります。これは、被積分関数の座標軸への投影から関数のピークを特定することと同じです。VEGAS の効率は、この仮定の妥当性に依存します。被積分関数のピークが適切に局所化されている場合に最も効率的です。被積分関数を近似的に分離可能な形式で書き直すことができれば、VEGAS による積分の効率が向上します。VEGAS には、いくつかの追加機能が組み込まれており、層別サンプリングと重要度サンプリングの両方が組み合わされています。^[9]$${\displaystyle g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots }$$

• 準モンテカルロ法
• 補助フィールドモンテカルロ
• 統計物理学におけるモンテカルロ法
• モンテカルロ法
• 差異の削減

• カフェ数学 : モンテカルロ積分 : モンテカルロ積分（原理、仮説、信頼区間）を説明するブログ記事
• Boost.Math : ナイーブ モンテ カルロ統合: C++ ナイーブ モンテ カルロ ルーチンのドキュメント
• 統計物理学の問題に応用されたモンテカルロアプレット