はじめに (対象読者・この記事でわかること)

この記事は、Javaでプログラミングを学び始めた初心者から、Monte Carlo法を利用した数値シミュレーションに興味がある中級者までを対象としています。
本稿を読むことで、モンテカルロ法の基本原理と円周率πを推定するアルゴリズムが理解でき、Javaの標準APIだけでシンプルに実装し、結果をコンソール上に出力できるようになることが目指せます。実装例と共に、よくある落とし穴やパフォーマンス改善のポイントも紹介するので、実務や学習プロジェクトでそのまま活用できるはずです。

前提知識

この記事を読み進める上で、以下の知識があるとスムーズです。
- Java SE 8 以降の基本的な文法（クラス・メソッド・ループ・例外処理）
- 基本的な確率・統計の概念（乱数、期待値）
- コマンドラインでのコンパイルと実行（javac、java）

モンテカルロ法による円周率推定の概要

モンテカルロ法は、確率的サンプリングを利用して数値的に積分や面積を求める手法です。円周率πは、半径1の円が単位正方形に内接する面積比（π/4）として表せます。
1. 単位正方形（座標 (0,0)〜(1,1)）内に乱数点 (x, y) を多数生成する
2. 生成した点が円の内部（x² + y² ≤ 1）に入っているか判定する
3. 円内部に入った点の数を全点数で割り、その比率を 4 倍するとπの近似値になる

この手法はアルゴリズムが極めて単純で、乱数生成と条件判定だけで実装できる点が魅力です。サンプル数が増えるほど誤差は √N のオーダーで減少しますが、実行時間とのトレードオフが必要です。

Javaで実装する具体的手順

ステップ1 乱数生成と点の判定ロジックを作る

まずは java.util.Random（もしくは java.security.SecureRandom）で [0,1) の実数を生成し、円内部かどうかを判定するメソッドを書きます。

Java

import java.util.Random;

public class MonteCarloPi {
    private static final Random RAND = new Random();

    /** 乱数点が単位円の内部か判定 */
    private static boolean isInsideCircle() {
        double x = RAND.nextDouble();   // 0.0 以上 1.0 未満
        double y = RAND.nextDouble();
        return x * x + y * y <= 1.0;
    }
}
  • nextDouble() は均等分布で実数を返すため、偏りのないサンプルが得られます。
  • isInsideCircle はシンプルに条件式だけで判定でき、インライン化による最適化が期待できます。

ステップ2 反復計算とπの推定

次に、サンプル数 iterations を引数に取り、円内部点数をカウントし、π を計算する estimatePi メソッドを実装します。

Java

public static double estimatePi(long iterations) {
    long insideCount = 0L;
    for (long i = 0; i < iterations; i++) {
        if (isInsideCircle()) {
            insideCount++;
        }
    }
    return 4.0 * insideCount / iterations;
}
  • long 型を使用することで、数億回規模のサンプルでもオーバーフローを防げます。
  • ループ内部は極めて軽量なので、JIT コンパイラが自動的にインライン展開やループ展開を行い、実行速度が向上します。

ハマった点やエラー解決

1. 乱数シードの固定忘れ

デバッグ時に結果が毎回変わると原因究明が大変です。Random にシードを与えて固定すると再現性が確保できます。

Java

private static final Random RAND = new Random(12345L); // デバッグ用シード

2. int でサンプル数を管理したため、1億回以上でオーバーフロー

int は約2^31‑1 までしか保持できません。大規模サンプルでは必ず long に変更しましょう。

3. nextDouble() の範囲が [0.0, 1.0) であることを忘れ、1.0 が出た際に x*x + y*y が 2 になるケースがあり、境界条件の判定が不正確に

→ nextDouble() は 1.0 を返さないため実際には問題なしですが、nextFloat() を使うと 1.0 が出る可能性がある点に注意。

解決策

  • デバッグ時はシード固定、リリース時はデフォルトシードで乱数の多様性を保つ。
  • ループ変数・カウンタはすべて long に統一し、型安全を徹底する。
  • 境界条件は <= 1.0 として、丸め誤差が入っても安全側に判定させる。
  • 大規模サンプルを実行したい場合は、マルチスレッド化（ForkJoinPool）や Stream.parallel() を活用し、CPU コア数に合わせて分割計算することが推奨されます。

マルチスレッド化例（簡易版）

Java

import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.RecursiveTask;

public class ParallelMonteCarloPi {
    private static final int THRESHOLD = 10_000_000; // タスク分割基準

    private static class PiTask extends RecursiveTask<Long> {
        private final long n;

        PiTask(long n) { this.n = n; }

        @Override
        protected Long compute() {
            if (n <= THRESHOLD) {
                long count = 0;
                for (long i = 0; i < n; i++) {
                    if (MonteCarloPi.isInsideCircle()) count++;
                }
                return count;
            } else {
                long half = n / 2;
                PiTask left = new PiTask(half);
                PiTask right = new PiTask(n - half);
                left.fork();
                long rightResult = right.compute();
                long leftResult = left.join();
                return leftResult + rightResult;
            }
        }
    }

    public static double parallelEstimatePi(long totalSamples) {
        ForkJoinPool pool = new ForkJoinPool();
        long inside = pool.invoke(new PiTask(totalSamples));
        return 4.0 * inside / totalSamples;
    }
}
  • RecursiveTask でサンプル数を分割し、CPU コア数に応じて自動的に並列実行。
  • THRESHOLD は実機のキャッシュサイズやスレッド数に合わせて調整すると最適性能が得られます。

まとめ

本記事では、モンテカルロ法を用いた円周率πの推定をテーマに、Java の標準ライブラリだけで実装する手順とデバッグ時に遭遇しやすい落とし穴、さらにマルチスレッド化による高速化までを網羅しました。

  • モンテカルロ法の原理と円周率への応用方法を理解
  • シンプルな Java 実装（乱数生成・点判定・π計算）をコード例で習得
  • 実装上の課題（シード管理、型選択、境界条件）とその解決策を把握

これにより、読者は 自分だけの数値シミュレーションツール を作成でき、統計的手法の感覚を身につけることができます。次回は、可視化ライブラリ（JavaFX や Plotly）を組み合わせて結果をグラフ化する方法や、他の確率問題への応用例を紹介する予定です。

参考資料