モンテカルロ法とビュフォンの針、ついでにp5.js

はじめに

　ここでは、ビュフォンの針という針を投げるだけで円周率が求まる不思議な実験についてお教えします。その後、実際に体験することができます。とりあえず体験だけしてみたいという方も、どういう理屈で円周率が求まるのか知りたい方も大歓迎です。専門知識が必要な箇所と不必要な箇所を分けているのでぜひ、ご覧ください。

• この記事の主な対象者
  　・小学生以上の方
  　・数学に興味のある方
• この記事の目標
  　・ビュフォンの針とモンテカルロ法の知識の習得
  　・算数、数学に興味を持ってもらう

ビュフォンの針について

• ビュフォンの針問題の歴史
  　ビュフォンの針は、フランスの博物学者、数学者のジョルジュ=ルイ・ルクレール・ド・ビュフォン伯爵(1707～1788)が提案した、平行線が引かれた平面の上から針(細い棒)をランダムな位置に投げた場合の線と針が交わる確率はどうなるかを求める問題のことです。ここを勘違いしている方も多いと思います。ビュフォンの針自体は、円周率を求める方法ではないんですよね。僕も調べてから初めて知りました。
  
  図1　ビュフォンの針問題のイメージ図
• ビュフォンの針実験について
  　こちらは、先程のビュフォンの針問題を利用した円周率を求める実験法についてになります。ネット上でよく見るビュフォンの針は、ほとんどがこれの事を指しています。具体的な実験法としては、針をランダムに投げて針と線とが交わった回数を求めるという作業を何千回も何万回も繰り返すという手法です。記事の後半で、この実験をシミュレーションするのでそこだけでも寄っていってください。

モンテカルロ法について

• モンテカルロ法の歴史
  　20世紀半ばに、アメリカの数学者のスタニスワフ・ウラムと友人のジョン・フォン・ノイマンによって考案され、賭博で有名な地名モンテカルロの名前をとって「モンテカルロ法」と命名されました。因みにモンテカルロは、フランスの下のモナコ公国にある都市の名前です。
  　モンテカルロ法とは「乱数を取り扱う技法」や「決定論的な数学の問題を乱数を用いて解くこと」の総称です。噛み砕いていうと、ランダムな要素が絡んでくるけれど最終的には一つの結果に到達するような実験の事です。モンテカルロ法の命名自体は、20世紀ですがその歴史的起源は、ビュフォンの針とされています。ビュフォンの針実験も何回もランダムに針を投げますしね。しかし近年では手法の変化によりモンテカルロ法の意味が変わってきており、現在ではコンピュータを用いたシミュレーションや疑似乱数などが主要な部分となっています。

πの導出

• 針と平行線が交わる確率
  　ここからは、数学的な知識を交えてビュフォンの針について説明していきます。積分と確率が分かっていたらなんとなくわかると思います。「そんなことはいいから、早く針投げてるところが見たい！」という方は、「実際にやってみよう‼」 まで飛ばしてください。

結論から言うと針と平行線が交わる確率は、針の長さ＝ℓ、平行線の間隔＝ｄと置き、ℓ≦ｄの場合、

\frac{2×l}{π×ｄ}

となり、図2は針の長さと平行線関係を表したものです。

図2　針の長さと平行線の間隔の関係

　この針と平行線が交わる確率を求めるために、まず針を投げることによって得られる全事象を表します。
　はじめに、投げた針の中心をＰ、点Ｐから最も近い平行線までの距離をｙ、針と平行線とのなす角をθとおきます。図3は、それぞれの関係を表す図の一例です。

図3　針と平行線との関係

➀ 考えうるｙの範囲は、平行線に接している場合から平行線同士のちょうど中心にある場合まで存在するので

０＜＝ｙ＜＝\frac{d}{2}

➁針と平行線のなす角θは、

０＜＝θ＜＝\frac{π}{2}

と言えます。
➀、➁より一般に確率変数ｙ、θは、範囲が０＜＝ｙ＜＝ｄ/2、０＜＝θ＜＝π/2となるようにどの値も同程度にとることができ、一様分布していると言えるので、P(θ、ｙ)のとれる全事象は、図4のようになります。

図4　P(y、θ)の全事象の領域を表したグラフ

次に針と平行線が交わった際に得られる事象とそれらが交わる条件を求めます。
③考えうるｙの範囲は、平行線に接している場合から直角に平行線と交わっている場合まで存在するので

０＜＝ｙ＜＝\frac{d}{2}

④針と平行線のなす角θは、

０＜＝θ＜＝\frac{π}{2}

と言える。
③、④より針と平行線が交わるようなｙの範囲は、

ｙ＜＝\frac{ℓ}{2sinθ}

と表すことができます。そしてこの領域を図4のグラフ上に表すと図４のようになります。
（斜線部：針と平行線が交わる範囲　ドット部：針が落ちる範囲の全事象）

図5　P(y、θ)の全事象の領域と針と平行線が交わる領域(斜線部)を表したグラフ

図5より、針と平行線が交わる確率は、斜線部の面積÷全体の面積で求められます。

