【Python】NumPyとmatplotlibでモンテカルロ法による円周率推定と収束過程の可視化を試してみた

PythonのNumPyとmatplotlibでモンテカルロ法による円周率推定と収束過程を可視化を題材に、matplotlibで乱数計算と集計を行い、matplotlibで結果を可視化しました。推定値だけでなく、サンプル数を増やしたときの収束、誤差の縮小、複数試行の分布まで実コードで確認します。

隔離した環境で実行した結果は「実際に4枚の図を生成できた」でした。コード、実測値、4枚のグラフを順に見ながら処理を分解します。

このプロンプトの核は「何のグラフを何枚描くか」を先に固定してしまうことです。モンテカルロ法の題材はAIに丸投げすると散布図1枚で終わりがちなので、原理(散布図)→収束(折れ線)→収束速度(両対数)→ばらつき(ヒストグラム)と、見たい観点を4つ明示しました。

再現性のために乱数シードの固定を最初に指定しています。乱数実験は実行のたびに結果が変わるため、シードを固定しないと記事の数値と読者の手元の数値が一致せず、「同じものを作れた」という体験になりません。

また「forループを使わずNumPyの配列演算だけで」と縛ることで、10万点でも一瞬で終わる実装をAIから引き出しています。収束過程の計算にnp.cumsumを名指ししたのも同じ理由で、実装方針まで指定するとコードの質が安定します。

生成プロンプトで工夫したポイント

• np.random.default_rng(42)でシードを固定した乱数生成器を使い：乱数実験は毎回結果が変わるため、シードを固定しないと記事と読者の手元で数値が一致しない。再現性の確保はモンテカルロ系の題材で最重要の指定。
• 点の判定にPythonのforループは使わず、NumPyの配列演算だけで行う：この一文がないとAIは1点ずつ判定するループを書きがち。ベクトル化を強制すると10万点でも一瞬で終わり、NumPyらしいコードが安定して得られる。
• np.cumsumで累積平均を計算し、点数Nを増やしたときの推定値の推移を横軸対数で描き：「収束過程」の計算方法を関数名まで名指しすることで、全Nの推定値を一度に求める効率的な実装に誘導できる。横軸対数の指定で序盤の激しい振れも見やすくなる。
• 絶対誤差|推定値−π|を両対数プロットにし、理論的な収束速度1/√Nに比例する目安線と比較する：両対数にすると1/√Nの減衰が直線に見え、実測誤差が理論通りに減っていることを1枚で示せる。単なる「誤差のグラフ」ではこの比較は出てこない。

同じコードを作れる生成プロンプト（コピペ用・全文）

実際の生成時には、この記事のシステム都合の指定（図の保存先など）を少しだけ足していますが、以下はそのままコピペして同じものが作れるように整えた全文です。

この記事のセットアップ手順と掲載コードは、Windows 11 Pro、PowerShell 5.1、Python 3.13.3で動作確認しています。仮想環境を有効化せず、その中のPythonを直接指定するため、以下のコマンドはPowerShellとコマンドプロンプト（cmd）の両方で使えます。

python -m venv .venv
.\.venv\Scripts\python.exe -m pip install numpy matplotlib

macOS・Linuxでは仮想環境内のPythonパスが異なります。今回の動作確認環境とは異なるため、以下は環境差分を補う参考手順です。

python3 -m venv .venv
./.venv/bin/python -m pip install numpy matplotlib

• グラフのタイトルや軸ラベルを日本語にするため、手元のWindowsで文字化けする場合はplt.rcParams["font.family"] = "Meiryo"のように日本語フォントを指定する。
• 図4は乱数を1,000万点分(2,000回×5,000点)まとめて生成するため150MB程度のメモリを一時的に使う。メモリが厳しい環境では試行回数を減らして調整できる。

NumPyのdefault_rng(42)で一様乱数10万点を一括生成し、x**2 + y**2 <= 1のブール配列だけで円内判定を済ませる(Pythonループなし)。収束過程はnp.cumsumでブール配列の累積和を取り、点数で割って「N点時点の推定値」を全Nについて一度に計算するのがポイント。

