昨日、モンテカルロ法にて円周率を求めるプログラムをつくった。

乱数シュミレーションとしては、まぁまぁの数値を導き出していると個人的には思った。

さらに精度を上げようとして何回か同じ処理をして平均値を取ればいいかな？と思い

次のようにプログラムを変更してみました。

はじめに二重ループを使い平均値を求めるプログラムを組んでしまって

非常に効率が悪いことに気が付いたのは秘密です。(；´Д｀)

相変わらず汚いコードですが素人のやることですからプロフェッショナルな人は

優しく解りやすくご指摘くださればうれしいです。

package jp.yucchi.montecarlo_pi;

import java.lang.management.ManagementFactory;
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.*;
import java.util.logging.Level;
import java.util.logging.Logger;

public class Montecarlo_PI implements Callable> {

    public static void main(String[] args) {
        double TOTAL_THREAD = 1_000_000.0;
        int procs = ManagementFactory.getOperatingSystemMXBean().getAvailableProcessors();
        ExecutorService executor = Executors.newFixedThreadPool(procs);
        List>> futures = new ArrayList<>();
        for (double i = 0; i < TOTAL_THREAD; i++) {
            futures.add(executor.submit(new Montecarlo_PI()));
        }
        executor.shutdown();
        List mPi = new ArrayList<>();
        for (Future> future : futures) {
            try {
                List pi = future.get();
                mPi.addAll(pi);
            } catch (InterruptedException | ExecutionException ex) {
                Logger.getLogger(Montecarlo_PI.class.getName()).log(Level.SEVERE, null, ex);
            }
        }
        double totalPi = 0.0;
        for (Double pi : mPi) {
            totalPi += pi;
        }
        double avgPi = totalPi / mPi.size();
        System.out.println("Ｍａｔｈクラスによる円周率は " + Math.PI);
        System.out.println("モンテカルロ法による円周率は " + avgPi);
    }

    @Override
    public List call() throws Exception {
        List pi = new ArrayList<>();
        double hitCount = 0.0;
        double NEEDLE = 100_000_000.0;
        for (double i = 0; i < NEEDLE; i++) {
            double x = Math.random();
            double y = Math.random();
            if (Math.pow(x, 2.0) + Math.pow(y, 2.0) < 1.0) {
                hitCount++;
            }
        }
        double _pi = 4.0 * hitCount / NEEDLE;
        System.out.println("Montecarlo PI : " + _pi);
        pi.add(_pi);
        return pi;
    }
}

一億個の針を落として計算する処理を百万回して平均値をとってみました。

って・・・計算処理が終わってない(￣。￣;)

しかたないから一億個の針を落とす計算処理を千回に変更してみました。

上記のソースコードの NEEDLE の値を 1000.0 に変更して実行してみました。

これでも終わりそうにないくらい時間がかかります。

不本意ですが計算処理回数を３２回に変更しました。(..;)

実行結果は次のようになりほんの気持ち程度精度は上がりました。

run:
Montecarlo PI : 3.14146332
Montecarlo PI : 3.1416326
Montecarlo PI : 3.1415482
Montecarlo PI : 3.141632
Montecarlo PI : 3.14127944
Montecarlo PI : 3.14178148
Montecarlo PI : 3.1416116
Montecarlo PI : 3.14141588
Montecarlo PI : 3.1416632
Montecarlo PI : 3.14117084
Montecarlo PI : 3.14161116
Montecarlo PI : 3.14168484
Montecarlo PI : 3.14166664
Montecarlo PI : 3.14151984
Montecarlo PI : 3.14131844
Montecarlo PI : 3.14161456
Montecarlo PI : 3.14160796
Montecarlo PI : 3.14157636
Montecarlo PI : 3.14174076
Montecarlo PI : 3.14137984
Montecarlo PI : 3.14161576
Montecarlo PI : 3.14160444
Montecarlo PI : 3.14165568
Montecarlo PI : 3.14140136
Montecarlo PI : 3.1414314
Montecarlo PI : 3.14134932
Montecarlo PI : 3.14164312
Montecarlo PI : 3.14159468
Montecarlo PI : 3.141492
Montecarlo PI : 3.14139636
Montecarlo PI : 3.14149788
Montecarlo PI : 3.1414862
Ｍａｔｈクラスによる円周率は 3.141592653589793
モンテカルロ法による円周率は 3.14153397375
構築成功 (合計時間: 22 分 47 秒)

時間に余裕のあるときに計算処理回数を増やしてみて結果を確認したいと思います。

ちなみにこの計算は Intel Core i7-2600K (3.4GHz) の CPU を積んだハードウェアでしました。

OS は、Windows 7 64bit 版です。

Dual Xeon 4 Core (3.0GHz) のシステムより２倍くらい速いです！

ハイエンドマシンをお持ちの方は是非試してみてくださいね。(ё_ё)

実行結果の合計時間: 22 分 47 秒ってのは使用した NetBeans 7.2 開発版での実行時の値です。

モンテカルロ法という乱数シュミレーションで円周率を導き出してみます。

半径 r の円の面積は PI * r^2 ですよね。

半径 r を 1 とした円を中心から十文字に四等分した右上の部分を使って円周率を求めます。

円の中心から辺の長さ 1 の正方形を重ねます。

正方形の中には４等分された右上の部分の扇形の円がきっちり半径 1 の状態であります。

この時 x 軸、y 軸の 1 の位置は重ねた正方形の辺の長さ 1 になりますね。

よって、正方形の面積 : 扇の面積 = 1 : PI / 4 となります。

したがって PI を求める式は次のようになります。

PI = 4 * 扇の面積/正方形の面積

プログラムは乱数を使ってこの正方形の中にたくさんの針を落とします。

針の落ちた位置が扇の中か外かで円周率を求めていくという力技の方法です。

落とす針の数が多いほど精度は上がる理屈ですが処理時間も増えてしまいます。

とりあえず試してみましょう。 

package jp.yucchi.montecarlo_pi;

import java.util.Random;

public class Montecarlo_PI {

    private static final int NEEDLE = 100000000;

    public static void main(String[] args) {
        int hitCount = 0;
        double x, y, pi;
        for (int i = 0; i < NEEDLE; i++) {
            x = Math.random();
            y = Math.random();
            if (Math.pow(x, 2.0) + Math.pow(y, 2.0) < 1.0) {
                hitCount++;
            }
        }
        pi = 4.0 * hitCount / NEEDLE;
        System.out.println("Ｍａｔｈクラスによる円周率は " + Math.PI);
        System.out.println("モンテカルロ法による円周率は " + pi);
    }
}

実行結果は次のようになりました。

run:
Ｍａｔｈクラスによる円周率は 3.141592653589793
モンテカルロ法による円周率は 3.14185944
構築成功 (合計時間: 9 秒)

何回かやって見たのですが小数点以下３位くらいまでの精度しかでませんでした。

精度を上げるにはこの処理を数回（数万回？）繰り返して平均をとればいいのかなぁ・・・？