Jserv's blog

Buffon 投針試驗：當圓周率計算遇上機率論

1999 年，一向以優雅高檔著稱的 [Givenchy]，推出由知名設計師 Serge Mansau 所設計的「π香水」(pi)，據稱，是一款令人精神充沛的木質味香水，瓶身造型展現π的精神：無止盡的學習與探索，而不規則的密度感、基於琥珀色調的雕刻圖案，使其得以在眾多男用香水中脫穎而出。是的，計算圓周率一直是令人著迷的議題，從古埃及至今，無數專家學者乃至業餘數學家前仆後繼地投入，1777 年，Georges-Louis Leclerc，Comte de Buffon (1707—1788) 提出嶄新的途徑，將圓周率這等幾何問題出發的計算，巧妙地以機率統計原理來表示，自此，開創使用隨機數值處理典型數學表示的先河，我們就來看看傳奇性的 Buffon 投針試驗。

法國數學家、科學家、《自然史》作者，也是風格家的 Comte de Buffon 在 1777 年某日，邀請賓客齊聚大廳，共襄盛舉一次試驗活動。古稀之年的 Buffon 鋪好一張白紙，其上預先畫好了一條條等距的平行線，接著取出一大把質量均等、長度為平行線間距一半的小針，待賓客就座後，Buffon 發言道：

「煩請各位將這些小針一根一根扔往白紙上，並且告知扔下的針是否與紙上平行線相交」

(示意圖)
客隨主意，雖摸不著頭緒，但也一個個加入了試驗的行列。一把小針扔完了，把它撿起來又扔，而 Buffon 則在一旁不停地記數著，忙碌了將近一個鐘頭。最後，Buffon 高聲宣佈：

「各位賓客，依據我的紀錄，剛才的投針結果，共投針 2212 次，其中與平行線相交有 704 次。而總數 2212 與相交數 704 的比值為 3.142。」

說到這裡，Buffon 故作停頓，神秘張望賓客，接著說：

「這就是圓周率π的近似值！」

Buffon 利用平凡不過的除法，計算出圓周率的近似值，並宣稱投針的數目越多，圓周率的近似值將會越精準，這就是數學史上著名的 Buffon 投針問題，記載於其著作《機率算術試驗》(1777 年)，此外，Buffon 也試著將機率應用於審判場合，比方說，若能對每個審判員規定某個足以理解真相或說出真相機會的數值，即可算出法庭作出正確判決的機會，換言之，就是「審判的概率」(Probabilite des jugements)。

圓周率π在這種看似雜亂的場合出現，實在出乎意料。一個直觀的理解途徑可透過物理上的對應，取一根鐵絲，將其彎成一個圓圈，適度剪裁使其直徑恰等於平行線間距離 d。於是乎，對於這個圓圈來說，無論如何扔下，都將和平行線有兩個交點。也就是說，若圓圈扔下的次數為 n，那麼，相交的交點總數必為 2n。接著，我們展開物理的形變，將圓圈拉開、拉直，這樣就成為長度為 πd 的鐵絲，再將這條鐵絲扔下，與平行線相交的情況就複雜許多，由於 1 < π < 4，我們可歸納交點數量為 4 個 (跨越三個間隔)、3 個 (跨越兩個間隔)、2 個 (跨越一個間隔)，1 個 (不跨越間隔，僅一端碰觸平行線)，以及 0 個 (幾乎與平行線平行)。注意，在形變的過程中，鐵線在圓圈和直線這兩個型態的長度同為 πd，根據機會均等的原理，在同等大量的投擲下，兩者與平行線組交點的總數應有一致的期望值，換言之，當長為 πd 的鐵絲扔下 n 次，其與平行線相交的交點總數應接近 2n。
(示意圖：考慮AB與XY這兩個間距)
當正多邊形之邊數趨於無限多時，極限就是正圓，在概念上來說，也就是上述產生形變的小針，而小針觸碰線的次數，依據機率統計原理，應比照特定之比例。接著，我們將鐵線切為長度為 (1/2)d 的小線，在物理上大量的扔放小線的行為來說，小線碰線的期望值 m 應保持等比例，也就是說：

πd : (1/2)d ≈ 2n : m

這也是 Buffon 投針試驗中所作的參數配置，約分後可得漂亮的式子：

π ≈ n / m

在古典數學中，求圓周率之值是幾何問題，而 Buffon 卻以此拍案叫絕的方式，以機率方法打通兩個看似風馬牛不相及的領域，成為幾何概率的典型例子。此外，Buffon 還翻譯了牛頓的著作，並探討牛頓和萊布尼茨發現微積分的歷史，集多年研究成果編成巨著《自然史》(四十四卷，1749—1808) ，原定出版五十卷，Buffon 生前僅出版三十六卷，後八卷由他的學生完成，其中包含物種進化的思想，推動了古生物學的發展。涉獵廣泛的 Buffon 也是第一個對地質史劃分時期的科學家，首次提出太陽與慧星碰撞產生行星的理論。