可視化は観点を変えた4枚構成——①原理が見える色分け散布図、②対数横軸の収束折れ線、③両対数で1/√Nの理論線と比べる誤差減衰、④2,000回の独立試行による推定値ヒストグラム——で、モンテカルロ法の「収束の遅さ」と「ばらつきの正規分布性」まで見せる。

乱数による推定、サンプル数を増やしたときの収束、誤差の縮小、複数試行の分布を別々の図で確認します。一つの数値だけを見るのではなく、推定が安定していく過程まで視覚化する構成です。

実行に成功したコードは、入力となる乱数の準備、判定と集計、グラフ作成、要約の返却というまとまりに分けています。コード全文に加え、重要な処理だけを実コードから抜き出して確認します。

配列演算を使う箇所とグラフを描く箇所を分離すると、計算結果と表示上の問題を切り分けやすくなります。

"""NumPyとmatplotlibでモンテカルロ法による円周率推定と収束過程を可視化する。"""

from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


def analyze(out_dir: str) -> dict:
    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)

    # 乱数点の生成(シード固定で再現性を確保)
    rng = np.random.default_rng(42)
    n_points = 100_000
    x = rng.random(n_points)
    y = rng.random(n_points)

    # 円周率の推定(単位円の1/4の面積比 = π/4 を利用)
    inside = x**2 + y**2 <= 1.0
    pi_estimate = 4.0 * inside.mean()

    # 図1: 散布図 — ランダムな点を円の内外で塗り分けて原理を見る
    n_show = 5_000
    xs, ys, ins = x[:n_show], y[:n_show], inside[:n_show]
    fig, ax = plt.subplots(figsize=(7, 7))
    ax.scatter(xs[ins], ys[ins], s=4, color="#1f77b4", alpha=0.6, label="円の内側")
    ax.scatter(xs[~ins], ys[~ins], s=4, color="#d62728", alpha=0.6, label="円の外側")
    theta = np.linspace(0.0, np.pi / 2.0, 300)
    ax.plot(np.cos(theta), np.sin(theta), color="black", linewidth=1.5, label="単位円の境界")
    ax.set_title(f"モンテカルロ法の原理(先頭{n_show:,}点を表示)")
    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_aspect("equal")
    ax.legend(loc="upper right")
    fig.tight_layout()
    fig.savefig(out / "scatter_points.png", dpi=120)
    plt.close(fig)

    # 図2: 収束過程 — 点を増やしながら推定値がπに近づく様子(累積平均)
    counts = np.arange(1, n_points + 1)
    running_estimate = 4.0 * np.cumsum(inside) / counts
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.plot(counts, running_estimate, color="#1f77b4", linewidth=1.0, label="推定値の推移")
    ax.axhline(np.pi, color="#d62728", linestyle="--", linewidth=1.2, label="真の円周率 π")
    ax.set_xscale("log")
    ax.set_ylim(2.6, 3.7)
    ax.set_title("円周率推定の収束過程(横軸は対数)")
    ax.set_xlabel("試行点数 N")
    ax.set_ylabel("推定値")
    ax.legend()
    fig.tight_layout()
    fig.savefig(out / "convergence.png", dpi=120)
    plt.close(fig)

    # 図3: 誤差の減衰 — 両対数プロットで理論的な収束速度 1/√N と比較する
    abs_error = np.abs(running_estimate - np.pi)
    checkpoints = np.unique(np.logspace(1, np.log10(n_points), 60).astype(int)) - 1
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.plot(counts[checkpoints], abs_error[checkpoints], marker="o", markersize=4,
            linewidth=1.0, color="#1f77b4", label="実測の絶対誤差")
    ax.plot(counts[checkpoints], 1.7 / np.sqrt(counts[checkpoints]),
            linestyle="--", color="#2ca02c", label="理論的な目安 1/√N に比例")
    ax.set_xscale("log")
    ax.set_yscale("log")
    ax.set_title("推定誤差の減衰(両対数)")
    ax.set_xlabel("試行点数 N")
    ax.set_ylabel("絶対誤差 |推定値 − π|")
    ax.legend()
    fig.tight_layout()
    fig.savefig(out / "error_loglog.png", dpi=120)
    plt.close(fig)

    # 図4: 推定値の分布 — 独立した試行を2,000回繰り返してばらつきを見る
    n_trials, n_per_trial = 2_000, 5_000
    tx = rng.random((n_trials, n_per_trial))
    ty = rng.random((n_trials, n_per_trial))
    trial_estimates = 4.0 * (tx**2 + ty**2 <= 1.0).mean(axis=1)
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.hist(trial_estimates, bins=40, color="#1f77b4", edgecolor="white")
    ax.axvline(np.pi, color="#d62728", linestyle="--", linewidth=1.5, label="真の円周率 π")
    ax.axvline(trial_estimates.mean(), color="#ff7f0e", linewidth=1.5, label="試行の平均")
    ax.set_title(f"推定値の分布({n_per_trial:,}点×{n_trials:,}回の独立試行)")
    ax.set_xlabel("円周率の推定値")
    ax.set_ylabel("回数")
    ax.legend()
    fig.tight_layout()
    fig.savefig(out / "estimate_histogram.png", dpi=120)
    plt.close(fig)

    # 結果の要約を返す
    return {
        "使用した点の数": n_points,
        "円の内側に入った点の数": int(inside.sum()),
        "円周率の推定値": round(float(pi_estimate), 6),
        "真の円周率": round(float(np.pi), 6),
        "絶対誤差": round(float(abs(pi_estimate - np.pi)), 6),
        "N=100時点の推定値": round(float(running_estimate[99]), 6),
        "N=10000時点の推定値": round(float(running_estimate[9_999]), 6),
        "独立試行の平均": round(float(trial_estimates.mean()), 6),
        "独立試行の標準偏差": round(float(trial_estimates.std()), 6),
        "保存した図": [
            "scatter_points.png",
            "convergence.png",
            "error_loglog.png",
            "estimate_histogram.png",
        ],
    }

