ギャンブルは宇宙の謎を解き明かすか 〜モンテカルロ法で円周率を求める〜

Pocket

今日は数学の話をします。

適当な数字を使ってシミュレーションや数値計算を行う手法があります。モンテカルロ法と言います。カジノで有名なモナコ公国にある地区の名前、モンテ・カルロから名づけられました。今回、このモンテカルロ法を使って円周率を求めてみます。

なお、今回は数学の話ですが、数学が苦手な人でもわかるように頑張って説明してみます。また、実際の計算はプログラミングを行ってコンピュータに行わせますが、こちらについてもプログラミングの知識がなくてもわかってもらえるように頑張ります。

どうやって?

まず、正方形とそれに内接する円を描きます。

この図の青い正方形に向かって「適当に」ボールを投げます。その時に、そのボールが赤い丸の中に入る確率は、面積の比率と一致します。

ボールが赤い丸の中に入る確率 = 赤い丸の面積 ÷ 青い正方形の面積

ボールを何回も投げた場合にボールが赤い丸の中に入る確率は、次の式となります。

ボールが赤い丸の中に入る確率 = ボールを投げて赤い丸に入った回数 ÷ ボールを投げた全部の回数

上の二つの式の左辺は一緒なので、次の式が成り立ちます。

赤い丸の面積 ÷ 青い正方形の面積 = ボールを投げて赤い丸に入った回数 ÷ ボールを投げた全部の回数

次に、説明を簡単にするために図形を4等分します。4等分したうちの右下だけ使うことにします。この時の正方形の1辺の長さを1とします。すると、円の半径も1になります。

この時、それぞれの面積は次の通りです。

赤い円の面積    = 1(半径)× 1(半径)× 円周率 ÷ 4
          = 円周率 ÷ 4

青い正方形形の面積 = 1(縦) × 1(横)
           = 1

切り取っても面積の比率は変わらないので、今までに上げた式を使って円周率を求める式に変形します。

ボールが赤い丸の中に入る確率 = 赤い丸の面積 ÷ 青い正方形の面積
               = 円周率 ÷ 4 ÷ 1
               = 円周率 ÷ 4

ボールが赤い丸の中に入る確率 = ボールを投げて赤い丸に入った回数 ÷ ボールを投げた全部の回数

円周率 ÷ 4         = ボールを投げて赤い丸に入った回数 ÷ ボールを投げた全部の回数

円周率            = ボールを投げて赤い丸に入った回数 ÷ ボールを投げた全部の回数 × 4

ようやく円周率を求める式が完成しました。
あとは、ボールを投げまくって、全部の回数と赤い丸に入った回数を数えれば良いわけです。

で、どうやって?

ここからは、コンピュータの力を借ります。
計算には、ブラウザ(皆さんが今このページを見るために使用しているアプリケーション)に搭載されているJavaScriptというプログラミング言語を使用します。JavaScriptには、0~1の範囲で適当な値を返してくれるMath.random()という命令があります。これを2回呼び出し、それぞれx座標、y座標とします。これを直角三角形の斜辺とみなし、ピタゴラスの定理で斜辺の長さを求めます。

計算式は、次の通りです。

長さ=√(xの2乗+yの2乗)

円の半径は1なので、この長さが1より小さい場合は扇型の中に(x,y)が入っている事になります。

今までの情報をまとめると、こんな感じのプログラムを作れば円周率が求まるはずです。
1.適当なx,yを求める
2.斜辺の長さを求める
3.斜辺の長さが1より小さい場合は「ボールを投げて赤い丸に入った回数」をカウントアップする
4.「ボールを投げた全部の回数」をカウントアップする
5.「ボールを投げて赤い丸に入った回数」÷「ボールを投げた全部の回数」×4を計算して円周率を求める
6.満足するまで1~5を繰り返す

レッツ・トライ

さっそく作ってみました。
開始ボタンを押すと円周率の算出が始まります。もう一度押すと、中断します。
では満足するまで繰り返してください。

