PHPでモンテカルロ法を使って円周率を求める

(1/1)
シミュレーションの一種であるモンテカルロ法を使って円周率を推計するプログラムを PHP を使って作ってみる。同時に、モンテカルロ法の様子を画像として表示できるようにした。

サンプル・プログラム:実行例

PHPでモンテカルロ法を使って円周率を求める

目次

サンプル・プログラム

圧縮ファイルの内容
montecarlo.phpサンプル・プログラム本体

モンテカルロ法とは

ある事象をモデル化した数式や関数があるとき、その定義域に含まれる値をランダムに大量生成して実際に計算し、得られた結果を統計処理することで、推計値が求まるという一種のシミュレーションである。
中性子の運動を探るためにアメリカの数学者スタニスワフ・ウラムが考案し、コンピュータの父フォン・ノイマンが、カジノで有名なモナコ公国の地区の名前をとって命名した。
PHPでモンテカルロ法を使って円周率を求める
モンテカルロ法で円周率を求めることができる。
正方形の紙に円が描かれていたとする。ここで、正方形の中にランダムに点を打ち、円の中に入っているかどうかで場合分けるとする。
PHPでモンテカルロ法を使って円周率を求める
一方、正方形の面積 As および 円の面積 Ac
 mimetex ‥‥(1)
 mimetex ‥‥(2)
である。
正方形の中にランダムに打った点のうち、円の中に入った個数 Nin と入らなかった個数 Nout の比は、面積の比になるから
 mimetex ‥‥(3)
となる。
式 3 に式 1,2 を代入して変形すると、
 mimetex ‥‥(4)
と円周率を求める式となる。

これを PHP プログラムで組んだのが "montecarlo.php" である。計算を簡単にするため、円を 4 分の 1 に切り取って行っている。打つ点の数を 10 万個以上にすると、だいぶ円周率に近くなってくる。

準備:初期値

0023: //リリース・フラグ(公開時にはTRUEにすること)
0024: define('FLAG_RELEASE', FALSE);
0025: 
0026: //リファラ・チェック(直リン防止用;空文字ならチェックしない)
0027: define('REFER_ON', 'www.pahoo.org');
0028: //define('REFER_ON', '');
0029: 
0030: //TrueTypeフォント・ファイル;各自の環境に合わせて要変更
0031: define('FILE_FONT', '../../../../common/font/ipagp.ttf');
0032: 
0033: //初期値
0034: define('COUNT',     10000);          //試行回数
0035: define('WIDTH',     500);            //直径(ドット)
0036: define('COLOR_IN', 'FF0000');        //描画カラー:円周の内側
0037: define('COLOR_OUT', '0000FF');        //描画カラー:円周の外側
0038: define('COLOR_STR', '00FF00');        //描画カラー:円周率
0039: define('SIZE_STR', 30);         //描画サイズ:円周率

上記の初期値は任意に変更が可能である。
FILE_FONT には、円周率の値を描画するために必要となる TrueType フォント・ファイルを、プログラムを配置したパスから見た相対パスで記述する。
描画カラーは、RGB の 16進表記である。

ランダムな点座標

0111: /**
0112:  * ランダムな点座標を1つ返す
0113:  * @param なし
0114:  * @return array(float, float) 座標
0115: */
0116: function getRandomDot() {
0117:     $x = mt_rand() / mt_getrandmax();
0118:     $y = mt_rand() / mt_getrandmax();
0119: 
0120:     return array($x$y);
0121: }

ユーザー関数 getRandomDot は、ランダムな点座標(XY 座標)を 1 つ返す。組み込み関数  mt_rand  と  mt_getrandmax  を利用し、0 以上 1 未満の小数乱数を発生させている。

ポイントを返す

0123: /**
0124:  * ポイントを返す
0125:  * @param $image = 画像リソースID(省略時はピクセル描画しない)
0126:  * @return int ポイント
0127: */
0128: function getPoint($image=NULL) {
0129:     list($x$y) = getRandomDot();
0130:     $res = sqrt(($x * $x) + ($y * $y)) <= 1 ? 1 : 0;
0131: 
0132:     if ($image != NULL) {
0133:         $color = ($res == 1) ? COLOR_IN : COLOR_OUT;
0134:         $arr = color2rgb($color);
0135:         $color = imagecolorallocate($image$arr[0]$arr[1]$arr[2]);
0136:         imagesetpixel($image$x * WIDTHWIDTH - $y * WIDTH$color);
0137:     }
0138:     return $res;
0139: }

ユーザー関数 getPoint は、原点からの距離が 1 以下なら 1 を、1 より大きいなら 0 を返す。
また、引数に画像リソース ID が指定されていれば、キャンパスに点を描画する。

ポイントを集計する

0141: /**
0142:  * ポイントを集計する
0143:  * @param int $n 試行回数
0144:  * @param $image = 画像リソースID(省略時はピクセル描画しない)
0145:  * @return int 集計ポイント
0146: */
0147: function sumPoints($n$image=NULL) {
0148:     $res = 0;
0149:     for ($i = 0; $i < $n$i++) {
0150:         $res += getPoint($image);
0151:     }
0152: 
0153:     return $res;
0154: }

ユーザー関数 sumPoints は、指定した試行回数だけ getPoint を繰り返し実行する。そして、返ってきたポイントを加算し、戻り値とする。

メイン・プログラム

0156: // メイン・プログラム ======================================================
0157: $n  = getParam('n', FALSECOUNT);       //試行回数
0158: 
0159: $image = @imagecreatetruecolor(WIDTHWIDTH);
0160: imagesavealpha($imageTRUE);
0161: 
0162: //モンテカルロ法で円周率を求める
0163: $res = sprintf('%.4f', 4 * sumPoints($n$image) / $n);
0164: 
0165: $arr = color2rgb(COLOR_STR);
0166: $color = imagecolorallocate($image$arr[0]$arr[1]$arr[2]);
0167: $font = dirname(__FILE__) . '/' . FILE_FONT;  //パスの通し方に注意
0168: imagettftext($imageSIZE_STR, 0, 10, SIZE_STR + 10, $color$font$res);
0169: 
0170: header('Content-type: image/png');       //MIMEPNG
0171: imagepng($image);                        //ブラウザ表示
0172: imagedestroy($image);                    //後処理

メイン・プログラムでは、モンテカルロ法で円周率を求め、変数 $res に格納する。
$res の値はキャンパスで描画するのだが、ここでフォント・ファイルへのパスの通し方に注意が必要だ。

参考サイト

(この項おわり)
header