コード全文は上の折り畳みに入れてあるので、全部を上から読む必要はありません。ここでは特に重要な部分だけを抜き出して、何をしているのか順番に見ていきます。

NumPyで再現可能な乱数を用意する

    # 乱数点の生成(シード固定で再現性を確保)
    rng = np.random.default_rng(42)
    n_points = 100_000
    x = rng.random(n_points)
    y = rng.random(n_points)

乱数生成器とサンプル数を明示し、同じ条件で結果を確かめられるようにしています。

累積和で収束過程をまとめて計算する

    counts = np.arange(1, n_points + 1)
    running_estimate = 4.0 * np.cumsum(inside) / counts
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.plot(counts, running_estimate, color="#1f77b4", linewidth=1.0, label="推定値の推移")
    ax.axhline(np.pi, color="#d62728", linestyle="--", linewidth=1.2, label="真の円周率 π")

各サンプル数での推定値を、Pythonの繰り返しではなくNumPy配列として一括計算しています。

matplotlibで検証用の図を保存する

    fig.tight_layout()
    fig.savefig(out / "scatter_points.png", dpi=120)
    plt.close(fig)

    # 図2: 収束過程 — 点を増やしながら推定値がπに近づく様子(累積平均)

実行結果を後から比較できるよう、作成したグラフを画像ファイルとして保存しています。

複数試行の分布を配列で集計する

    ty = rng.random((n_trials, n_per_trial))
    trial_estimates = 4.0 * (tx**2 + ty**2 <= 1.0).mean(axis=1)
    fig, ax = plt.subplots(figsize=(9, 5))
    ax.hist(trial_estimates, bins=40, color="#1f77b4", edgecolor="white")
    ax.axvline(np.pi, color="#d62728", linestyle="--", linewidth=1.5, label="真の円周率 π")

複数回の試行をまとめて計算し、単一の推定値だけでなく結果のばらつきも確認します。

参考：

©NumPy公式ドキュメント: Random Generator(default_rng)

サンプル数ごとの判定結果を累積し、その時点までの推定値を計算します。横軸を対数にした折れ線にすると、少ない試行で大きく動く段階と、試行を増やして変動が小さくなる段階を同じ図で比べられます。

基準となる値を水平線で重ねることで、推定値がどちら側へずれ、どのように近づくかを読み取れます。