図6　針と平行線が交わる確率を求める数式

以上により平行線が引かれた平面の上から針(細い棒)をランダムな位置に投げた場合、針と平行線が交わる確率は、針の長さ＝ℓ、平行線の間隔＝ｄ、ℓ≦ｄの場合

\frac{2ℓ}{πd}

であると言えます。

• πの導出
  　ここまで見て感の良い方は、お気づきかと思いますが先程求めた針と平行線が交わる確率の
  　平行線の間隔が針の2倍、つまりd=2ℓだった場合

\frac{2ℓ}{π×2ℓ}=\frac{1}{π}

となります。これを逆数にするとπが求まります。
　
　つまり、平行線の間隔が針の長さの２倍だった際『「線と交わった針の数」÷「投げた針の総数」の逆数』をすることでπを求めることができます。

• モンテカルロ法を用いてπを求める際の注意点
  　ここまで説明してきたビュフォンの針実験ですが、これにはいくつかの欠点が存在します。主に下記の三つです。
  　・高い精度でπを出そうとすると少なくとも10万回以上針を投げる必要があり、現実的ではないこと。
  　・求められるπの値は、飽くまで近似値であること。
  　・πと求められた値を比較した際、小数点第２位以下の精度が低く求めることができる桁数も10桁程度であること。
  一つ目の問題点については、コンピュータを用いることで大幅な時間の短縮とランダム性の証明をまとめて行うことが可能ですが、二つ目と三つ目はビュフォンの針実験における大きな問題点であると言えます。

使用する言語

　使用する言語は、p5.jsです。この言語は、ProcessingというJAVAを単純化させたプログラミング言語をさらにjavascriptに移植したものです。そのためHTMLやCSSの知識がなくてもプログラムを作成することができ、コード自体も簡略化されているので、非常にソースコードが書きやすいです。また、ブラウザ上で利用することができダウンロードが必要無いため、プログラミング初心者の方にもおすすめです。（完全無料でアカウント登録しなくても利用できる。）

実際にやってみよう‼

• 実際にやってみよう‼
  　図6は、p5.jsにて実際にビュフォンの針実験をシミュレーションしている様子です。このシミュレーションでは、10000本の針を投げ、設定された数のの１０％針を投げる度にその本数とその時点で求まった円周率を表示しています。
  
  図6　実際にシミュレーションをやってみた結果

• ソースコードの解説
  　これがp5.jsで行うビュフォンの針実験のソースコードです。下記の「p5.js」を開いたときに表示されるソースコードを全て消した後、そのままコピー＆ペーストして、左上の実行ボタンを押せば実際に体験することができます。（※実行した後エラーみたいなのが出ますが、仕様なので気にしないでください。）

p5.js

ビュフォンの針

let h=0,l=20,a=0,p=0,c,loopkaisu=0, y_val=40, flag=100,kaisu=100;
//  hit,length,all,pai
function setup() {
  c=400;//実験のキャンパスの大きさ
  createCanvas(1500,1500);
  background("0000");

  line(0,40,400,40);
  line(0,80,400,80);  
  line(0,120,400,120);    
  line(0,160,400,160);   
  line(0,200,400,200);    
  line(0,240,400,240);    
  line(0,280,400,280);    
  line(0,320,400,320);    
  line(0,360,400,360);    
  line(0,400,400,400);  

  line(400,0,400,400);
  line(0,400,400,400);
  textSize(20);
}

function draw() {

  for(j=0;j<10;j++){
    //↑時間短縮のため一回で１０^n本を表示する
      var sx = random(c),ex = random(c),sy = random(c),ey = random(c);
      var val;
      ex = (ex-sx);
      ey = (ey-sy);
      val = sqrt((ex*ex) + (ey*ey));//長さ
    
  line(sx,sy,sx+l*(ex/val),sy+l*(ey/val));
 　   //始点ｘ，始点ｙ，終点ｘ，終点y　（20は適当な長さ）
      //line(0,60,400,60)後はsy,eyを変更するとよい
      a++;
  
    for(i=0;i<=20;i++){
      //線と針が交わっているかを上から順に一本ずつ測定している。（実験用キャンパスの大きさ）÷（針の長さ）＝（iの上限値）
      if(sy <= i*l && i*l <= sy+l*(ey/val))
          h++;
      }
    }
    
  loopkaisu++;
//何本落としたかと求まったπを表示するnoloopで強制終了(一回で１０本落としているので表示する際にflagを１０倍する)
  if(loopkaisu == 1000){
      p=a/h
      textSize(40);
      text(flag*10+"本:"+p,400,y_val); 
    　print(p)
      noloop()
    }
    else if(flag == loopkaisu){
      p=a/h;
      textSize(40);
      text(flag*10+"本:"+p,420,y_val); 
    　print(p)
      y_val+=40;
      flag+=kaisu;
    }
}

おまけ(ビュフォンの針の精度測定)

現在編集中

• どのくらいの精度があるのか
• πを測る際の前提条件
• 実際に計測した結果