社会工学実習
経営工学分野 第１回
– ビュッフォンの針で円周率を –
担当：八森正泰・山本芳嗣・後藤順哉
1
ビュッフォンの針の問題の理論（１時間目）
図 1 のように針を等間隔に引かれた平行線上にでたらめに落として、針が平行線と交差するかどうかを
観測する。この実験はビュッフォン (Buffon) の針の実験と呼ばれている。脚注に引用したように彼は肩越
しに棒切れを投げたと伝えられている。でも、いったに何のためにビュッフォンさん1 はこのようなことを
したのか。針の長さを 、平行線の間隔を L（簡単のために < L）として説明しよう。落とした針の一端
図 1: 平行線と針達
は 2 本の平行線 D1 と D2 にあるとして、針が平行線 D1 と交わる確率 P1 を計算してみよう。平行線と交
わる確率を知りたければそれを 2 倍すればよい。そのため針と平行線とのなす角度を θ、針の指定された端
（たとえば糸通し穴の有る方）と D1 との距離を y とする。まず、角度 θ と y は、それぞれ
−π ≤ θ ≤ π; 0 ≤ y ≤ L
(1.1)
1 Georges Louis Leclerc Comte de Buffon
Born: 7 Sept 1707 in Montbard, C¯
ote d’Or, France Died: 16 April 1788 in Paris, France
At the age of 20 Georges Buffon discovered the binomial theorem. He corresponded with Cramer on mechanics, geometry,
probability, number theory and the differential and integral calculus. His first work Sur le jeu de franc-carreau introduced
differential and integral calculus into probability theory. He next wrote Th?orie de la terre and became the most important
natural historian of his day having great influence across a wide scientific field. He is remembered most in mathematics
for a probability experiment which he carried out calculating π by throwing sticks over his shoulder onto a tiled floor and
counting the number of times the sticks fell across the lines between the tiles. This experiment caused much discussion
among mathematicians which helped towards an understanding of probability.
1
D1
l
θ
lsinθ
y
L
D2
図 2: 平行線と１本の針
の範囲を動く。さらに、
針が D1 と交差する
⇔ y ≤ sin θ
(1.2)
である。従って投げられた針の θ と y が上の矩形の範囲 (1.1) に一様な確率で落ちると仮定すると、それが
D1 と交わる確率はこの矩形と、矩形の中で (1.2) を満たす (θ, y) の作る領域の面積比
π
sin θdθ
0
π
Ldθ
−π
となる。ここで分子の積分範囲が [0, π] であることに注意。図 3 参照。sin の積分が − cos であることを思
い出せば、分子分母はそれぞれ
π
π
sin θdθ =
sin θdθ
0
0
= [− cos θ]π0 = − [cos θ]π0
= − ((−1) − (1)) = 2
π
π
Ldθ = L
−π
1dθ = 2πL
−π
となる。したがって、
P1 =
さらに、針が平行線と交わる確率 P は
P = 2P1 =
2
(1.3)
πL
2
πL
(1.4)
L
lsinθ
−π
π
0
図 3: ２つの領域
となる。(1.4) より、針を N 回投げて n 回平行線上に落ちたときは、n/N はこの確率 P を近似しているは
ずであるから
n
2
≈
N
πL
となる。ここで ≈ は近似を意味する。だから、N と n から円周率 π を次式のように推定することができる。
π≈
2 N
L n
(1.5)
ビュッフォンは、このようにして π の値を実験により求める方法（今日の言葉ではモンテカルロ法) を示唆
した。下の表はいにしえの忍耐強い人々の実験結果を示している。([1], [2])
表 1: ビュッフォンの針の実験による π の推定値
忍耐強い人
N
やった年
n
推定値 2 N/Ln
1850
0.8 5000
2532 3.1596
1855
0.6 3204 1218.5 3.1553
ドゥ・モルガン
1860
1
600
382.5 3.137
フォックス大尉
1864
0.75 1030
489 3.1595
ラッツァリーニ
1901
0.83 3408
1808 3.1415929
レイナ
1925
0.5419 2520
859 3.1795
注：平行線の間隔 L = 1.0、交わっているか曖昧な場合に 0.5 とカウント
ウルフ
スミス・ダベルディーン
ちょっとラッツァリーニの結果は良すぎるような気がするので、思考実験をしてみよう。今 L = 10cm、
= 7.85398cm として、２回針を投げて１回だけ交わったら、π の推定値として
2 N
= 2 × 0.785398 × 2 = 3.141592
L n
が得られることになる。そんな安易でいいのかなあ∼。
3
2
もっと簡単なやり方
針を投げないでもビュッフォンさんと同じようなことができる。今、n を適当な自然数（つまり、1, 2, 3, . . . ）
とし、x 座標も y 座標も −n と n の間にある整数点（x 座標も y 座標も整数の点）を考える。−n と n の間に挟ま
れる矩形領域の一辺の長さは 2n だから、植木算によってこのような点は全部で (2n+1)×(2n+1) = (2n+1)2
個ある。さて、そのような点 (p, q) を１つ取り出したときそれが、原点 (0, 0) を中心とした半径 n の円の内
部に入るかどうかは
p2 + q 2 < n 2
(2.1)
が成り立つかどうかを調べればよい。しかも、このような点の総数と (2n + 1)2 の比は、おおよそ矩形領
n
-n
n
-n
図 4:
域の面積 (2n)2 と、原点 (0, 0) を中心とした半径 n の円の面積 πn2 の比に近いはずである。つまり、上式
を満たす整数点の総数を sn と書くと
sn
πn2
π
≈
=
2
(2n + 1)
(2n)2
4
であるから、π の近似として
π≈
4sn
(2n + 1)2
が得られる。直感的には n を大きくして矩形領域を拡大すればこの近似は良くなると考えられる。これで、
痛い思いをしながら針を投げなくても、(2n + 1)2 個の点一つ一つについて (2.1) を計算すればよいことに
なった。
しかし、何かを投げたいひとにはそのようなやり方もある。今度は、針ではなく小さな小さなビーズを投
げることにしよう。紙の上に正方形とその内接円を描く。そこにビーズを投げて、正方形に入った総数と円
に入った総数の比を計算する。明らかにその比の４倍は π の近似となる。でも、正方形が小さいと、たく
さんのビーズがその外に飛び出してしまってくたびれ儲けとなる。正方形が大きすぎると、その上に均等に
ビーズを投げるのが難しい。実際にビーズを投げて π を求めるのは余り賢くなさそう。
4
3
えーっ？これ、みんな π なの？
以上述べたビュッフォンの方法は π を求める方法としては最悪かも。π を求めるには三角関数の無限級数
展開に基づいたものなど多くの公式が知られており、モンテカルロ法によって π の多くの桁を計算しよう
という酔狂で忍耐強い御仁は今や見あたらない。下にその公式のいくつかと π を記憶するための符丁を挙
げておく。美しさを味わってほしい。インドの天才ラマヌジャンなんかはも∼ぉ最高。ただし、arctan を
用いた公式はこれ以外にもたくさんたくさんある。
4
n2 − 12 + n2 − 22 + · · · + n2 − n2
n2
1 1 1
π = 4 1 − + − + . . . ライプニッツ
3 5 7
24n × (n!)4
ウォリス
π = lim
n→∞ n((2n)!)2
π = lim
n→∞
1
π =3+
ブラウンカー
1
7+
1
15 +
1
1+
292 +
1 + ···
1
1
1
+
+
+ ···
1 × 3 5 × 7 9 × 11
π=8
∞
n=0
√
ライプニッツ
4(−1)n
(−1)n
−
2n+1
(2n + 1) × 5
(2n + 1) × 2392n+1
π=4
π=
1
6×
π =2 1+
1
10
1
π =2+
3
π =3+
1 1
1
+ +
+ · · · オイラー
4 9 16
1 1×2 1×2×3 1×2×3×4
+
+
+
+ · · · オイラー
3 3×5 3×5×7 3×5×7×9
1
1
1
1+
4+
1+
5 + ···
一体誰なの？
10
10
10
2
3
4
2+
2+
2+
2 + ···
オイラー
5
7
9
1+
∞
9801
π= √
8
(4n)!(1103 + 26390n)
(n!)4 3964n
n=0
−1
ラマヌジャン
∞
π=
(−1)n (6n)!(13591409 + 545140134n)
12
(3n!)(n!)3 6403203n+3/2
n=0
π = 16 arctan
∞
π=
マチン
1
n
16
n=0
−1
1
1
− 4 arctan
マチン
5
239
4
2
1
1
−
−
−
8n + 1 8n + 4 8n + 5 8n + 6
チュドゥノフスキー兄弟
ベイリー、ボールウェイン、プラウフ
π ≈ Yes, I have a number.
π ≈ May I have a large container of coffee?
π ≈ Wie? O! Dies π macht ernstlich so vielen M¨
uh!
π ≈ 産医師異国に向かう産後厄なく産婦みやしろに虫散々闇に鳴
5
4
計算の原理
たとえばライプニッツの公式を使って π を 1000 桁まで計算したければ、
4 = +4.000000 . . . . . . . . . . . . . . . . . . . . . . . . . . . 000000
−4/3 = −1.333333 . . . . . . . . . . . . . . . . . . . . . . . . . . . 333333
4/5 = +0.800000 . . . . . . . . . . . . . . . . . . . . . . . . . . . 000000
−4/7 = −0.571428 . . . . . . . . . . . . . . . . . . . . . . . . . . . 571428
とすべての項を（桁上がりを考えて）1000 桁より少し多く書くことから始めなければならない。だから、単
に公式をそのままプログラムすればすむというものではなく、そこがプログラマの腕の見せ所。下に 2400
桁を計算してくれるおそらく最も短いであろう C のプログラムを挙げておく。改行くらい入れてよという
声が聞こえそうだが。
int a=10000,b,c=8400,d,e,f[8401],g;main(){for(;b-c;)
f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),
e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}
5
π 計算の変遷
これまで計算された π を次の表に振り返ってみよう。日本の小学校で教える π が３になるということは2 、
紀元前 1200 年の中国の時代に戻るということ？そんなんでいいのかなあ。
世界記録保持者の金田 (かなだ) 康正東大教授達がまたしても記録を更新した。2002 年 12 月 06 日付の日
経新聞3 は次のように報じている。
円周率１兆 2000 億けた計算、東大が世界新記録「３. １４…」と数が無限に続く円周率を、小
数点以下 1 兆 2411 億けたまで計算し、円周率計算の世界新記録を作ったと、金田康正東大教授
らが 6 日、発表した。従来の記録は、やはり金田教授らが 1999 年 9 月に作った 2061 億けた。今
回は、新しい計算方法を使い、けた数を一気に約６倍に伸ばした。小数点以下 1 兆けた目の数
字は 2、1 兆 2411 億けた目の数字は 5 だった。金田教授は、東大情報基盤センターの並列スー
パーコンピューター「日立ＳＲ８０００／ＭＰＰ」を使い、日立製作所の社員らと今年 9 月か
ら計算に着手。主計算と検証計算などを合わせ、約 600 時間かけて新記録を作った。途中、計
算を止めたこともあり、すべての計算が終わったのは 11 月 24 日だった。
6
忍耐強く実験してみよう（２時間目）
用具：方眼紙、針、電卓
手順 1 約 500 本 (これを N で表す) の針を用意する。針の長さ
を測定する。平行線の間隔 L（ただし
L > ）を何種類か設定する（８の課題１を考慮せよ）。以降の手順は各 L それぞれに対して行なう。
手順 2 方眼紙に間隔 L の平行線を引く。
手順 3 用意した N 本の針を適当な高さから方眼紙の上に落とす。
2 計算が大変なときは
3 を使っても良いというのが指導要領だとのこと
3 http://satellite.nikkei.co.jp/news/main/
6
表 2: π 計算の変遷
バビロニア
BC2000
古代エジプト
BC2000
中国
BC1200
BC250
130
アルキメデス
後漢書
計算値
1
3 + = 3.125
82
16
= 3.16049
9
3
3.14185
√
10 = 3.1622
正確な小数部分の桁数
1
1
0
3
1
1000 年ほどすっ跳ばす
1220
1665
1700
1706
1719
1874
フィボナッチ
ニュートン
関 孝和
マチン
ドゥ・レニュイ
ウィリアム・シャンクス
3.141818
3
16
10
100
112
527
ここから 20 世紀
ライトウィーズナー達
ジェニューイ
シャンクス、レンチ
ギュー、ブイエ
金田、吉野、田村
金田、田村、久保、小林、花村
チュドゥノフスキー兄弟
金田、高橋
金田、高橋
1949
1958
1961
1973
1982
1987
1989
1997
1999
2037
10000
100265
1001250
16777206
134217700
1011196691
51539600000
206158430000
2002
1241100000000
ここから君たちの世紀
金田 ∗
手順 4 平行線と交差している針の本数 n を数える。
手順 5 π の推定値 π
ˆ を求める。
手順 6 手順 3 から手順 5 を 10 回繰り返す。その結果、π の 10 個推定値 π
ˆ1 , π
ˆ2 , · · · , π
ˆ10 が得られる。
手順 7 π
ˆ1 , π
ˆ2 , · · · , π
ˆ10 の平均、標準偏差、95%区間推定を求める。
7
平均・標準偏差・区間推定（３時間目）
m 個のサンプルデータ x1 , x2 , · · · , xm が得られたとする。(本実験では、例えば、π
ˆ1 , π
ˆ2 , · · · , π
ˆ10 がサン
プルデータとなる。)
1. 平均 x
¯
x
¯=
1
(x1 + x2 + · · · + xm )
m
7
2. 偏差平方和
m
m
(xi − x
¯)2 =
i=1
x2i −
(
m
i=1
m
i=1
3. 標本（不偏）分散 s2
xi )2
m
x2i − m¯
x2
=
i=1
m
i=1 (xi
−x
¯)2
m−1
s2 =
4. 標準偏差 s
s=
√
s2
5. 母平均 µ の 95%区間推定
上記の x
¯ と s と m から母平均 µ が 95% の確からしさで存在する区間を作ることができる。最後のメ
モに示したように、t 分布の上側確率 2.5% 点 t(m − 1, 0.025) を用いると
s
s
P x
¯ − t(m − 1, 0.025) √ < µ < x
¯ + t(m − 1, 0.025) √
m
m
= 0.95
(7.1)
√
が分かる。ここで P [ ] は [ ] 内の起こる確率を示す。つまり、未知の母平均 µ が x
¯ −t(m−1, 0.025)s/ m
√
とx
¯ + t(m − 1, 0.025)s/ m との間にある確率は 0.95 であることになる。従って、t(9, 0.025) = 2.262
であるから本実験から得られる円周率 π の 95% 区間推定は m = 10 の場合
s
s
x
¯ − 2.262 × √ , x
¯ + 2.262 × √
10
10
(7.2)
となる。
0.3
0.2
0.025
0.025
0.1
-4
-2
2
4
図 5: 自由度 9 の t-分布と 2.5%点
8
課題
1. 平行線の間隔 L を適当に３つ以上取り上げて、L の違いによる π の推定精度への影響について考察し
なさい。
8
2. 上記の実験で設定した針の長さ l、平行線の間隔 L について、試行回数を N = 100, 500, 1000 と変化
させた場合に、π
ˆ が π の値に最も近くなる交差回数 n∗ を求めよ。また、n = n∗ ± 1 となった場合の
π の推定値をそれぞれの場合について求めよ。これらの結果に基づき、ビュッフォンの針の実験によ
る π の値の推定精度について考察しなさい。
3. その他、気がついた点について述べなさい。
✓
✏
メモ
平均 µ、分散 σ 2 の正規分布 (これを N (µ, σ 2 ) を書く)
正規分布と Student の t 分布
とはその確率密度関数 f (x) が
f (x) = √
1
(x − µ)2
exp −
2σ 2
2πσ
で与えられる分布であるa 。また、自然数 n に対して確率密度関数
1 Γ((n + 1)/2)
1
gn (x) = √
2
Γ(n/2) (1 + x /n)(n+1)/2
nπ
を持つ分布を自由度 n の t 分布b という。ここで、Γ(·) は
∞
Γ(x) =
tx−1 exp(−t)dt
0
で定義されるガンマ関数である。両者の間に以下の関係が知られている。
「X1 , . . . , Xm が互いに独立に正規分布 N (µ, σ 2 ) に従うとき、その平均と分散を
¯=
X
m
i=1
Xi
m
とすると
,
S=
√
¯ 2
− X)
m−1
m
i=1 (Xi
¯ − µ)
m(X
S
は自由度 m − 1 の t 分布に従う。」
a この分布は偶然誤差の分布法則として数学者
Johann Carl Friedrich Gauss (April 30, 1777 - February 23, 1855) によっ
て提案された。ユーロ導入前のドイツの 10 マルク紙幣には、彼の肖像の横に正規分布の関数形とグラフが印刷されていた。そ
れを手にして喜んでいた私に向かって友人は「10 マルクもらってこんなに喜んだ人を見たことがない」とあきれていた。
b この分布は Guinness で働いていた統計技師 William S. Gosset によって 1908 年に発見された。会社との契約で本名で発
表できなかった彼はペンネームとして Student を選んだため、“Student’s t distribution” とも呼ばれる。しかしなぜ “t” が
えらばれたのかはよくわからない。
✒
✑
9
✓
✏
メモ
平均の区間推定
自由度 m − 1 の t 分布 gm−1 (x) に対して
t
+∞
gm−1 (x)dx = 0.975 あるいは同じことだが
−∞
gm−1 (x)dx = 0.025
t
なる値 t (この値を通常 t(m − 1, 0.025) と書き、上側確率 2.5% 点と呼ぶ) を決めると、t 分布が y 軸に
対して左右対称であるa ことから
t(m−1,0.025)
gm−1 (x)dx = 1 − (0.025 + 0.025) = 0.95
−t(m−1,0.025)
となるので、
√
¯ − µ)
m(X
< t(m − 1, 0.025)
S
となる確率が 0.95 であることになる。この条件は簡単な式変形により
−t(m − 1, 0.025) <
¯ − t(m − 1, 0.025) √S < µ < X
¯ + t(m − 1, 0.025) √S
X
m
m
¯ − t(m − 1, 0.025)S/√m と X
¯ + t(m −
と同値である。つまり、95%の確からしさで、未知の平均 µ は X
√
1, 0.025)S/ m の間にあることが言える。t(m − 1, 0.025) の値は統計学の本に載っているので、実際
に上式の 2 つの値を計算することができる。ほんの一部を次のページに載せておいた。このように未
知のパラメータを区間で推定することを区間推定といい、統計学の基本的な考え方の一つである。
a 偶関数であるという
✒
✑
参考文献
[1] N.T. Gridgeman, Geometric Probability and the number π. Scripta Mathematica, 25 (1960), 3, 183–
195.
[2] ジャンポール・ドゥラエ（畑政義訳）「π–魅惑の数」朝倉書店、2001 年
[3] ビュッフォンの針を投げてくれる Java Applet がある：
http://www.mste.uiuc.edu/reese/buffon/buffon.html
10
表 3: t 分布パーセント点
自由度
0.10
0.05
0.025
0.01
0.005
0.001
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
3.078
1.886
1.638
1.533
1.476
1.440
1.415
1.397
1.383
1.372
1.363
1.356
1.350
1.345
1.341
1.337
1.333
1.330
1.328
1.325
6.314
2.920
2.353
2.132
2.015
1.943
1.895
1.860
1.833
1.812
1.796
1.782
1.771
1.761
1.753
1.746
1.740
1.734
1.729
1.725
12.706
4.303
3.182
2.776
2.571
2.447
2.365
2.306
2.262
2.228
2.201
2.179
2.160
2.145
2.131
2.120
2.110
2.101
2.093
2.086
31.821
6.965
4.541
3.747
3.365
3.143
2.998
2.896
2.821
2.764
2.718
2.681
2.650
2.624
2.602
2.583
2.567
2.552
2.539
2.528
63.657
9.925
5.841
4.604
4.032
3.707
3.499
3.355
3.250
3.169
3.106
3.055
3.012
2.977
2.947
2.921
2.898
2.878
2.861
2.845
318.313
22.327
10.215
7.173
5.893
5.208
4.782
4.499
4.296
4.143
4.024
3.929
3.852
3.787
3.733
3.686
3.646
3.610
3.579
3.552
11
© Copyright 2024 ExpyDoc
