こんにちは、JS2IIUです。
プログラミングの世界では、確定的なロジックを積み上げることが基本ですが、時には「デタラメな数（乱数）」が強力な武器になることがあります。その代表例がモンテカルロ法です。
一見すると無秩序な乱数を使って、円周率（\(pi\)）のような普遍的な定数を近似できるというのは、直感的でありながら非常に数学的な面白さを含んでいます。
本記事では、Pythonと数値計算ライブラリNumPy、そして描画ライブラリMatplotlibを用いて、モンテカルロ法による円周率の計算プロセスを可視化・動画化するプログラムを作成します。単なる計算だけでなく、「どのように値が収束していくか」を目の当たりにすることで、確率的アルゴリズムの本質を理解していきましょう。今回もよろしくお願いします。
モンテカルロ法の基本概念
モンテカルロ法（Monte Carlo method）とは、乱数を用いた試行を繰り返すことで、数値的な解を得る手法の総称です。解析的に（数式変形だけで）解くことが難しい問題でも、シミュレーションを行うことで近似解を得ることができます。
なぜ乱数で円周率がわかるのか？
円周率を求めるために、以下の図形をイメージしてください。
- 一辺の長さが \(1\) の正方形。
- その正方形の中に、半径 \(1\) の扇形（円の4分の1）がぴったり収まっている状態。
ここで、この図形に向かってランダムに無数の点を打ち込みます（ダーツを投げるイメージです）。このとき、数学的に以下の面積の関係が成り立ちます。
- 正方形の面積 \(S_{square} = 1 \times 1 = 1\)
- 扇形の面積 \(S_{quarter} = \pi \times 1^2 \times \frac{1}{4} = \frac{\pi}{4}\)
もし、正方形の中に一様にランダムな点を打ったなら、「扇形の中に入った点の数（\(N_{in}\)）」と「打った点の総数（\(N_{total}\)）」の比率は、面積の比率とほぼ等しくなるはずです。
$$ \frac{N_{in}}{N_{total}} \approx \frac{S_{quarter}}{S_{square}} = \frac{\pi / 4}{1} = \frac{\pi}{4} $$
この式を変形すると、円周率 \(\pi\) を求める式が得られます。
$$ \pi \approx 4 \times \frac{N_{in}}{N_{total}} $$
つまり、「ランダムに点を打ち、扇形に入った割合を4倍する」だけで、円周率が近似できるのです。これが今回実装するアルゴリズムの全貌です。
Pythonによる実装：円周率シミュレーション
それでは、実際にコードを書いていきましょう。今回は計算過程が見えるように、Matplotlibのアニメーション機能を使って、mp4動画（またはGIF）として出力するプログラムを作成します。
以下のコードは、Python環境（Jupyter Notebookやローカル環境）で動作します。ffmpegがインストールされていればmp4で、そうでなければGIFで保存されるように設計しています。
ソースコード
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# --- 設定パラメータ ---
TOTAL_FRAMES = 200 # アニメーションの総フレーム数
POINTS_PER_FRAME = 50 # 1フレームあたりに追加する点の数
TOTAL_POINTS = TOTAL_FRAMES * POINTS_PER_FRAME # 最終的な点の総数
# --- データの初期化 ---
# 点の座標（x, y）を保持するリスト。
# 内側(in)と外側(out)で色を分けるために別々に管理します。
x_inside, y_inside = [], []
x_outside, y_outside = [], []
# 円周率の推定履歴と、その時の試行回数を保持するリスト
pi_history = []
counts_history = []
# カウンター変数の初期化
total_count = 0
inside_count = 0
# --- グラフのセットアップ ---
# 1つのウィンドウに2つのグラフ（散布図と推移図）を並べます
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
fig.suptitle('Monte Carlo Simulation for Pi', fontsize=16)
# 左側のグラフ：散布図（点のプロット）の設定
ax1.set_title("Distribution of Points")
ax1.set_xlim(0, 1)
ax1.set_ylim(0, 1)
ax1.set_aspect('equal') # 縦横比を1:1にして円を歪ませない
# 基準となる円の境界線（半径1の円弧）を描画
theta = np.linspace(0, np.pi/2, 100)
ax1.plot(np.cos(theta), np.sin(theta), 'g--', linewidth=2, label='Boundary')
# プロット用のオブジェクトを初期化（データは空）
scat_in, = ax1.plot([], [], 'r.', markersize=2, label='Inside') # 赤色：円内
scat_out, = ax1.plot([], [], 'b.', markersize=2, label='Outside') # 青色：円外
ax1.legend(loc='upper right')
# 右側のグラフ：収束グラフ（円周率の推定値推移）の設定
ax2.set_title("Convergence of Pi")
ax2.set_xlim(0, TOTAL_POINTS)
ax2.set_ylim(2.5, 4.0) # 円周率3.14が見やすい範囲に設定
# 真の値（3.14159...）を緑の破線で表示
ax2.axhline(y=np.pi, color='g', linestyle='--', label='True Pi')
ax2.set_xlabel("Number of Points")
ax2.set_ylabel("Estimated Pi")
# 折れ線グラフとテキスト表示用のオブジェクト初期化
line_pi, = ax2.plot([], [], 'b-', linewidth=1.5, label='Estimated Pi')
text_pi = ax2.text(0.05, 0.9, '', transform=ax2.transAxes, fontsize=12)
ax2.legend(loc='upper right')
# --- アニメーション更新関数 ---
# この関数がフレームごとに呼び出されます
def update(frame):
global total_count, inside_count
# 1. ランダムな点を生成 (0.0 〜 1.0 の範囲)
# NumPyを使って一括生成することで計算を高速化します
x = np.random.rand(POINTS_PER_FRAME)
y = np.random.rand(POINTS_PER_FRAME)
# 2. 原点(0,0)からの距離の二乗を計算 (x^2 + y^2)
# 距離が1以下なら円の内側です
distances_sq = x**2 + y**2
# 3. 円の内側か外側かを判定するブールマスクを作成
mask_in = distances_sq <= 1.0
# 4. 座標データをリストに追加（描画用）
# mask_inを使って、内側の点と外側の点を振り分けます
x_inside.extend(x[mask_in])
y_inside.extend(y[mask_in])
x_outside.extend(x[~mask_in]) # ~ は否定（NOT）演算子
y_outside.extend(y[~mask_in])
# 5. カウントの更新
count_in_this_frame = np.sum(mask_in)
inside_count += count_in_this_frame
total_count += POINTS_PER_FRAME
# 6. 円周率の推定値を計算: pi = 4 * (内側の数 / 全体数)
current_pi = (inside_count / total_count) * 4 if total_count > 0 else 0
# 履歴データの更新
pi_history.append(current_pi)
counts_history.append(total_count)
# --- グラフ描画オブジェクトの更新 ---
# 左：散布図
scat_in.set_data(x_inside, y_inside)
scat_out.set_data(x_outside, y_outside)
# 右：収束グラフ
line_pi.set_data(counts_history, pi_history)
# テキスト表示（現在の試行回数と推定値）
text_pi.set_text(f'Points: {total_count}\nEstimated Pi: {current_pi:.5f}')
# 更新したオブジェクトを返す（blit=Trueのため）
return scat_in, scat_out, line_pi, text_pi
# --- アニメーション生成と保存 ---
print("アニメーションを作成中...（計算量により時間がかかります）")
# FuncAnimationの作成
ani = FuncAnimation(fig, update, frames=TOTAL_FRAMES, interval=50, blit=True)
# 保存処理：ffmpegがあればmp4、なければgifで保存
output_file = 'monte_carlo_pi.mp4'
try:
# fps=30で滑らかな動画を作成
ani.save(output_file, writer='ffmpeg', fps=30)
print(f"保存完了: {output_file}")
except Exception as e:
print(f"MP4保存エラー: {e}")
print("GIFとして保存を試みます...")
output_file = 'monte_carlo_pi.gif'
try:
ani.save(output_file, writer='pillow', fps=30)
print(f"保存完了: {output_file}")
except Exception as e2:
print(f"GIF保存エラー: {e2}")
# グラフを表示（インタラクティブ環境の場合）
# plt.show()
コードの技術的詳細解説
このプログラムには、Pythonで数値シミュレーションを行う際に役立ついくつかの重要なテクニックが含まれています。
1. NumPyによるベクトル演算（高速化）
モンテカルロ法は「数打ちゃ当たる」戦法なので、試行回数（点の数）が精度に直結します。しかし、Pythonの通常の for ループは計算速度があまり速くありません。
そこで、NumPyのベクトル演算を使用します。
# 遅い書き方（イメージ）
# for i in range(POINTS_PER_FRAME):
# x = random.random()
# ...
# 速い書き方（今回の採用手法）
x = np.random.rand(POINTS_PER_FRAME)
y = np.random.rand(POINTS_PER_FRAME)
distances_sq = x**2 + y**2
mask_in = distances_sq <= 1.0np.random.rand(50) のように書くことで、50個の乱数を一度に生成し、その後の二乗計算や判定も配列全体に対して一括で行われます。これにより、Pythonのループ処理のオーバーヘッドを回避し、高速なシミュレーションが可能になります。
2. ブールインデックス参照
データの振り分けにおいて、以下の記述は非常にPython的で強力です。
x_inside.extend(x[mask_in])
x_outside.extend(x[~mask_in])mask_in は [True, False, True, ...] のような真偽値の配列です。これをインデックスとして渡すことで、True の位置にある要素だけを瞬時に抽出できます。~（チルダ）はビット反転演算子で、ここでは False（円の外側）を抽出するために使用しています。
3. Matplotlibのアニメーション
FuncAnimation は、指定した関数（ここでは update）を繰り返し呼び出し、グラフを少しずつ書き換えることでアニメーションを実現します。
frames:update関数を呼び出す回数。blit=True: グラフ全体を再描画せず、変化があった部分だけを更新する設定。描画速度が向上します。
実行結果と考察
プログラムを実行すると、徐々に点が打たれ、右側のグラフで推定値が変化していく様子が動画として保存されます。
グラフから読み取れること
- 初期段階（試行回数が少ない時）：
右側のグラフを見ると、推定値（青線）が \(3.14\) から大きく外れて、\(3.0\) になったり \(3.4\) になったりと激しく振動します。これはサンプルの偏りが結果に大きく影響しているためです。 - 後半（試行回数が増えた時）：
点が数千、数万と増えるにつれて、振動の幅が小さくなり、真の値（緑の点線）にピタリと吸い寄せられていく様子がわかります。
なぜ収束するのか？（大数の法則）
これは統計学における「大数の法則」によるものです。試行回数を無限に増やしていけば、経験的確率（シミュレーションの結果）は理論的確率（真の面積比）に限りなく近づくことが保証されています。
ただし、モンテカルロ法の収束速度は \(\frac{1}{\sqrt{N}}\) に比例するという特性があります。つまり、精度を1桁上げる（誤差を1/10にする）ためには、計算回数を100倍にする必要があります。これが「計算コストがかかる」と言われる理由です。
モンテカルロ法の応用分野
今回扱った円周率の計算は、あくまで教育的な例題です。しかし、この「ランダムな試行で全体像を把握する」というアプローチは、実社会の極めて重要なシステムで活用されています。
- 金融工学（オプション価格決定）：
株価の動きはランダムな要素を含みます。何万通りもの将来の株価変動シナリオを乱数で生成し、その平均を取ることで、「リスク」や「適正価格」を算出します。 - 強化学習（AlphaGoなど）：
囲碁や将棋のAIでは、現在の局面からランダムに手を打ち進める「プレイアウト」を大量に行い、その勝率が高い手を選択する「モンテカルロ木探索」という手法が使われています。 - 数値積分・物理シミュレーション：
複雑な形状の体積を求めたり、高次元の積分計算を行ったりする場合、規則的に計算するよりもランダムに点を打ったほうが効率が良い場合があります。
まとめ
本記事では、モンテカルロ法を用いて円周率を近似するPythonプログラムを作成し、その仕組みを解説しました。
- モンテカルロ法は、乱数を用いて近似解を求める強力な手法である。
- PythonとNumPyを組み合わせることで、大量の試行を高速に処理できる。
- 可視化することで、試行回数と精度の関係（収束の様子）が直感的に理解できる。
「正確な答えを出すために、あえてデタラメな数を使う」。この逆説的な面白さがモンテカルロ法の魅力です。ぜひ、コードのパラメータ（点の数など）を変更して、収束のスピードがどう変わるか実験してみてください。
最後まで読んでいただきありがとうございました。
コメント