モンテカルロ積分

通常の数値積分手法（例：台形則）は決定的アルゴリズムだが、モンテカルロ積分は確率的手法であり、各試行によって異なる結果が得られる。最終的な結果は推定値であり、誤差範囲が付随する。

多重積分$${\displaystyle I=\int _{\Omega }f({\overline {\mathbf {x} }})\,d{\overline {\mathbf {x} }}}$$に対し、領域 Ω の体積を$${\displaystyle V=\int _{\Omega }d{\overline {\mathbf {x} }}}$$

とする。単純なモンテカルロ法^[3]ではΩ上で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 N\to \infty }$ で収束する。$${\displaystyle \lim _{N\to \infty }Q_{N}=I.}$$

積分値 Iの近似値として Q_N を用いた場合、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].}$$ これにより、Q_Nの分散は次のように与えられる：$${\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}}.}$$

列$${\displaystyle \left\{\mathrm {E} (\sigma _{1}^{2}),\mathrm {E} (\sigma _{2}^{2}),\mathrm {E} (\sigma _{3}^{2}),\ldots \right\}}$$ は、それぞれがVar(f)に等しいため、有界であるとされ、Var(f)が有限であると仮定すれば、この分散は1/Nに従って漸近的にゼロへと減少する。

したがって、Q_Nの誤差の推定は以下のようになる：$${\displaystyle \delta Q_{N}\approx {\sqrt {\mathrm {Var} (Q_{N})}}=V{\frac {\sqrt {\mathrm {Var} (f)}}{\sqrt {N}}},}$$

この式が示す通り、モンテカルロ積分の誤差は ${\displaystyle {\tfrac {1}{\sqrt {N}}}}$ に比例して減少する。つまり、サンプル数 N を増やせば増やすほど精度が向上する。

この収束速度は積分の次元数に依存しない。この点が、従来の決定的手法と比較した場合のモンテカルロ積分の大きな利点である^[4]。多くの決定的数値積分手法では、次元数の増加に対して計算コストが指数関数的に増加する「次元の呪い」が存在するが、モンテカルロ法ではそれをある程度回避できる。

ただし、モンテカルロ積分では、誤差の見積もりが「確率的な推定値」であり、厳密な誤差上限ではない点に注意が必要である。ランダムサンプリングによって、被積分関数の重要な特徴が十分に捉えられないことがあり、その場合は誤差の過小評価につながる可能性がある。

一様サンプリングによるモンテカルロ法は、高次元で複雑な関数に対してはあまり有効でない。問題に特化したサンプリング分布を用いた改良手法を採用することで初めて決定的手法よりも優れた性能を発揮する。

多くの高次元の被積分関数は、積分領域内のごく狭い領域だけが積分値に大きく寄与する、いわゆる「局在化」した性質を持っている^[5]。このような場合には、適切なサンプル分布を用いることで、必要なサンプル数を大幅に削減しつつ、精度の高い近似が可能となる。

問題に特化したサンプリング技法としては、例えば「層化抽出法」や「重点サンプリング法」などがあり、広く研究されている。

モンテカルロ積分の代表的な例として、円周率 π の推定がある。次の関数を考える：

$${\displaystyle H\left(x,y\right)={\begin{cases}1&{\text{if }}x^{2}+y^{2}\leq 1\\0&{\text{else}}\end{cases}}}$$領域 Ω = [−1, 1] × [−1, 1]、体積 V = 4 とする。このとき、

$${\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})}$$

右図では、相対誤差 ${\displaystyle {\tfrac {Q_{N}-\pi }{\pi }}}$を N の関数として測定しており、${\displaystyle {\tfrac {1}{\sqrt {N}}}}$に比例することがわかる。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
int main() {
    // Initialize the number of counts to 0, setting the total number to be 100000 in the loop.
    int throws = 99999, insideCircle = 0;
    double randX, randY, pi;

    srand(time(NULL));

    // Checks for each random pair of x and y if they are inside circle of radius 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++;
        }
    }

    // Calculating pi and printing.
    pi = 4.0 * insideCircle / throws;
    printf("%lf\n", pi);
}

import numpy as np

rng = np.random.default_rng(0)

throws = 2000
radius = 1

# Choose random X and Y data centered around 0,0
x = rng.uniform(-radius, radius, throws)
y = rng.uniform(-radius, radius, throws)

# Count the times (x, y) is inside the circle,
# which happens when sqrt(x^2 + y^2) <= radius.
inside_circle = np.count_nonzero(np.hypot(x, y) <= radius)

# Calculate area and print; should be closer to Pi with increasing number of throws
area = (2 * radius)**2 * inside_circle / throws
print(area)

再帰的層化抽出法は、一次元の適応求積法（英語版）を多次元積分へと拡張した手法である。

各再帰ステップにおいて、積分値とその誤差は単純なモンテカルロ法によって推定される。もし誤差推定値が指定された精度より大きければ、積分領域は複数の部分領域に分割され、再帰的にこの手順が各部分領域に適用される。

