複雑なシミュレーションや数値計算を乱数を用いて近似解を求める方法

モンテカルロ法

モンテカルロ法は、複雑な数値計算を、乱数を用いて近似解を求める方法の総称である。

モンテカルロ法を説明するのに、よく円周率などを求める例が使われる。半径 1 の単位円を考える。その単位円の面積は、π × 12 = π と求められる。また、その単位円を囲む正方形の面積は 2 × 2 = 4 である。つまり、正方形の面積と単位円の面積の比は 4 : π である。いま、正方形の中に、ランダムに点を n 個プロットしたとき、単位円の中にプロットされるのは、n × π / 4 個と期待できる。

ここで、R を利用して、-1 ≤ x ≤ 1 および -1 ≤ y ≤ 1 を満たす点を多数生成させて、これらの点のうち x2 + y2 ≤ 1 を満たす点の数をカウントする。単位円の中にプロットされた点の数を利用して、円周率を計算する。

n <- 1e4
x <- runif(n, min = -1, max = 1)
y <- runif(n, min = -1, max = 1)

sum(x * x + y * y <= 1)
## [1] 7852

4 * sum(x * x + y * y <= 1) / n
## [1] 3.1408
モンテカルロ法を利用して円周率を計算する

このような作業を繰り返していくと、円周率を毎回計算することができる。 例えば 10 万回繰り返したときに、各回で計算した π の平均を計算すると、より正確な円周率を求められるようになる。

p <- rep(NA, length = 1e5)
n <- 1e4

for (i in 1:length(p)) {
    x <- runif(n, min = -1, max = 1)
    y <- runif(n, min = -1, max = 1)
    p[i] <- 4 * sum(x * x + y * y <= 1) / n
}

mean(p)
## [1] 3.14158
モンテカルロ法を利用して円周率の分布を計算する