Lesson 5: モンテカルロ法

Lesson 5-1: モンテカルロ法の例: 円周率の計算

この節では、モンテカルロ法について学びます。モンテカルロ法は、ランダムな数(乱数)を用いて問題を解く数理的な手法で、物理学、工学、金融など幅広い分野で応用されています。

以下は、円周率（π）の値を推定するためにモンテカルロ法を適用するFortranコードの例です。この方法は統計的確率の原理に基づき、ランダムサンプリングを使用して数値結果を導きます。

(lesson5_1.f90)

program main
  implicit none
  integer :: isample, nsample, ncount
  real(8) :: x1,x2,r

  nsample = 10000

  open(30,file="montecarlo.out")
  ncount = 0
  do isample = 1, nsample
     call random_number(x1)
     call random_number(x2)

     r = sqrt(x1**2+x2**2)
     if(r < 1d0)ncount = ncount + 1
     write(30,"(I10,2x,I10,2x,e26.16e3)")isample,ncount,4d0*ncount/dble(isample)
     
  end do
  close(30)

end program main

提供されたFortranコードは、モンテカルロ法を用いてπの値を推定します。これは、単位正方形内でランダムな点を生成し、その中に内接する単位円内の点の比率を求めることによって達成されます。以下はコードの詳細です：

1. implicit none: この文は、変数の明示的な宣言を強制し、コードの明瞭性を向上させ、未宣言の変数に関連する潜在的な問題を防ぎます。
2. 変数の宣言:
   • isample: モンテカルロシミュレーションの現在の繰り返しを表すループ変数。
   • nsample: 生成されるモンテカルロサンプルの数（10,000に設定）。
   • ncount: 単位円内の点の数を追跡するためのカウンタ。
   • x1, x2: 0から1までの間で生成されたランダムな数値を格納する変数。
   • r: 点（x1、x2）から原点までのユークリッド距離。
3. open(30, file="montecarlo.out"): "montecarlo.out" という名前のファイルを書き込み用にユニット番号30で開きます。このファイルにはシミュレーションの結果が保存されます。
4. do isample = 1, nsample: 指定されたサンプル数（nsample）のループを開始します。
5. call random_number(x1) および call random_number(x2): ランダムな数値 x1 および x2 を生成します。両方とも0から1の間です。
6. r = sqrt(x1**2 + x2**2): 点（x1、x2）から原点までのユークリッド距離を計算します。
7. if(r < 1d0) ncount = ncount + 1: 点（x1、x2）が単位円内にあるかどうかを確認します（原点からの距離が1未満）。Trueの場合、カウンタ ncount を増やします。
8. write(30, "(I10,2x,I10,2x,e26.16e3)") isample, ncount, 4d0 * ncount / dble(isample): ファイルにフォーマットされた出力を書き込みます。現在のサンプリングに基づく、現在のシミュレーションによる π の推定値が含まれます。
9. end do: ループを終了します。
10. close(30): 出力ファイルを閉じます。

要約すると、このプログラムは単位正方形内でランダムな点を生成し、単位円内の点の割合に基づいて π を推定するモンテカルロシミュレーションを実施します。シミュレーションの結果は "montecarlo.out" という名前のファイルに記録されます。

以下の図は、前述のコードを使用して推定された π の値をサンプリングポイントの数の関数として視覚的に表現しています。図から明らかなように、サンプリングポイントの数が増加するにつれて、その値は正確な π の値（おおよそ3.14）に収束し、大きなサンプリングの極限での収束を示しています。

Lessono5-2：積分のためのモンテカルロ法

このレッスンでは、関数の積分に対するモンテカルロ法を探求します。積分を考えてみましょう： \begin{align} S = \int_a^b dx f(x). \tag{1} \end{align} ここで、積分の値 \(S\) は、下の図(a)で示された赤く塗られた領域に対応します。重要なのは、この積分の値 \(S\) が図(b)で示された青く塗られた領域とも同等であるということです。具体的には、次のように表現されます： \begin{align} S = \int_a^b dx f(x) = (b-a) \times \bar{f}, \tag{2} \end{align} ここで \(\bar{f}\) は区間 \((a \leq x \leq b)\) 上の関数 \(f(x)\) の平均値です。この観察から、積分 \(S\) は積分範囲内の関数 \(f(x)\) の平均値 \(\bar{f}\) を用いて計算できることが示唆されます。

