モンテカルロ法で円周率をPythonで求める実装入門

モンテカルロ法で円周率をPythonで求める仕組みと実装

Pythonでモンテカルロ法を10万回試行しても、円周率の誤差が0.01以下になるとは限りません。

モンテカルロ法の基本原理と円周率との関係

モンテカルロ法とは、乱数を使ったシミュレーションによって数値的な答えを求める手法です。名前の由来はモナコのカジノで有名な「モンテカルロ」であり、賭け事に使われる確率論から着想を得ています。

円周率（π）をモンテカルロ法で求める場合、考え方はシンプルです。一辺が1の正方形の中に、半径1の四分円を描きます。この正方形の中にランダムに点を打ち続けたとき、四分円の内側に入った点の割合は「π/4」に近づきます。

つまり、次の式が成り立ちます。

π ≈ 4 × (四分円の内側に入った点の数) ÷ (全体の点の数)

これが基本です。

点が円の内側にあるかどうかは、原点からの距離が1以下かどうかで判定します。座標 (x, y) に対して $$x^2 + y^2 \leq 1$$ を満たせば内側、満たさなければ外側です。この判定を何万回・何十万回と繰り返すことで、円周率の近似値が得られます。

試行回数が多いほど精度は上がります。ただし収束の速さは「試行回数の平方根」に比例するため、精度を10倍にしたければ試行回数を100倍にする必要があります。これは覚えておけばOKです。

金属加工の現場でいえば、金型の断面積計算や円形部品の面積近似など、幾何学的な計算が日常的に発生します。そうした場面でモンテカルロ法の「乱数で面積を推定する」という発想は直感的に理解しやすいはずです。

Pythonでのモンテカルロシミュレーションのコードとrandomモジュール

まずはPython標準ライブラリの `random` モジュールだけを使ったシンプルな実装を見てみましょう。

```python
import random

def monte_carlo_pi(n):
inside = 0
for _ in range(n):
x = random.uniform(0, 1)
y = random.uniform(0, 1)
if x2 + y2 <= 1:
inside += 1
return 4 * inside / n

# 試行回数ごとの結果を比較
for n in 1000, 10000, 100000, 1000000:
pi_estimate = monte_carlo_pi(n)
print(f"試行回数: {n:>8} 円周率の推定値: {pi_estimate:.6f}")
```

このコードを実行すると、試行回数が増えるにつれて 3.14159... に近づく様子が確認できます。これは使えそうです。

`random.uniform(0, 1)` は0以上1以下の一様乱数を返す関数で、ここではx座標とy座標をそれぞれ独立に生成しています。1回のループで1点をプロットし、`x2 + y2 <= 1` の条件を満たしたときだけ `inside` をカウントアップします。

実行環境によって多少異なりますが、試行回数100万回での円周率の推定誤差は平均で ±0.002 程度です。これは小数点以下2桁の精度に相当します。

注意点があります。`random.uniform` はPythonの擬似乱数であるため、シードが同じなら毎回同じ結果が出ます。再現性が必要なデバッグ時は `random.seed(42)` のように固定するとよいでしょう。

Numpyを使った高速なモンテカルロ法の実装と試行回数の比較

`random` モジュールのループ処理は試行回数が増えると処理時間がかかります。そこで実用的には `numpy` を使ったベクトル演算で高速化するのが一般的です。

```python
import numpy as np

def monte_carlo_pi_numpy(n):
x = np.random.uniform(0, 1, n)
y = np.random.uniform(0, 1, n)
inside = np.sum(x2 + y2 <= 1)
return 4 * inside / n

# 速度比較
import time

n = 1_000_000

start = time.time()
monte_carlo_pi(n) # randomモジュール版
print(f"randomモジュール: {time.time() - start:.3f}秒")

start = time.time()
monte_carlo_pi_numpy(n) # numpy版
print(f"numpy版: {time.time() - start:.3f}秒")
```

典型的な環境では、試行回数100万回においてrandomモジュール版が約2〜3秒かかるのに対し、numpy版は0.05秒程度で完了します。約40〜60倍の速度差が出ることになります。

速度が違います。

この速度差の理由は、numpyがCで書かれた内部処理でまとめて演算を行うからです。Pythonのforループは1回ごとに型チェックや変数参照が発生しますが、numpyは配列全体をまとめて処理するため、オーバーヘッドが大幅に削減されます。

試行回数別の精度の目安を以下に示します。

試行回数	推定誤差の目安	所要時間（numpy）
1,000回	±0.05程度	0.001秒未満
10,000回	±0.02程度	0.001秒未満
100,000回	±0.006程度	0.005秒程度
1,000,000回	±0.002程度	0.05秒程度
10,000,000回	±0.0006程度	0.5秒程度

