suzu6眺めていられるモンテカルロ法で円周率の近似を求めるシミュレーション|JavaScript
私は書いたコードのログが流れてるところをぼーっと眺めるのが好きです。 シミュレーション途中のグラフや画像が変わっていくのも大好きです。
そんな訳で、前回の眺めていられるモンテカルロ法で円周率の近似を求めるシミュレーションに引き続きぼーっとできるシミュレーションを作りました。 今回はモンテカルロ法で円周率の近似を求めるシミュレーションです。
ブラウザはChromeを推奨！
モンテカルロ法
モンテカルロ法とはシミュレーションや数値計算を乱数を用いて行う手法の総称です。ランダム法とも呼ばれる。 例えば、数値積分や統計学のブートストラップ法もモンテカルロ法の一つとされています。
乱数って何だろう？と疑問に思ったら【Unity道場スペシャル 2017札幌】乱数完全マスターのスライドを読んでください。わかりやすくてかなり面白いです！
円周率の求め方
モンテカルロ法を用いて円周率を求める方法は、円と外接する正方形の面積比を使います。
面積をS、半径をrとすると、 です。 この円に外接する正方形の面積は となります。 よって、円と外接する正方形の面積比は であることがわかりました。
ここで、モンテカルロ法を使います。 外接する正方形の範囲に一様に点をプロットした時、円の中にプロットされる確率は面積比と同様に になります。 プロットした点の数(N)と円の中の点の数(x)の比から、πは で求まるようになる。
たくさん点をプロットしていけば、πにどんどん近似していくはずです。
円周率を求める方法は他にも、ライプニッツの方法やガウス・ルジャンドルの方法などがあります。
シミュレーション
原理はなんとなく伝わったはずなので、実際に計算するところを眺めて見ましょう。 計算を簡単にするため円と外接する四角の にプロットします。面積比は変わらず です。
円の中の点を緑、外の点を白としています。 試行回数が点をプロットした回数に当たる。 大体４万回ぐらいから3.14...と少数２桁が合うようになります。 案外、収束しないもんなんですね。
下のグラフがπの理論値と計算した値の比較になります。
コードの解説
単純に計算結果だけ求める場合は以下のようにします。
JavaSript
function calcPi(){
var trials = 10000; // 試行回数
var t_inner = 0; // 円内の点の数
for(var i = 0; i < trials; i++){
var x = Math.random(); // 0 ~ 1.0の乱数
var y = Math.random(); // 0 ~ 1.0の乱数
var r = Math.sqrt(x * x + y * y);
if(r < 1){
// 円の内側
t_inner++;
}
}
// 計算結果
var result = 4 * t_inner / trials;
console.log(result);
}
円の内側かどうかは で判定することができます。 内側の点を数えて、計算します。
Python
import random
inner = 0
for i in range(10000):
x, y = random.random(), random.random()
if x**2 + y**2 < 1:
inner += 1
print(inner *4 / 10000) # 3.1446
コード引用：Remrinのpython攻略日記
シミュレーションのコードについて
基本は上記のアルゴリズムを元に、計算過程を眺められるようグラフなど追加していきました。
利用したライブラリ
点のプロットにはp5js、グラフはamChartsを使います。 使い方は、それぞれ公式かこのブログのタグの記事を読んでで見てください。p5jsとamCharts。
html
<!-- p5js -->
<script src="https://cdnjs.cloudflare.com/ajax/libs/p5.js/0.7.3/p5.min.js"></script>
<!-- amCharts -->
<script src="https://www.amcharts.com/lib/4/core.js"></script>
<script src="https://www.amcharts.com/lib/4/charts.js"></script>
<script src="https://www.amcharts.com/lib/4/themes/material.js"></script>
<script src="https://www.amcharts.com/lib/4/themes/animated.js"></script>
<!-- 計算 -->
<script src="/images/posts/59/montecarlo.js"></script>
<div>
<!-- シミュレーションの表示 -->
<div id="canvas" class="has-text-centered"></div>
<button type="button" class="button is-primary" onclick="restart();">再計算</button>
<div class="m-left-20">πの演算結果：<span id="result">0</span></div>
<div class="m-left-20">試行回数：<span id="trials">0</span> 内側の点の数：<span id="t_inner">0</span> </div>
<div id="chartdiv" class="chart has-text-centered"></div>
</div>
<style>
.chart {
width: 100%;
height: 400px;
}
.m-left-20{
margin-left: 20px;
}
</style>
JavaScript
全文はこちら
var t_inner = 0;
var trials = 0;
// p5js
function setup() {
let canvas = createCanvas(400, 400);
canvas.parent('canvas');
background(50);
}
function draw() {
const width = 400;
while (true) {
monteCarlo(width);
trials++;
// results
if (trials % 50 == 0) {
// 1フレームに50回計算する。
break;
}
}
let pi = t_inner * 4 / trials;
// グラフにデータを追加
addData(pi);
// 結果表示
document.getElementById('result').innerHTML = pi;
document.getElementById('trials').innerHTML = trials;
document.getElementById('t_inner').innerHTML = t_inner;
}
p5jsではsetup()でキャンバスを作成して、draw()でフレーム処理を行います。
演算速度を早めるため、１フレームに50個の点をプロットしています。
function monteCarlo(width) {
var x = Math.random(); // 0 ~ 1.0の乱数
var y = Math.random(); // 0 ~ 1.0の乱数
var r = Math.sqrt(x * x + y * y);
let isInner = false;
if (r < 1) {
isInner = true;
t_inner++;
}
plotEllipse(x * width, y * width, isInner);
}
モンテカルロ法のアルゴリズム部分。 内側の点の数を数えます。
キャンバスの大きさwidthをかけて、点の座標と内側かどうかを渡す。
function plotEllipse(x, y, isInner) {
// color of outer points
let c = color(255);
if (isInner) {
// color of inner points
c = color(25, 254, 0);
}
fill(c);
noStroke();
ellipse(x, y, 1, 1);
}
点を描画する関数。円の内側の点を緑、外側を白とする。
点を描くときは、point()ではなくellipse()を使いましょう。
var series_calc = null;
var series_pi = null;
function createResultChart() {
var chart = am4core.create("chartdiv", am4charts.XYChart);
// X軸
var valueAxisX = chart.xAxes.push(new am4charts.ValueAxis());
valueAxisX.title.text = 'number of trials';
valueAxisX.renderer.minGridDistance = 50;
// スケール固定
valueAxisX.min = 0;
// valueAxisX.max = trials - 1;
var valueAxisY = chart.yAxes.push(new am4charts.ValueAxis());
valueAxisY.title.text = 'π';
// スケール固定
valueAxisY.min = 3.1;
valueAxisY.max = 3.2;
var data = [{ trials: 0, calculated: 0, pi: Math.PI }];
series_calc = createCalculatedSeries(chart);
series_calc.data = data;
series_pi = createPiSeries(chart);
series_pi.data = data;
chart.cursor = new am4charts.XYCursor();
chart.cursor.xAxis = valueAxisX;
chart.legend = new am4charts.Legend();
// chart.scrollbarY = new am4core.Scrollbar();
return;
}
function createCalculatedSeries(chart) {
var series = chart.series.push(new am4charts.LineSeries());
series.name = "calculated π";
series.dataFields.valueY = "calculated";
series.dataFields.valueX = "trials";
series.tooltipText = "{valueY}";
return series;
}
function createPiSeries(chart) {
var series = chart.series.push(new am4charts.LineSeries());
series.name = "theoretical π";
series.dataFields.valueY = "pi";
series.dataFields.valueX = "trials";
series.tooltipText = "{valueY}";
return series;
}
function addData(pi) {
let data = { trials: trials, calculated: pi, pi: Math.PI };
series_calc.addData(data, 0); // 第二引数はpush
series_pi.addData(data, 0); // 第二引数はpush
}
グラフの部分。 グラフの作成とシミュレーション結果をグラフに追加していきます。 大体眺めていられる遺伝的浮動(genetic draft)のシミュレーションと同じ。
addData()はdraw()の中で呼び、フレームを合わせています。
function restart(){
trials = 0;
t_inner = 0;
series_calc = null;
series_pi = null;
chartInit();
setup();
}
再計算するリセットの関数。
終わりに
少数２桁の精度でも意外と収束に時間がかかるんですね。