多次元空間においては、単純な「半分に分ける」戦略は適さない。なぜなら、部分領域の数が指数的に増加して追跡が困難になるからである。その代わりに、どの次元で分割すれば最も効果的かを見積もり、その次元のみに沿って領域を分割する。

このアルゴリズムは、関数の分散（ばらつき）が大きい領域にサンプル点を集中させることにより、全体の分散を削減し、サンプリング効率を向上させる。この特徴は、上記の図にも示されている。

代表的な実装例として、類似のアルゴリズムを実装しているMISER法が広く利用されている。

MISER法は、再帰的層化抽出法に基づく手法である。この技法は、分散が最も大きい領域に積分点（サンプル点）を集中させることで、全体の積分誤差を削減することを目的としている^[6]。

層化抽出法の考え方は、以下のような観察に基づいている：2つの互いに交わらない領域 a と b に対して、モンテカルロ積分の推定値がそれぞれ ${\displaystyle E_{a}(f)}$, ${\displaystyle E_{b}(f)}$、対応する分散が ${\displaystyle \sigma _{a}^{2}(f)}$, ${\displaystyle \sigma _{b}^{2}(f)}$ であるとする。このとき、以下のように積分の合成推定値 ${\displaystyle E(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つの部分領域に分割することで進行する。どの軸に沿って分割するかは、d 次元のすべての可能な分割方向を試し、分割後の2領域における分散の合計が最小になる方向を選択する。各部分領域の分散は、そのステップに利用可能なサンプル点の一部を使って推定される。その後、2つの半空間に残りのサンプル点を割り当てる。このサンプル割り当ては前述の比率式に従う。

この再帰的な積分点の割り当ては、ユーザーが指定した再帰の深さまで続けられる。各サブ領域は、単純なモンテカルロ推定によって積分され、最終的にすべてのサブ領域の推定値と誤差が合成されて、全体の積分結果とその誤差見積もりが得られる。

重点サンプリング（Importance Sampling）は、モンテカルロ積分を実行するための非常に重要な手法である^[2][7]。この手法の本質は、領域 Ω 上での一様サンプリングを、より一般的な分布 ${\displaystyle p({\overline {\mathbf {x} }})}$ によるサンプリングへと拡張する点にある。重要度サンプリングでは、積分の推定量 Q_N の分散を小さくするために、関数 f にとって「重要な領域」からサンプルを選ぶ確率を高める分布を選ぶ。

例えば、中心 0、標準偏差 1 の正規分布を持つガウス関数を、範囲 [-1000, 1000] で数値積分したいとする。このとき、一様分布からサンプルを取ると、有効なサンプル（積分に寄与する点）はほとんど得られない。代わりに、同じく中心 0、標準偏差 1 の正規分布に従ってサンプリングすれば、積分値の推定ははるかに効率的になる。この「適切な分布」は被積分関数 f に依存する。

形式的には、以下のようにして推定量を構成する：

与えられたサンプル集合が、$${\displaystyle p({\overline {\mathbf {x} }}):\qquad {\overline {\mathbf {x} }}_{1},\cdots ,{\overline {\mathbf {x} }}_{N}\in V,}$$という分布 ${\displaystyle p({\overline {\mathbf {x} }})}$から得られたとき、積分 I の推定量 Q_N は次のように表される^[2]：$${\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} }})}$ が一定である一様分布の場合にも成り立つ。

実際のサンプル生成には、主にメトロポリス・ヘイスティングス法が使われる^[2]。これは、指定された分布 ${\displaystyle p({\overline {\mathbf {x} }})}$に従って効率的にサンプルを生成でき、積分計算に役立つ。

「ラスベガス法」とは異なります。

VEGASアルゴリズムは、関数 f のヒストグラムを作成し、それに基づいてサンプリング分布を改善することで、目的の分布に近づくよう設計された反復的な重要度サンプリング手法である^[8]。

VEGAS法は以下のように動作する：

まず、積分領域上で複数回パスを行い、被積分関数のヒストグラムを構築する。各パスによって得られたヒストグラムは、次のパスで使用するサンプリング分布を定義する。この操作を繰り返すことで、サンプリング分布が最適化されていく。

VEGASでは、ヒストグラムビン数が次元数 d に対して指数関数的（K^d）に増加するのを避けるために、確率分布を分離可能な関数として近似する：$${\displaystyle g(x_{1},x_{2},\ldots )=g_{1}(x_{1})g_{2}(x_{2})\ldots }$$この近似により、必要なビンの数は Kd に抑えられる。この近似は、関数のピークを各座標軸方向に投影して得るという意味を持つ。

VEGAS法の効率は、この近似の妥当性に依存する。特に、被積分関数のピークが空間的に集中している場合には効率が高くなる。逆に、関数がこの形に書き換えられない（すなわち、軸方向に明確なピークが存在しない）場合には効率が低下する。

VEGAS法は、層化抽出と重点サンプリングの両方の特徴を組み合わせている^[8]。

• モンテカルロ法
• 分散減少法