お使いのブラウザはJavaScriptが動作しないため、表示できません。

どうでしょう。実際の円周率にそこそこ近い値が算出されたのではないでしょうか。

なんで円周率ちょうどにならないの?

いくつか原因がありますが、原因のひとつとしてMath.random()の精度があります。コンピュータは命令した通りに正確に動くのは得意ですが、「適当に」動くのは苦手なのです。今回使ったMath.random()が、コンピュータにとって苦手な「適当な」作業をさせています。適当な値の事を「乱数」と言うのですが、Math.random()が返す値は「疑似乱数」と言って、文字の通り完璧な乱数ではありません。
・・・いや実は、円周率は無理数であることが証明されているため、Math.random()が完璧な乱数を返却したとしてもこの方法でちょうどの値は求まりません。

ちなみに、上で動いているプログラムのソースコードはこんな感じです。
お見せするのも恥ずかしいものですが・・・

(function() {
//-----------------------------------------------
//-----------------------------------------------
//-----------------------------------------------
var divs = document.getElementsByTagName("div");
var addDiv = null;
for(var i = 0; i < divs.length; i++) {
	if(divs[i].getElementsByTagName("div").length != 0) continue;
	if(divs[i].innerHTML.indexOf("<!-- piMonteCarlo -->") >= 0) {
		addDiv = divs[i];
		break;
	}
}

if(addDiv == null) {
	return;
}

addDiv.style.width = "400px";
addDiv.style.border = "solid 1px";
addDiv.style.padding = "5px";
addDiv.style.backgroundColor = "ivory";
addDiv.innerHTML = "";

var disableSVG = false;
try {
	if (document.implementation.hasFeature("http://www.w3.org/TR/SVG11/feature#Structure", "1.1")) {
		// Enable
	} else {
		disableSVG = true;
	 }
} catch(e) {
	disableSVG = true;
}

if(disableSVG) {
	var spanMessage = document.createElement("span");
	spanMessage.innerHTML = "お使いのブラウザはSVGに対応していないため、動作しません。";
	addDiv.appendChild(spanMessage);
	return;
}

var NS = "http://www.w3.org/2000/svg";
var svg = document.createElementNS(NS, "svg");
svg.setAttribute("width", "200");
svg.setAttribute("height", "200");

var svgBox = document.createElementNS(NS, "polyline");
svgBox.setAttribute("points", "0,0 200,0 200,200 0,200 0,0");
svgBox.setAttribute("fill", "white");
svgBox.setAttribute("stroke", "black");

var fanShape = document.createElementNS(NS, "path");
fanShape.setAttribute("d", "M 0 0 L 0 200 A 200 200 0 0 0 200 0 Z");
fanShape.setAttribute("fill", "none");
fanShape.setAttribute("stroke", "green");

svg.appendChild(svgBox);
svg.appendChild(fanShape);

var piButton = document.createElement("input");
piButton.type = "button";
piButton.value = "開始";

var piButtonClear = document.createElement("input");
piButtonClear.type = "button";
piButtonClear.value = "クリア";
piButtonClear.disabled = true;

var piCheckbox = document.createElement("input");
piCheckbox.type = "checkbox";
piCheckbox.id = "piCheckbox";
piCheckbox.checked = true;

var piLabel = document.createElement("label");
piLabel.htmlFor = "piCheckbox";
piLabel.appendChild(document.createTextNode("描画する"));

var piContinue = false;
var piCount = 0;
var piInCount = 0;

var spanBegin = document.createElement("span");
spanBegin.style.fontSize = "8pt"

var spanPi = document.createElement("span");
spanPi.style.fontSize = "8pt"
spanPi.style.color = "red"

var spanEnd = document.createElement("span");
spanEnd.style.fontSize = "8pt"

var piHead = document.createElement("div");
piHead.appendChild(piButton);
piHead.appendChild(piButtonClear);
piHead.appendChild(spanBegin);
piHead.appendChild(spanPi);
piHead.appendChild(spanEnd);

addDiv.appendChild(piHead);
addDiv.appendChild(svg);
addDiv.appendChild(piCheckbox);
addDiv.appendChild(piLabel);

var blinkCount = 0;
var blinkTimerID = null;

blinkStart();

addEventListener(piButton, "click", function() {
	if(piButton.value == "中断") {
		piButton.value = "再開";
		piContinue = false;
		piButtonClear.disabled = false;
		return;
	}
	piContinue = true;
	piButton.value = "中断";
	piButtonClear.disabled = true;
	blinkStop();

	var timerID = setInterval( function() {
		if(!piContinue) {
			clearInterval(timerID);
			return;
		}
		piCount++;
		
		var randomX = Math.random();
		var randomY = Math.random();

		var result = Math.sqrt(randomX * randomX + randomY * randomY);

		if(result < 1) {
			piInCount++;
		}

		drawInfo();

		if(!piCheckbox.checked) return;

		if(result < 1) {
			var circle = document.createElementNS(NS, "circle");
			circle.setAttribute("cx", randomX * 200);
			circle.setAttribute("cy", randomY * 200);
			circle.setAttribute("r", 2);
			circle.setAttribute("stroke", "red");
			circle.setAttribute("fill", "none");
			circle.setAttribute("stroke-width", 1);
			svg.appendChild(circle);
		} else {
			var line1 = document.createElementNS(NS, "line");
			line1.setAttribute("x1", randomX * 200 - 2);
			line1.setAttribute("y1", randomY * 200 - 2);
			line1.setAttribute("x2", randomX * 200 + 2);
			line1.setAttribute("y2", randomY * 200 + 2);
			line1.setAttribute("stroke", "blue");
			line1.setAttribute("stroke-width", 1);
			svg.appendChild(line1);
			var line2 = document.createElementNS(NS, "line");
			line2.setAttribute("x1", randomX * 200 + 2);
			line2.setAttribute("y1", randomY * 200 - 2);
			line2.setAttribute("x2", randomX * 200 - 2);
			line2.setAttribute("y2", randomY * 200 + 2);
			line2.setAttribute("stroke", "blue");
			line2.setAttribute("stroke-width", 1);
			svg.appendChild(line2);
		}

	}, 1);

});

addEventListener(piButtonClear, "click", function() {
	while (svg.lastChild) {
		svg.removeChild(svg.lastChild);
	}

	piCount = 0;
	piInCount = 0;

	spanBegin.innerHTML = "";
	spanPi.innerHTML = "";
	spanEnd.innerHTML = "";

	svg.appendChild(svgBox);
	svg.appendChild(fanShape);

	piButton.value = "開始";
	piButtonClear.disabled = true;
	blinkStart();
});

function drawInfo() {
	spanBegin.innerHTML = "[ " + piInCount + " / " + piCount + " * 4 = ";
	spanPi.innerHTML = piInCount / piCount * 4;
	spanEnd.innerHTML = " ]";
}

function blinkStart() {
	blinkTimerID = setInterval( function() {
		blinkCount++;
		if(blinkCount < 8) {
			piButton.style.color = "orange";
		} else if(blinkCount <= 8) {
			piButton.style.color = "black";
		} else if(blinkCount >= 10) {
			blinkCount=0;
		}
	}, 100);
}

function blinkStop() {
	if(blinkTimerID == null) return;
	blinkCount=0;
	piButton.style.color = "black";
	clearInterval(blinkTimerID);
	blinkTimerID = null;
}

function addEventListener(element, type, func) {
    if(element.attachEvent) {
        element.attachEvent("on" + type, func);
    } else {
        element.addEventListener(type, func, true);
    }
}
//-----------------------------------------------
//-----------------------------------------------
//-----------------------------------------------
})();

ではでは。

Pocket

コメントを残す

メールアドレスが公開されることはありません。