モンテカルロ法で面積比による
円周率を求めるシミュレーションをしてみました。
が、同様に、
ビュフォンの針実験という試行でも円周率を求められるとのことを知恵袋で教わる。
dの間隔の罫線を引いて、罫線間隔の半分（d/2）の長さの針を彫り投げる。針が、罫線に触れる数を数える。
1/Π＝罫線に触れた数／針を投げた回数
なんとも不思議な計算法だと思わん？ 平行罫線と直線の針を投げると直線しか出てけぇへんのに、どうして円周率が出てくる？？誰が考えたんやろねぇ（ビュフォンという仏物理学者さん）。
罫線に対して、落ちた針の角度は一様であるはず。
罫線に対する角度を横軸（一様発生）に、縦軸に、罫線間隔Dに対するどの範囲で針が触れるかを考える。
針は中点を基準にして、Dという間隔のＤ／２に基準線をもうけて、これに触れるかどうかで、判定することにする。
両側二辺のどちらかに触れるという判定より、Ｄ／２の片側で検討して、その両側で確率２倍というのが考えやすい。
角度０の場合、罫線上の１ポイント（無限小はゼロ？）。９０度の場合は、0～Ｄ／２のどの位置でも触れるので１００％。
４５度の場合、（Ｄ／２）／√２の、正弦波関数になるのが解ると思う。
面積はS=∫（D/2）sinθdθ（0<=θ<=90度）=D/2
全体の面積は同一次元で、縦＊横＝（D/2）＊パイ／４＝D*パイ/2
だから確率は1/パイ となる。
シミュレーションは、針が落ちた場所の始点（X1)と、角度θを乱数で求め、終点X2は、ｘ1＋sinθ/2で計算し、ｘ1＝0，X2<=0, X2>2を、針が触れたと判定した。
100万回リピート： 3.14192713242594
RANDOMIZE
OPTION ANGLE DEGREES
LET en=0
LET rep = 1000000
FOR i=1 TO rep
LET X1=RND
REM LET Y1=RND
LET S=360*RND
LET x2=x1+0.5*SIN(s)
REM LET y2=y1+0.5*COS(s)
IF x1=0 THEN LET en=en+1
IF x2<=0 THEN LET en=en+1
IF x2>=1 THEN LET en=en+1
REM IF en>1 THEN PRINT i/en
NEXT i
IF en>1 THEN PRINT rep/en
END
100回リピートで、針と罫線との関係を図示してみたけど、おもろうない。
この後、面積法の時の様に、真値との収束具合の図でもかいてみようかな？
PR