積分のためのモンテカルロ法は、平均値 \( \bar{f} \) を計算することを含みます。モンテカルロ法の積分の中核アルゴリズムでは、式 (\(1\)) を積分するために区間 \((a \le x \le b) \) 内で多くのランダムな数 \(x_i\) が生成されます。そして、平均値 \( \bar{f} \) は以下のように近似されます： \begin{align} \bar{f} \approx \frac{1}{N} \sum_{i=1}^N f(x_i), \tag{3} \end{align} ここで \(N\) は生成されたランダムな数の総数を表します。したがって、式 (\(1\)) の積分は以下のように近似されます： \begin{align} S = \int_a^b dx \, f(x) = (b-a) \times \bar{f} \approx \frac{b-a}{N} \sum_{i=1}^N f(x_i) \tag{4} \end{align} 区間 \((a \le x \le b)\) 内で生成されたランダムな数の集合 \(\{x_i \}\) を用いて。

モンテカルロ法に慣れるために、以下の積分を評価するFortranプログラムを開発してみましょう： \begin{align} S = \int_0^1 dx \frac{4}{1+x^2}. \tag{5} \end{align} 式 (\(5\)) を解析的に積分すると、\(S = \pi\) となります。

以下のPythonコードは、モンテカルロ法を使用して式 (\(5\)) の値を推定する実装の例です。

(lesson5_2.f90)

program main
  implicit none
  integer :: isample, nsample
  real(8) :: x, f, s

  nsample = 10000

  s = 0d0
  open(30,file="pi_montecarlo.out")
  do isample = 1, nsample
     call random_number(x)
     f = 4d0/(1d0+x**2)
     s = s + f
     write(30,"(I10,2x,e26.16e3)")isample,s/isample
  end do
  close(30)

  
end program main

このFortranプログラムは、πの値を推定するためのモンテカルロシミュレーションを実行します。以下はコードの詳細です：

1. implicit none: この文は、すべての変数の明示的な宣言を強制し、コードの明瞭性を向上させ、未宣言の変数に関連する潜在的な問題を防ぎます。
2. 変数の宣言:
   • isample: モンテカルロシミュレーションの現在の繰り返しを表すループ変数。
   • nsample: 生成されるモンテカルロサンプルの数（10,000に設定）。
   • x: 0から1までの間で生成されるランダムな数値を格納する変数。
   • f: ランダムな数値に基づいて計算される関数値を格納する変数。
   • s: 関数値の合計を累積する変数。
3. nsample = 10000: モンテカルロサンプルの数を10,000に設定します。
4. s = 0d0: 変数 s をゼロに初期化します。この変数はシミュレーション中に関数値の合計を累積するために使用されます。
5. open(30, file="pi_montecarlo.out"): "pi_montecarlo.out" という名前のファイルを書き込み用にユニット番号30で開きます。このファイルにはシミュレーションの結果が保存されます。
6. do isample = 1, nsample: 指定されたサンプル数（nsample）のループを開始します。
7. call random_number(x): ランダムな数値 x を0から1まで生成します。
8. f = 4d0 / (1d0 + x**2): 生成されたランダムな数値 x に基づいて関数値 f を計算します。この関数は、モンテカルロ法を使用してπを推定するためのものです。
9. s = s + f: 関数値を変数 s に累積します。
10. write(30, "(I10,2x,e26.16e3)") isample, s/isample: ファイルにフォーマットされた出力を書き込みます。現在のサンプル番号（isample）と累積された関数値の実行平均が含まれており、これを使用してπを推定します。
11. end do: ループを終了します。
12. close(30): 出力ファイルを閉じます。

要約すると、このプログラムはπの値を推定するためのモンテカルロシミュレーションを実行し、ランダムな点が生成され、それらが特定の数学的な関数への寄与が累積されます。累積値の実行平均を使用してπの値を近似し、結果は "pi_montecarlo.out" という名前のファイルに書き込まれます。

以下の図は、上記のコード (lesson5_2.f90) を用いたπの推定値をサンプリングポイントの数の関数として視覚的に示しています。図から明らかなように、サンプリングポイントの数が増加するにつれて、その値は正確なπの値（おおよそ3.14）に収束しています。大きなサンプリングの極限での収束が確認できます。