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

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

(2023年1月8日)getRandomDot()改良--PHP8.2対応

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

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

目次

サンプル・プログラム

圧縮ファイルの内容
montecarlo.phpサンプル・プログラム本体
montecarlo.php 更新履歴
バージョン 更新日 内容
1.2.0 2023/01/08 getRandomDot()改良--PHP8.2対応
1.1 2022/04/29 PHP8対応,リファラ・チェック改良
1.0 2019/12/07 初版

モンテカルロ法とは

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

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

準備:初期値

  36: //TrueTypeフォント・ファイル;各自の環境に合わせて要変更
  37: define('FILE_FONT', '../../../../common/font/ipagp.ttf');
  38: 
  39: //初期値
  40: define('COUNT',     100000);        //試行回数
  41: define('WIDTH',     500);           //直径(ドット)
  42: define('COLOR_IN',  'FF0000');      //描画カラー:円周の内側
  43: define('COLOR_OUT', '0000FF');      //描画カラー:円周の外側
  44: define('COLOR_STR', '00FF00');      //描画カラー:円周率
  45: define('SIZE_STR',  30);            //描画サイズ:円周率

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

ランダムな点座標

 130: /**
 131:  * ランダムな点座標を1つ返す
 132:  * @param   なし
 133:  * @return  array(float, float) 座標
 134: */
 135: function getRandomDot() {
 136:     //PHP8.2以上はRandom\Randomizerを使う
 137:     if (isPHPversion(8, 2)) {
 138:         $n = 100000000;
 139:         $random = new Random\Randomizer();
 140:         $x = (double)$random->getInt(0, $n) / $n;
 141:         $y = (double)$random->getInt(0, $n) / $n;
 142:         $random = NULL;
 143:     //PHP8.2未満はmt_randを使う
 144:     } else {
 145:         $x = mt_rand() / mt_getrandmax();
 146:         $y = mt_rand() / mt_getrandmax();
 147:     }
 148: 
 149:     return array($x, $y);
 150: }

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

なお、乱数ジェネレーター問題を解決したPHP8.2以上では、Random\Randomizer クラスを使って乱数を求めるようにした。

ポイントを返す

 152: /**
 153:  * ポイントを返す
 154:  * @param   $image = 画像リソースID(省略時はピクセル描画しない)
 155:  * @return  int ポイント
 156: */
 157: function getPoint($image=NULL) {
 158:     list($x, $y) = getRandomDot();
 159:     $res = sqrt(($x * $x+ ($y * $y)) <1 ? 1 : 0;
 160: 
 161:     if ($image !NULL) {
 162:         $color = ($res == 1? COLOR_IN : COLOR_OUT;
 163:         $arr = color2rgb($color);
 164:         $color = imagecolorallocate($image, $arr[0], $arr[1], $arr[2]);
 165:         imagesetpixel($image, (int)($x * WIDTH), (int)(WIDTH - $y * WIDTH), $color);
 166:     }
 167:     return $res;
 168: }

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

ポイントを集計する

 170: /**
 171:  * ポイントを集計する
 172:  * @param   int $n 試行回数
 173:  * @param   $image = 画像リソースID(省略時はピクセル描画しない)
 174:  * @return  int 集計ポイント
 175: */
 176: function sumPoints($n, $image=NULL) {
 177:     $res = 0;
 178:     for ($i = 0$i < $n$i++) {
 179:         $res +getPoint($image);
 180:     }
 181: 
 182:     return $res;
 183: }

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

メイン・プログラム

 185: // メイン・プログラム ======================================================
 186: $n  = getParam('n', FALSE, COUNT);      //試行回数
 187: 
 188: $image = @imagecreatetruecolor(WIDTH, WIDTH);
 189: imagesavealpha($image, TRUE);
 190: 
 191: //モンテカルロ法で円周率を求める
 192: $res = sprintf('%.4f', 4 * sumPoints($n, $image) / $n);
 193: 
 194: $arr = color2rgb(COLOR_STR);
 195: $color = imagecolorallocate($image, $arr[0], $arr[1], $arr[2]);
 196: $font = dirname(__FILE__. '/' . FILE_FONT;    //パスの通し方に注意
 197: imagettftext($image, SIZE_STR, 0, 10, SIZE_STR + 10, $color, $font, $res);
 198: 
 199: header('Content-type: image/png');      //MIMEはPNGで
 200: imagepng($image);                       //ブラウザ表示
 201: imagedestroy($image);                   //後処理

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

参考サイト

(この項おわり)
header