眺めていられるモンテカルロ法で円周率の近似を求めるシミュレーション|JavaScript

私は書いたコードのログが流れてるところをぼーっと眺めるのが好きです。 シミュレーション途中のグラフや画像が変わっていくのも大好きです。

そんな訳で、前回の眺めていられるモンテカルロ法で円周率の近似を求めるシミュレーションに引き続きぼーっとできるシミュレーションを作りました。 今回はモンテカルロ法で円周率の近似を求めるシミュレーションです。

ブラウザはChromeを推奨!

モンテカルロ法

モンテカルロ法とはシミュレーションや数値計算を乱数を用いて行う手法の総称です。ランダム法とも呼ばれる。 例えば、数値積分や統計学のブートストラップ法もモンテカルロ法の一つとされています。

乱数って何だろう?と疑問に思ったら【Unity道場スペシャル 2017札幌】乱数完全マスターのスライドを読んでください。わかりやすくてかなり面白いです!

円周率の求め方

モンテカルロ法を用いて円周率を求める方法は、円と外接する正方形の面積比を使います。

面積をS、半径をrとすると、Scircle=πr2S_{circle} = \pi r^2 です。 この円に外接する正方形の面積は Ssquare=4r2S_{square} = 4r^2となります。 よって、円と外接する正方形の面積比は Scircle/Ssquare=π/4S_{circle} / S_{square} = \pi / 4 であることがわかりました。

説明

ここで、モンテカルロ法を使います。 外接する正方形の範囲に一様に点をプロットした時、円の中にプロットされる確率は面積比と同様に π/4\pi / 4 になります。 プロットした点の数(N)と円の中の点の数(x)の比から、πは π=4x/N\pi = 4x / N で求まるようになる。

たくさん点をプロットしていけば、πにどんどん近似していくはずです。

円周率を求める方法は他にも、ライプニッツの方法ガウス・ルジャンドルの方法などがあります。

シミュレーション

原理はなんとなく伝わったはずなので、実際に計算するところを眺めて見ましょう。 計算を簡単にするため円と外接する四角の 1/41/4 にプロットします。面積比は変わらず π/4\pi / 4 です。

円の中の点を緑、外の点を白としています。 試行回数が点をプロットした回数に当たる。 大体4万回ぐらいから3.14...と少数2桁が合うようになります。 案外、収束しないもんなんですね。

下のグラフがπの理論値と計算した値の比較になります。

πの演算結果:0
試行回数:0 内側の点の数:0 

コードの解説

単純に計算結果だけ求める場合は以下のようにします。

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);
}

円の内側かどうかはx2+y2<r\sqrt{x^2 + y^2} < r で判定することができます。 内側の点を数えて、計算します。

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を使います。 使い方は、それぞれ公式かこのブログのタグの記事を読んでで見てください。p5jsamCharts

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>&nbsp;内側の点の数:<span id="t_inner">0</span>&nbsp;</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()でフレーム処理を行います。 演算速度を早めるため、1フレームに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();
}

再計算するリセットの関数。

終わりに

少数2桁の精度でも意外と収束に時間がかかるんですね。

参考

 関連記事

 新着 or 更新記事