IA1 - 布豐投針，圓周率\(\pi\)的估計
 |
| Mathematica®生成400根針的隨機投放示意圖 |
布豐投針問題 (Buffon's needle problem) 是這樣敘述的：
有一房間的地板由數片平行並列的木板所構成，板與板之間溝紋的距離皆相同，如上圖棕色隔線的間距。試問隨機扔一針，該針觸及溝紋的機率為何？
這項問題主要由18世紀的布豐伯爵所提出的，並由義大利數學家Mario Lazzarini在1901年實地進行拋針實驗，共拋出了3,408根針(🥴)並計算出\(\pi\)的近似值約為\(355/113\)，參見維基百科。在我們繼續追本溯源這個歷史問題以前，讓我們先看看當代最常見估計圓周率的方法。
現在我們相信這些處於單位圓內點的數量和全部點的數量的比值應等於那\(1/4\)圓面積和正方形面積之比，即 \[ \frac{n}{N}=\frac{A_圓}{A_方}.\tag{1} \] 而\(A_正 = 1^2\)和\(A_圓=\pi 1^2/4\)，所以我們便得到了 \[ \pi = 4\frac{n}{N},\tag{2} \] 正是\(\pi\)的估計。我們預測若\(N\)愈大，則其估計值愈接近真正的圓周率。另外\(\pi\)可以解析表示成 \[ \pi = \int_{-1}^{1}\frac{dx}{\sqrt{1-x^2}}. \tag{3} \] 上式是由現代分析學之父卡爾·魏爾斯特拉斯 (Karl T. W. Weierstraß, 1815-1897) 在1841年所定義的。
我們接著可以透過rand來實現這個簡單的蒙地卡羅模擬，腳本如下：
所以不同\(N\)值隨機產生的點如圖2，可以發現當灑的點愈多的時候，圓內和圓外之間所形成的邊界便愈平滑，所以內外點數的比值也會接近兩者面積之比(1)，故估計出來的\(\pi\)也會預期會愈精細。
而在圖3我們把不同\(N\)所估計出來的值和解析的值放在一起做比較，當\(N<\mathcal{O}(100)\)的時候MC \(\pi\)的結果震盪極大，而且也與真實的值不穩合度很高。這是由於隨機灑點的結果可能有時有較多的點落於圓內，而估計出來的\(\pi\)便太大，反之則太小，但當\(N>\mathcal{O}(100)\)時，大部分的\(\pi\)估計值皆與真實的值相去不遠了，這也可以從圖2最右邊的子圖中看出來。現在我們已經對隨機灑點估計圓周率的方法有一個概念性的了解了，這個方法其實非常的直觀，所以也是最常被提出來說明的方法，接下來我們要繼續看看歷史上著名布豐投針！
在圖4的示意圖裡，方便起見我們取針的長度和溝紋間距都為\(b\)，我們只需要計算第一個溝和第二個溝也就是在\(0\)到\(b\)之間的結果即可，其它溝之間的情況完全是一致的。我們先看針的落點\(x\)的位置，以藍色的頭為基準點，它會在\(0\leq x\leq b\)之間，所以針落在這個區間任意其中一值的機率是 \[ p_x = \frac{1}{b}. \] 另外針與水平線夾\(\phi\)角，這裡的\(\phi\)是銳角，因為大於\(\pi/2\)的情況與小於\(\pi/2\)的情況是對稱的，不需要重複計算，故\(0\leq \phi \leq \pi/2\)，所以每次針落下處在這個角度的任意機率是 \[ p_\phi = \frac{2}{\pi}. \] 總括來說，當針落下時，它處的位置與夾角的組合會是上述二者的乘積 (落下的距離與夾角是獨立事件) \[ p_{針}(x,\phi)=p_x p_\phi =\frac{2}{\pi b}.\tag{4} \] 另外針在橫軸上的投影是\(b\cos\phi\)，如果要碰到位在0處的溝，我們必須要符合以下條件 \[ x-b\cos\phi<0 \tag{5} \] 才行。亦即如果\(\phi\)是可以任意在\([0,\pi/2]\)之間取值的話，則落點\(x\)必須滿足\(0\leq x\leq b\cos\phi\)才能保證針會碰到位於\(0\)處的溝紋。我們最後寫下針會碰到溝的機率是 \[ p_{碰溝}=\int_0^{\pi/2} \int_0^{b\cos \phi} p_針 dx d\phi=\int_0^{\pi/2} \int_0^{b\cos \phi}\frac{2}{\pi b}dx d\phi=\frac{2}{\pi}.\tag{6} \] 所以若我灑了\(N\)根針，計算後發現有\(n\)根與溝相交，則\(p_{碰溝}=n/N\)，利用(6)式很快得到了 \[ \pi = \frac{2}{p_{碰溝}}=2\frac{N}{n}\tag{7} \] 便將之與\(\pi\)連接起來了，這便是布豐投針估計圓周率的方法，當知道有多少針與溝相碰，則\(\pi\)的值便可估計出來了，這真的是歷史上的巧合！
另外不像本篇專題一開始那張Mathematica生成的示意圖，實際利用投針法估計\(\pi\)時針在縱軸\(y\)上的位置完全無關緊要，這也是合理的，因為針在\(y\)軸上的投影\(b\sin\phi\)永遠不會與溝相交，兩者互為的平行線。還有就是在這個例子裡我們取了針的長度\(\ell\)和溝的間距\(b\)是相等的假設，實際上它們兩者是可以不相等的，這時(5)便要改寫成\(x<\ell \cos\phi\)且須分別討論\(\ell < b\)或\(\ell>b\)的情況，但這只具數學上的趣味，對算法而言沒有任何提昇效率與簡化的幫助，我們在架構算法時，取\(\ell=b=1\)是最簡單的作法。
所以讓我們先來看看如何產生\(\phi\)，或者說是\(\cos\phi\)的值。考慮一單位圓，若任意生成圓內的點座標\((dx,dy)\)我們立即有\[\cos\phi=\frac{dx}{\ell}\]其中\(\ell=\sqrt{dx^2+dy^2}\)。因為要求要銳角故讓\(dx\)和\(dy\)都集中在第一向象限的地方產生，這就是一個製造\(\cos\phi\)的簡單方法，而且完全不需要預先知道\(\pi\)的值。所以投針算法改寫如下：
我們將布豐投針法估計的\(\pi\)值作於圖6，顯然地與灑點法相同，當針的數量愈多，估計出來的\(\pi\)值便愈精細，也愈接近真實的值。
一、灑點估計法
這個方法概念非常的簡單，我們將產生\(N\)個座標\((x,y)\)隨機散佈在\(0\)和\(1\)之間的點，所這些點將被包在以原點為中心，邊長為2的正方形的第一象限中 (first quarter)。接著我們計算出每個點距原點的長度\(d=\sqrt{x^2+y^2}\)，若\(d<1\)則紀錄，反之則丟棄，最後我們紀錄了\(n\)個點，這些點都被包在以原點為中心，半徑為\(1\)的\(1/4\)圓內，見圖1。 |
| 圖1. 在第一象限中藍色正方形的面積為\(1\)而紅色區域的面積為圓的\(1/4\) |
現在我們相信這些處於單位圓內點的數量和全部點的數量的比值應等於那\(1/4\)圓面積和正方形面積之比，即 \[ \frac{n}{N}=\frac{A_圓}{A_方}.\tag{1} \] 而\(A_正 = 1^2\)和\(A_圓=\pi 1^2/4\)，所以我們便得到了 \[ \pi = 4\frac{n}{N},\tag{2} \] 正是\(\pi\)的估計。我們預測若\(N\)愈大，則其估計值愈接近真正的圓周率。另外\(\pi\)可以解析表示成 \[ \pi = \int_{-1}^{1}\frac{dx}{\sqrt{1-x^2}}. \tag{3} \] 上式是由現代分析學之父卡爾·魏爾斯特拉斯 (Karl T. W. Weierstraß, 1815-1897) 在1841年所定義的。
我們接著可以透過rand來實現這個簡單的蒙地卡羅模擬，腳本如下：
import numpy as np
import matplotlib.pyplot as plt
def mcpi(N):
"""
Define Monte-Carlo pi estimation mcpi(N)
N: number of sampling points
"""
# Randomly generate N points with coordinate (x,y) between (0,1)
X,Y = np.random.uniform(size=(2,N))
# Counting how many points lie within the unit circle
NinCircle = np.sum(X**2+Y**2<1)
# Returning the estimation of pi
return 4*NinCircle/N
%config InlineBackend.figure_format = 'retina'
所以不同\(N\)值隨機產生的點如圖2，可以發現當灑的點愈多的時候，圓內和圓外之間所形成的邊界便愈平滑，所以內外點數的比值也會接近兩者面積之比(1)，故估計出來的\(\pi\)也會預期會愈精細。
 |
| 圖2. 不同\(N\)生成的點處在\(1/4\)單位圓內的點 |
而在圖3我們把不同\(N\)所估計出來的值和解析的值放在一起做比較，當\(N<\mathcal{O}(100)\)的時候MC \(\pi\)的結果震盪極大，而且也與真實的值不穩合度很高。這是由於隨機灑點的結果可能有時有較多的點落於圓內，而估計出來的\(\pi\)便太大，反之則太小，但當\(N>\mathcal{O}(100)\)時，大部分的\(\pi\)估計值皆與真實的值相去不遠了，這也可以從圖2最右邊的子圖中看出來。現在我們已經對隨機灑點估計圓周率的方法有一個概念性的了解了，這個方法其實非常的直觀，所以也是最常被提出來說明的方法，接下來我們要繼續看看歷史上著名布豐投針！
 |
| 圖3. \(\pi\)的估計值，藍點為mcpi的結果，紅線為精確值\(3.14159\cdots\) |
二、布豐投針估計法
布豐投針估計\(\pi\)的方法是一個歷史上的巧合，原來布豐伯爵僅僅是想要知道針與溝紋相交的機率罷了，至於為何會與\(\pi\)的估計扯上關係？讓我們先看看圖4的示意圖。 |
| 圖4. 布豐投針示意圖 |
在圖4的示意圖裡，方便起見我們取針的長度和溝紋間距都為\(b\)，我們只需要計算第一個溝和第二個溝也就是在\(0\)到\(b\)之間的結果即可，其它溝之間的情況完全是一致的。我們先看針的落點\(x\)的位置，以藍色的頭為基準點，它會在\(0\leq x\leq b\)之間，所以針落在這個區間任意其中一值的機率是 \[ p_x = \frac{1}{b}. \] 另外針與水平線夾\(\phi\)角，這裡的\(\phi\)是銳角，因為大於\(\pi/2\)的情況與小於\(\pi/2\)的情況是對稱的，不需要重複計算，故\(0\leq \phi \leq \pi/2\)，所以每次針落下處在這個角度的任意機率是 \[ p_\phi = \frac{2}{\pi}. \] 總括來說，當針落下時，它處的位置與夾角的組合會是上述二者的乘積 (落下的距離與夾角是獨立事件) \[ p_{針}(x,\phi)=p_x p_\phi =\frac{2}{\pi b}.\tag{4} \] 另外針在橫軸上的投影是\(b\cos\phi\)，如果要碰到位在0處的溝，我們必須要符合以下條件 \[ x-b\cos\phi<0 \tag{5} \] 才行。亦即如果\(\phi\)是可以任意在\([0,\pi/2]\)之間取值的話，則落點\(x\)必須滿足\(0\leq x\leq b\cos\phi\)才能保證針會碰到位於\(0\)處的溝紋。我們最後寫下針會碰到溝的機率是 \[ p_{碰溝}=\int_0^{\pi/2} \int_0^{b\cos \phi} p_針 dx d\phi=\int_0^{\pi/2} \int_0^{b\cos \phi}\frac{2}{\pi b}dx d\phi=\frac{2}{\pi}.\tag{6} \] 所以若我灑了\(N\)根針，計算後發現有\(n\)根與溝相交，則\(p_{碰溝}=n/N\)，利用(6)式很快得到了 \[ \pi = \frac{2}{p_{碰溝}}=2\frac{N}{n}\tag{7} \] 便將之與\(\pi\)連接起來了，這便是布豐投針估計圓周率的方法，當知道有多少針與溝相碰，則\(\pi\)的值便可估計出來了，這真的是歷史上的巧合！
另外不像本篇專題一開始那張Mathematica生成的示意圖，實際利用投針法估計\(\pi\)時針在縱軸\(y\)上的位置完全無關緊要，這也是合理的，因為針在\(y\)軸上的投影\(b\sin\phi\)永遠不會與溝相交，兩者互為的平行線。還有就是在這個例子裡我們取了針的長度\(\ell\)和溝的間距\(b\)是相等的假設，實際上它們兩者是可以不相等的，這時(5)便要改寫成\(x<\ell \cos\phi\)且須分別討論\(\ell < b\)或\(\ell>b\)的情況，但這只具數學上的趣味，對算法而言沒有任何提昇效率與簡化的幫助，我們在架構算法時，取\(\ell=b=1\)是最簡單的作法。
2.1 投針算法
所以讀者可能有想法要如何利用上面的介紹來架構投針算法，如果你的想法是像下面這樣：- 產生針頭的落點 (圖4中的藍點)：\(x= {\rm \tt rand} (0,1)\)
- 隨機產生針與橫軸的夾角：\(\phi={\rm \tt rand} (0,\pi/2)\)
- 計算(5)式：\( x^\prime = x-\cos\phi\)
- \({\tt if}~(x^\prime<0):~{\tt return}~1\)
\({\tt else}:~{\tt return}~0\)
 |
| 圖5. 單位圓內的點 |
所以讓我們先來看看如何產生\(\phi\)，或者說是\(\cos\phi\)的值。考慮一單位圓，若任意生成圓內的點座標\((dx,dy)\)我們立即有\[\cos\phi=\frac{dx}{\ell}\]其中\(\ell=\sqrt{dx^2+dy^2}\)。因為要求要銳角故讓\(dx\)和\(dy\)都集中在第一向象限的地方產生，這就是一個製造\(\cos\phi\)的簡單方法，而且完全不需要預先知道\(\pi\)的值。所以投針算法改寫如下：
- 產生針頭的落點：\(x= {\rm \tt rand} (0,1)\)
- 產生單位圓第一象限內的點：
\(dx= {\rm \tt rand} (0,1)\)
\(dy= {\rm \tt rand} (0,1)\)
\(\ell=\sqrt{dx^2+dy^2}\) - \({\tt if}~(\ell>1\)): 回到第2步 (產生的點不在單位圓內)
- 計算(5)式：\( x^\prime= x-(dx/\ell)\)
- \({\tt if}~(x^\prime<0):~{\tt return}~1\)
\({\tt else}:~{\tt return}~0\)
def buffonpi(x):
"""
Define pi estimation by Buffon's needle buffonpi(x).
This function starts from Step 2 of the algorithm.
x: Initial x-position of the needle between (0,1)
"""
# Step 2
dx, dy = np.random.uniform(0, 1, 2)
l = np.sqrt(dx**2 + dy**2)
# Step 3
if l > 1:
"""
python does not support goto natively!
In order to go back to Step 2, we have to call itself again and
feed it with the same input x.
"""
return buffonpi(x)
# Step 4
else:
xprime = x - (dx/l)
# Step 5
if xprime < 0:
return 1
else:
return 0
我們將布豐投針法估計的\(\pi\)值作於圖6，顯然地與灑點法相同，當針的數量愈多，估計出來的\(\pi\)值便愈精細，也愈接近真實的值。
 |
| 圖6. 投針法估計\(\pi\)值 |
1 網路上似乎有模擬goto的模組，不過python在開發時就有意拿掉goto，目的是為了避免程式運作的流程過於混亂
三、小結
即使布豐伯爵一開始在想投針問題的時候，他本人也不會預期到這個問題的解答與圓周率的估計有著緊密的連結。因為估計針與溝紋相交的機率仰賴知道\(\pi\)的值，而在布豐的年代，\(\pi\)的值是毋庸置疑的，只是後人發現，如果我可以實作這個問題，從投針的結果我可以逆回去推出\(\pi\)的值。從觀測出發，得出待定參數\(\pi\)的估計，這也與另一個系列P&S的概念息息相關。
而在第二節裡討論的兩種算法當中，我們也提到了，真正的蒙地卡羅方法估計未知參數時，是完全不需要仰賴這個參數的任何資訊，否則就不算是真正的在估計未知的參數，因為我們作弊了！所謂的蒙地卡羅如同賭博一樣，每次產生結果都是隨機的，但僅有當符合條件的結果才會被記錄下來，正如同僅有賭贏的人才能進到下一個關卡，否則便是賠光籌碼掃地出門！在灑點估計裡，紀錄那些落在圓內的點和所有的點相比較，我們便能體現出自然的規律 (即圓周與半徑之比\(\pi\))，即便我們對這個規律在開始時一無所知，故蒙地卡羅法是非常符合自然界運作方式的一種逼近。
Source on github: I&A1_buffon.ipynb
而在第二節裡討論的兩種算法當中，我們也提到了，真正的蒙地卡羅方法估計未知參數時，是完全不需要仰賴這個參數的任何資訊，否則就不算是真正的在估計未知的參數，因為我們作弊了！所謂的蒙地卡羅如同賭博一樣，每次產生結果都是隨機的，但僅有當符合條件的結果才會被記錄下來，正如同僅有賭贏的人才能進到下一個關卡，否則便是賠光籌碼掃地出門！在灑點估計裡，紀錄那些落在圓內的點和所有的點相比較，我們便能體現出自然的規律 (即圓周與半徑之比\(\pi\))，即便我們對這個規律在開始時一無所知，故蒙地卡羅法是非常符合自然界運作方式的一種逼近。
Source on github: I&A1_buffon.ipynb
留言
張貼留言