誤差が1/10になるには試行回数を100倍にする必要があります。これが基本です。

モンテカルロ法の結果をmatplotlibで可視化するPythonコード

数字だけでは直感に訴えにくいため、実際に点をプロットして可視化するとモンテカルロ法の仕組みが一目でわかります。

```python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches

def visualize_monte_carlo(n=5000):
x = np.random.uniform(0, 1, n)
y = np.random.uniform(0, 1, n)
inside = x2 + y2 <= 1

pi_estimate = 4 * np.sum(inside) / n

fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(xinside, yinside, color='steelblue', s=1, label='内側')
ax.scatter(x~inside, y~inside, color='salmon', s=1, label='外側')

# 四分円の描画
theta = np.linspace(0, np.pi / 2, 300)
ax.plot(np.cos(theta), np.sin(theta), color='black', linewidth=2)

ax.set_aspect('equal')
ax.set_title(f'モンテカルロ法による円周率推定\n試行回数: {n} 推定値: {pi_estimate:.5f}')
ax.legend(loc='lower left', markerscale=5)
plt.tight_layout()
plt.savefig('monte_carlo_pi.png', dpi=150)
plt.show()

visualize_monte_carlo(5000)
```

青い点が四分円の内側、赤い点が外側です。試行回数5,000回でも、円弧の形がはっきり浮かび上がります。

この可視化コードはJupyter Notebookでも動作します。Google Colabであれば環境構築不要で即座に実行できるため、Pythonに不慣れな方でも試しやすいです。

可視化で見ると理解が深まります。点の密度が均一であるほど推定精度が上がることも、図を見るとよくわかります。

参考として、matplotlibの公式ドキュメントで散布図の詳細なカスタマイズ方法が確認できます。

matplotlib scatter公式ドキュメント（英語）

金属加工現場でのモンテカルロ法活用：公差シミュレーションへの応用

ここが独自視点です。モンテカルロ法は円周率の計算だけでなく、金属加工の品質管理にも直結する手法です。

金属部品の加工では、寸法公差が積み重なることによる「累積公差（スタックアップ）」の問題があります。たとえば、複数の部品を組み合わせるとき、それぞれの部品が公差の上限・下限のどこに収まるかによって、組み合わせ後の寸法が変わります。

最悪ケースを想定した「ワーストケース解析」では公差をすべて足し合わせますが、実際にはすべての部品が同時に最悪値をとる確率は非常に低いです。そこでモンテカルロ法を使うと現実的な分布が得られます。

```python
import numpy as np

# 3つの部品の公差シミュレーション（各部品±0.05mmの正規分布と仮定）
n = 100000
part_a = np.random.normal(loc=10.0, scale=0.025, size=n) # 目標10mm
part_b = np.random.normal(loc=20.0, scale=0.025, size=n) # 目標20mm
part_c = np.random.normal(loc=15.0, scale=0.025, size=n) # 目標15mm

total = part_a + part_b + part_c

print(f"合計寸法の平均: {total.mean():.4f} mm")
print(f"合計寸法の標準偏差: {total.std():.4f} mm")
print(f"45.15mm以上になる確率: {(total > 45.15).mean() * 100:.2f}%")
```

このシミュレーションを使えば、組み立て後の寸法が許容範囲内に収まる確率を事前に推定できます。ワーストケース解析が「±0.15mm」を想定するのに対し、モンテカルロ法では「99%の確率で±0.07mm以内」といった確率的な評価が可能です。

これは実務で使えます。

実際、自動車部品メーカーや航空機部品の製造現場では、設計段階でモンテカルロシミュレーションを使った公差解析が標準的に行われています。SolidWorksやCATIAなどのCADソフトにも公差シミュレーション機能が搭載されており、その内部でもモンテカルロ法が使われているケースがあります。

Pythonでこの手法を理解しておくことで、CADソフトの出力結果を正しく読み解く基礎力が身につきます。公差解析の理解が条件です。

公差のスタックアップ解析についての詳細な解説は、以下の参考リンクでも確認できます。（英語・学術リソース）

Wikipedia: Tolerance analysis（公差解析の概要）

また、Pythonの統計解析ライブラリとして `scipy.stats` を使えば、正規分布以外の分布（一様分布・対数正規分布など）を使ったより高度なシミュレーションも可能です。加工ばらつきが正規分布に従わないと判断した場合に有効です。

SciPy公式ドキュメント: scipy.stats（各種確率分布の一覧）

モンテカルロ法は「答えがわかっているもの（円周率）」で練習してから「答えがわからないもの（公差の累積）」に応用するという学び方が最も効率的です。まずPythonで円周率を求めるコードを動かし、シミュレーションの感覚をつかんでから業務応用に進む流れが推奨されます。

この順番が原則です。