前回までは頻度による確率、ラプラスによる確率の定義を確認した上で、両者を結ぶビュッフォンの針と呼ばれる実験について解説しました。今回は Python を使ってビュッフォンの針の実験を実際に行い、円周率を推計します。

ビュッフォンの針 Python

[toc]

2019年4月7日:公開

床に落とした針を特定する

引き続き以下の問題について考えます。

問題(ビュッフォンの針)
床に間隔\(d\)で平行な直線を引き、そこに長さ\(L\)の針を落とす。ただし、\(L<d\)である。このとき針がいずれかの直線と交差する確率は?

床に描かれた平行な線の間隔\(d\)よりも短い長さ\(L\)の針を用意しているため、針が2本以上の直線と交わることはありません。

床に落とした個々の針は、その針の中心から最も近い直線までの最短距離\(x\)(下図のオレンジの線)と、直線に対する針の傾き\(\alpha \)(下図の紫の角度)の2つの情報によって特定できます。つまり、値の組\(\left( \alpha ,x\right) \)が与えられれば床に落とした1本の針を特定できますし、逆に、針を1本落とすと値の組\(\left( \alpha ,x\right) \)が1つだけ得られます。

図:針の位置を特定する情報

では、\(\alpha \)と\(x\)はそれぞれどのような値を取り得るでしょうか。\(x\)は針の中心から最も近い線までの最短距離ですので、針の中心が何らかの直線の上にくる場合に\(x=0\)で最小になり、針の中心が隣り合う2本の直線のちょうど真ん中にくる場合に\(x=\frac{d}{2}\)で最大になります。つまり、\(x\)の変域は、\begin{equation*}
0\leq x\leq \frac{d}{2}
\end{equation*}です。続いて\(\alpha \)の範囲ですが、最小値は針と線が平行の場合の\(\alpha =0\)であり、最大値は針と線が垂直の場合の\(\alpha =\pi \)です。つまり、\(\alpha \)の変域は、\begin{equation*}
0\leq \alpha \leq \frac{\pi }{2}
\end{equation*}です。

以上の情報をPhthonを使ってモデル化しましょう。上と同様に、床に描いた線どうしの間隔は d という変数で、針の中心から最も近い線までの最短距離は x という変数でそれぞれ表現します。また、直線に対する針の角度\(\alpha \)を alpha という変数で表現します。d の値はユーザーが自由に設定する一方で、xalpha はそれぞれランダムに抽出するものとします。ただし、先に解説したように、x は\(0\)以上\(\frac{d}{2}\)以下、alpha は\(0\)以上\(\frac{\pi }{2}\)以下でなければなりません。

具体的には以下の通りです。

ここで利用しているモジュールやメソッドについて簡単に解説します。

random モジュールをインポートすることで利用できる uniform メソッドは指定した範囲に属する値をランダムに生成します。具体的には、\begin{equation*}
\text{random.uniform(a, b)}
\end{equation*}と記述することで、a 以上 b 以下の値を浮動小数点数としてランダムに返します。

math モジュールをインポートすることで利用できる pi メソッドは円周率を返します。具体的には、\begin{equation*}
\text{math.pi}
\end{equation*}と記述すると円周率を得ます。

 

針が床の線と交わることの判定

下図において\(x\)と\(AC\)の長さを比べたとき、\(x\leq AC\)が成り立つ場合には針と直線が交わります。

図:針が線と交わるための条件

では、\(AC\)の長さをどのように求めればよいでしょうか。三角形\(ABC\)の辺である\(AB\)の長さは針の長さ\(L\)の半分\(\frac{L}{2}\)です。さらに、角\(ABC\)の大きさは\(\alpha \)と同じですので、\(\sin \alpha =\frac{AC}{AB}\)という関係から、\begin{eqnarray*}
AC &=&\sin \alpha \times AB \\
&=&\sin \alpha \times \frac{L}{2}
\end{eqnarray*}となります。したがって、針と直線が交わる条件\(x\leq AC\)は、\begin{equation*}
x\leq \sin \alpha \times \frac{L}{2}
\end{equation*}となります。

以上の情報をPhthonを使ってモデル化しましょう。上と同様に、針の長さは L という変数で表現します。L の値もユーザーが自由に設定できますが、その値は d 以下でなければならないことに注意する必要があります。針の長さは床に描いた線どうしの間隔以下ということです。

具体的には以下の通りです。

ここで利用しているモジュールやメソッドについて簡単に解説します。

range 関数はループ処理における繰り返し回数を指定する際に利用します。具体的には、整数 a を用いて\begin{equation*}
\text{range(a)}
\end{equation*}と記述すると、0から始まり1ずつ増加し a 個の整数からなるインデックスを取得します。利用例は以下の通りです。

 

円周率を実験で求める

前回解説したように、ラプラスの確率のもとでは、針と床の線が交わる確率\(p\)は、\begin{equation*}
p=\frac{2L}{d\pi }
\end{equation*}となります。ただし、\(L\)は針の長さ、\(d\)は床に引かれた線の間隔であり、\(L<d\)という関係が成り立ちます。

興味深いのは、上の確率には円周率\(\pi \)が含まれているという点です。そこでこれを円周率について解くと、\begin{equation*}
\pi =\frac{2L}{dp}
\end{equation*}となります。

針の長さ\(L\)と床に引かれた平行線の間隔\(d\)は自分で決められます。さらに、針と床の線が交わる確率\(p\)を頻度による確率として求めます。つまり、針を繰り返し投げて\(p\)を求めるということです。すると上の式の右辺の要素がすべて明らかになるため、結果として円周率\(\pi\)が求まります。つまり、ビュッフォンの針の解を頻度による確率として求めることは、同時に円周率を頻度による確率として求めることにもなります。針を繰り返し投げることで円周率を経験的に求めることができるのです。

Pythonを使って実験を行ってみましょう。

実行結果の例は以下の通りです。

次回は Python を使ってビュッフォンの針の実験を実際に行い円周率を推計します。
次へ進む 演習問題(プレミアム会員限定)