推定値と基準値の差を両対数グラフへ描き、試行数に対して誤差が小さくなる傾向を確認します。さらに同じ条件の推定を複数回実行し、ヒストグラムで結果の中心とばらつきを見ます。

収束グラフは一回の試行の時間変化、ヒストグラムは複数試行の分布を示します。役割の違う図を組み合わせることで、偶然得られた一つの値だけで判断することを避けられます。

隔離した実行環境でコードを動かした結果は「実際に4枚の図を生成できた」でした。掲載しているグラフと表は、その実行で保存された成果物と返り値だけを使っています。

実行時に返された要約は10項目です。記事内の表は可読性のため先頭5項目までとし、4枚の実キャプチャで各グラフを確認できます。

配列の形、保存先、グラフの軸設定を順番に確認すると、計算と描画のどちらで問題が起きたかを絞り込めます。

同じ考え方は、数値積分、確率シミュレーション、リスク評価など、解析的に求めにくい値を試行から近似する場面へ応用できます。

乱数の生成から判定、累積集計、誤差確認、複数試行の分布までを一つのPythonスクリプトで検証しました。役割の異なるグラフを並べることで、推定値だけでなく収束過程とばらつきも確認できます。

この記事の結果は、実行済みコードと「実際に4枚の図を生成できた」という実測結果を根拠にしています。サンプル数やシードを変える場合も、同じ4つの観点で比較すると変化を追いやすくなります。

参考にした一次情報

1. ^ NumPy公式ドキュメント: Random Generator(default_rng). https://numpy.org/doc/stable/reference/random/generator.html, (参照26-07-15).
2. ^ NumPy公式ドキュメント: numpy.cumsum. https://numpy.org/doc/stable/reference/generated/numpy.cumsum.html, (参照26-07-15).
3. ^ matplotlib公式ドキュメント: matplotlib.pyplot.scatter. https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.scatter.html, (参照26-07-15).
4. ^ matplotlib公式ドキュメント: axesのスケール設定(set_xscale/set_yscale). https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.set_xscale.html, (参照26-07-15).

※内容は執筆時点のものです。ライブラリやサイトの仕様は変わる可能性があるため、公式ドキュメントもあわせてご確認ください。

ITやプログラミングに関するコラム

• 
  【CSS】notで複数の件を除外する方法
• 
  【Git】remote設定を変更する方法
• 
  【VBA】コメントアウトを設定する方法
• 
  x86とx64の違いを分かりやすく解説
• 
  マークダウンで改行する方法
• 
  CapsLockキーを解除する方法
• 
  GitLabとGitHubの違いを解説
• 
  Linuxで環境変数を確認する方法
• 
  UbuntuのIPアドレスを確認する方法
• 
  パソコンのメモリの目安を用途別に選ぶ方法

ITやプログラミングに関するニュース

• 
  VercelがAI GatewayにSeedream 5.0 Proを追加、AI SDKのモデル指定で画像生成と編集が可能に
• 
  AWSがAmazon LocationのPlaces APIを強化、住所表記の指定と移動手段別の検索が可能に
• 
  VercelがトレースにTree・Waterfallビューを追加、ログ画面で処理の階層と所要時間を確認可能に
• 
  Googleがエージェント評価の再考を提唱、難易度を情報量で測るDiscovery Benchを解説
• 
  Google CloudがCloud Runサンドボックスを公開プレビューで提供、サービスヘルスは一般提供に
• 
  Google Cloud EMEAが英国金融の重要第三者に指定、イングランド銀行・PRA・FCAの直接監督下に
• 
  AWS DMS Schema ConversionがSQL Serverのオフライン変換に対応、ソースDBへ接続せずスキーマを変換可能に
• 
  EC2 G7インスタンスが米国東部（バージニア北部）で利用可能に、G6比でAI推論性能が最大4.6倍
• 
  SageMaker HyperPodが継続プロビジョニングでのAMIベース構成に対応、S3のスクリプト管理なしでSlurmクラスターを作成可能に
• 
  AWSがEMR on EKSでSparkトラブルシューティングエージェントに対応、失敗ジョブの原因分析を自然言語で依頼可能に