# Buffon's Needle 和蒙地卡羅模擬
###### tags: `others`
Buffon's Needle是一個源自18世紀的著名機率問題。而蒙地卡羅模擬是一種統計方法，用於通過隨機抽樣和重複實驗來估計各種領域中的數量或現象。它以蒙地卡羅賭場為名，因為這種模擬方法的核心思想是使用隨機數來模擬現實世界的隨機性。其基本思路是通過生成大量的隨機樣本，根據這些樣本進行數學模型的計算和模擬實驗。模擬的結果可以用於估計目標事件的機率、計算數值積分、進行風險分析等。一般步驟如下：
1. 定義模型和需要估計的目標。
2. 隨機生成符合模型分佈的樣本數據。
3. 對每個樣本進行計算或模擬實驗。
4. 根據樣本結果進行統計分析，得到目標的估計值。
5. 通過增加樣本數量來提高估計的準確性。
6. 進行敏感性分析和結果解釋。
## 問題描述
想像一個地板上有無限多條平行且等距離的線，彼此之間相距 $\ell$。
現在我們將一根長度為 $L$ （小於 $\ell$）的針隨機投擲到地板上，它會以均勻隨機的位置和均勻隨機的方向落下。Buffon's Needle 問題是要確定這根針與地板上的線相交的機率。這個機率取決於針的長度 $L$ 與地板上線的距離 $\ell$ 的比值。
我們現在就來使用Python程式語言進行蒙地卡羅模擬，估計針穿過線的機率。
## 對地板建模
地板的部分很簡單，只有一個參數，也就是 $\ell$。我們將地板視為一個二維的平面，其中這些直線都和 $y$ 軸平行，並和 $x$ 軸垂直。
我們可以寫一個class來建立地板物件：
```python=
class Floor:
def __init__(self, l):
self.l = l
```
## 對針建模
針在地板上的位置和方向可以用其中心的 $(x, y)$ 坐標和與垂直軸所成的角度 $\theta$ 來描述。
需要注意的是，針的 $y$ 坐標並不會影響它是否與地板上的線相交。因此，在計算中我們不需要考慮 $y$ 坐標，只需在模擬過程中追蹤針的 $x$ 坐標和角度 $\theta$。
在這邊我們為針也建立一個類別，包含的參數有其座標、長度和角度：
```python=
class Needle:
def __init__(self, L, theta, x, y=0):
self.L = L
self.theta = theta
self.x = x
self.y = y
```
## 模擬的狀態空間
為了進行模擬，創建一個無限大的地板是不實際的。相反，我們只在地板的兩條相鄰線之間投擲針，即 $0\leq x\leq \ell$。由於地板的對稱性，針穿過線的機率估計與在無限地板上的情況相同。
對於針與 $y$ 軸所成角度的範圍可以有 $0\leq\theta\leq 2\pi$，理論上是合理的。然而，我們可以利用一些旋轉對稱性。具有角度 $\theta$ 的針看起來等價於具有相同坐標和角度 $\theta + \pi$ 的針（假設針沒有頭尾之分）。因此，在我們的計算中，只允許 $0\leq\theta\leq\pi$。
針的狀態空間 $\mathcal{R}$ 定義為：$\mathcal{R} = \{(x, \theta)|0\leq x\leq \ell, 0\leq\theta\leq\pi\}$。在我們的模擬中，隨機投擲一根針到地板上等價於從狀態空間 $\mathcal{R}$ 中均勻抽取一個樣本。
## 投針的模擬
為了進行模擬，我們可以建立一個名為`toss_needle`的函數，該函數回傳一個在狀態空間中均勻分布的`Needle`物件。該函數的參數包括要投擲的針的長度`L`，以及針被投擲到的`Floor`物件：
```python=
def toss_needle(L: float, floor: Floor) -> Needle:
x = random.uniform(0, floor.l)
theta = random.uniform(0, math.pi)
return Needle(L, theta, x)
```
在這個函數中，我們使用了`random.uniform`函數來生成在指定範圍內均勻分布的隨機數。`x`變數代表針在地板上的 $x$ 坐標，取值範圍是從`0`到地板的寬度`floor.l`。`theta`變數代表針的角度，取值範圍是從`0`到`pi`。最後，函數回傳一個 Needle 物件。
## 確認針和地板上的線是否相交
如果針的一部分與 $x=0$ 或 $x=\ell$ 重疊，代表針與地板上的線相交了。這要怎麼判斷呢？
我們可以通過檢查針的端點的 $x$ 坐標來檢查重疊。如果針的左端點小於 $0$ 或右端點大於 $\ell$，則針與線相交。
針的兩個端點的 $x$ 坐標可以如下計算。首先，取針的中心點的 $x$ 坐標。然後，將針在 $x$ 軸上的投影的長度的一半，也就是 $(L/2)\sin\theta$ 加上和減去，得到一個範圍，如果這個範圍包含了 $0$ 或 $L$，則發生了相交。
因此，檢查針是否與地板上的線相交的程式碼如下：
```python=
def check_intersection(needle, floor):
x_left = needle.x - (needle.L/2) * math.sin(needle.theta)
x_right = needle.x + (needle.L/2) * math.sin(needle.theta)
if x_left <= 0 or x_right >= floor.l:
return True # 針與線相交
else:
return False # 針未與線相交
```
## 蒙地卡羅模擬
現在我們將使用蒙地卡羅模擬來估計投擲的針穿過地板上線的機率。我們將使用`toss_needle`函數在地板上投擲大量的針，並使用`cross_line`函數來跟踪有多少針與地板上的線相交。我們對針穿過地板上線的機率的估計，即投擲的針中穿過線的比例，可以經由解析解算出應該要是 $\frac{2L}{\pi \ell}$。
```python=
def estimate_probability(L, floor, num_needles):
count_crossed = 0 # 記錄穿過線的針的數量
for _ in range(num_needles):
needle = toss_needle(L, floor) # 投擲針到地板上
if check_intersection(needle, floor): # 檢查針是否與地板上的線相交
count_crossed += 1 # 如果相交，增加計數
probability = count_crossed / num_needles # 穿過線的針的機率
expected_probability = 2 * L / (math.pi * L) # 預期的穿過線的機率
return probability, expected_probability
```
Sign in with Wallet
Connect another wallet