近代科學的發展下，原本壁壘分明的個別人文、科學、哲學思想領域走向空前的大融通，匯流而成當代種種巨大變革，一如 Buffon 首次打破機率與幾何學的藩籬。數學領域的變遷也受到這等啟蒙，1904 年，R·Chartres 甚至提出另一種表示法：若寫下任意兩個整數，測這兩者互質的機率為 6 / π^2。

我們以更嚴謹的態度分析這個實驗結果，分析投下的小針與平行線相交的關係圖，如下：

其中，針與平行線的銳角交角為θ，針的中點到最近一條直線的距離為 Y。我們可得之，若 Y ≦ d/2 * sinθ，則小針壓線，反之，若 Y ＞ d/2 * sinθ，則小針不壓線。進一步來看，在上圖下方中，點 (θ, Y) 若落於矩形中的曲線下部分，針就壓線，反之就是不壓線，因此，針壓線的機率就是曲線下部分在長方形中所占的比例。經過分析後，我們就可著手以機率統計理論來探討。

先來看看 pdf (probability density function) 的表示。Y 在 0 到 d/2 之間的 pdf 為：

而 θ 在 0 與 π/2 之間的 pdf 為：

既然 Y 與 θ 為獨立之隨機數值，換言之，獨立事件之機率為其乘積，表示為：

而由前述分析可得，針與平行線相交的充要條件為：

我們知道，曲線下的面積就是針與平行線相交的機率，透過積分可知：

所以當 n 個細線被扔下，而其中 h 個與平行線相交，其機率為：

簡化後即可求得圓周率：

經過幾百年的演繹與探討，Buffon 投針試驗逐漸演化為一種數值方法的前身：「蒙地卡羅方法」(Monte Carlo method)，也就是透過利用亂數取樣 (random sampling) 模擬來解決數學問題。第二次世界大戰期間，Monte Carlo 方法被系統性地應用於科學研究中，誕生了 MANIAC (Mathematical Analyzer, Numerical Integrator and Computer)，而 Stanislaw Ulam、John von Neumann、Nicholas Metropolis、Enrico Fermi 等人發展法一種基於樣本統計的方法，來解決關於在原子彈設計中，中子隨機擴散問題和 Schrodinger 等式的特徵值估計問題。該方法的原理最初是 Stanislaw Ulam 闡述的，後來由 John von Neumann 深入研究，於 1949 年發表一篇名為 "The Monte Carlo method" 的論文而聞名，當然，到了進入電腦時代，這個方法才得以由原本手動產生亂數來解決問題，變成實際性的數值方法。

Monte Carlo 方法是由 Nicholas Metropolis 所命名，取自其亂數機率有如賭博一般，而恰似北非最西側的摩洛哥首都 Monte Carlo，也就是知名賭城，種種奇豔動人的賭場生活寫照。所有具有隨機效應的過程，均可能以 Monte Carlo 方法大量模擬單一事件，並藉由統計上平均值，獲得某設定條件下實際最可能測量值，更廣泛來說，自然界裏的布朗運動、電波的噪音、基因的突變、交通即時路況等等，無處不含有隨機的變化，均有可適用的場合。

最後，且讓我們實地模擬 Monte Carlo 方法始祖的 Buffon 投針，以下是小弟的 C 語言小實驗：[buffon.c]

/**
 * buffon.c - Simulating Buffon's Needle Problem
 *
 * Authored by Jim Huang <jserv.tw@gmail.com>
 * Public Domain.  Runtime confirming to POSIX.1-2001.
 */
#include <stdio.h>
#include <stdlib.h>

#include <math.h>
#include <time.h>
 
int main( int argc, char **argv )
{
	int i;
	int ntrial = 99991; /* the largest prime number < 100000 */

	int crossing = 0;
	srand48( time( NULL ) );
 
	/*
	 * ----------------------------------------
	 * ||
	 * || 1                         
	 * ||                / ) Theta   ||x
	 * -----------------/----------------------
	 * x <= 1/2 * sin( Theta ), when needles cross a crack.
	 */
	for (i = 0; i < ntrial; i++) {
		double theta = M_PI * drand48(); /* 0 <= Theta <= PI */

		double x = drand48(); /* 0 <= x < 1 */
		if (x <= 0.5 * sin(theta))
			crossing++;
	}
	printf("Needles dropped on floor = %8d \n", ntrial);
	printf("Those that cross a crack = %8d \n", crossing);
	printf("\tPI Estimate: ");
	printf("%0.5f \n", (double)ntrial / (double)crossing);
	return 0;
}

我得到最接近的輸出為：

Needles dropped on floor =    99991 
Those that cross a crack =    31830 
        PI Estimate: 3.14141

其中 99991 是五位數中最大的質數，在隨機試驗中沒有特別的意義。行文至此，彷彿聞到「π香水」給予人充沛與活力的美妙體驗，謎樣的數學，謎樣的π。

參考資料：

• Wikipedia [Buffon's needle]
• George Reese - [Buffon's Needle : An Analysis and Simulation]

由 jserv 發表於 April 4, 2008 04:32 PM