如何用一個公平的硬幣模擬 1/Pi 的概率？

11-25

示例：
1、模擬1/2的概率：拋擲1次硬幣即可
2、模擬1/4的概率：拋擲2次硬幣即可，取2個正面的結果
3、模擬1/5的概率：拋擲3次硬幣即可，取3個正面的結果，若出現2個正面則再擲3次硬幣重複之前判定過程
4、依此類推我們可以在有限期望拋擲次數的條件下模擬所有有理數的概率

現在我們需要要用這枚公平的硬幣模擬1/pi的概率，具體要求如下：
1、我們可以保證拋擲正反概率相同（公平的硬幣），但不能保證拋擲落點均勻分布
2、要求概率模擬十分精確，並且拋擲硬幣次數的期望為有限次。（由於pi是無理數所以不能通過逼近的方法模擬十分精確概率了）

(此答案於2014年10月8日更新，添加了一種做法)

statement of mission:
1. 只利用一個能生成Bernoulli(0.5)的隨機數發生器「拋硬幣」。
2. 設計一個with probability 1會在有限次拋硬幣後終止的遊戲，對於任意給定的 $epsilon$ &>0，都能調整遊戲規則，令該遊戲中玩家勝利的概率與 $frac{1}{pi}$ 的誤差小於 $epsilon$ 。

part 1: 設計一個almost surely會有限步結束的子遊戲，令玩家在該遊戲中獲勝的概率為 $frac{m-1}{m}$ ，其中 $m$ 是任意事先給定的正整數， $mgeq 2$ 。

solution to part 1: 首先拋硬幣 $m-1$ 次，如果全為正面那麼玩家失敗，如果恰好有一次正面那麼玩家獲勝，其他情形遊戲重置，再拋硬幣 $m-1$ 次。
玩家失敗的概率比較好算，每一次遊戲玩家失敗的概率為 $frac{1}{2^{m-1}}$ ，需要進行下一次遊戲的概率為 $1-frac{1+(m-1)}{2^{m-1}}=1-frac{m}{2^{m-1}}$ ，所以玩家失敗的概率是 $frac{1}{2^{m-1}}Big/ left( 1-left[1-frac{m}{2^{m-1}} ight] ight) = frac{1}{m}$

part 2: 做有限個子遊戲，依次設置它們的參數 $m$ 為 $2, 4 imes 1^2, 4 imes 2^2, 4 imes 3^2, ldots$ ，最後根據 $frac{1}{pi} = frac{1}{2} imes prod_{n=1}^{infty} frac{4n^2-1}{4n^2}$ (參考Infinite product)，這樣的遊戲序列隨著子遊戲數目的增加，其玩家獲勝的精確概率越來越逼近 $frac{1}{pi}$ ，並且逼近的誤差可以容易地計算出來。可以選取足夠大的遊戲數目，使得逼近誤差小於任意事先給定的值。

================================================================
(以下為2014年10月8日更新：)

這個問題還可以變難一點：如果硬幣不是公平的，即正面朝上的概率為 $p$ 並且 $p$ 未知，怎樣設計符合要求的遊戲？

根據分割線以上的解答，問題仍然可以約化為要求設計一個子遊戲，使得玩家獲勝的概率為 $frac{s}{t}$ ，其中 $0<s<t$ ， $s$ 和 $t$ 是兩個正整數。

遊戲設計：可以認為手中有 $t$ 個這樣的硬幣(因為總可以拋 $t$ 次，分別認為是第 $1,ldots,t$ 個硬幣拋出的結果)。遊戲規則如下：

1. 每一輪遊戲都同時拋出這 $t$ 個硬幣。
2. 如果恰有一枚硬幣正面朝上，則終止並判定勝負；否則進行下一輪拋擲。
3. 勝負判定方法：如果正面朝上的硬幣在第 $1, ldots, s$ 號硬幣之中，則玩家獲勝，否則失敗。

容易知道在給定條件「所有硬幣中恰有一枚正面朝上」時，每一枚特定的硬幣正面朝上的條件概率是相等的，所以都為 $1/t$ ，所以玩家獲勝的概率是 $s/t$ 。

(1) 如果題主要的只是拋擲硬幣次數的期望為有限次的話，那麼無論目標概率是有理數還是無理數，都可以用對半砍區間的方法模擬。

具體方法如下：
初始區間為[0,1)，目標概率為1/pi（約等於0.3183）。
擲一次硬幣。若為反，則砍去區間的左半邊，剩下[0.5, 1)。目標概率完全在此區間左邊，演算法結束，輸出「反」。若為正，則砍去區間的右半邊，剩下[0, 0.5)。目標概率仍在區間內，演算法繼續。
再擲一次硬幣。若為反，砍去區間的左半邊，剩下[0.25, 0.5)，繼續。若為正，砍去右半邊，剩下[0, 0.25)，結束，輸出「正」。
……

可以證明，無論目標概率是多少，拋擲硬幣次數的期望等於2。
但在極端情況下，可能要拋擲很多次硬幣，拋擲硬幣的次數沒有上限。
（當然，當目標概率為無理數時，你需要有一個演算法把它計算到所需的精度）

(2) 如果題主要的是拋擲硬幣的次數有上限的話，那麼無解。
因為有限次的拋擲只能獲得有限個等概率的結果，無論取哪些作為「正」，哪些作為「反」，都只能得到有理數的概率。
（事實上，這樣得到的概率只能是分母為2的整數次冪的有理數。回想一下，模擬1/3概率時有一種情況是「重來」，也就是說，拋擲硬幣的次數沒有上限！）

(3) 當然，對半砍區間法還有一個困難就是怎樣把1/pi算到足夠的精度。
所以，題主想要的還可能是這樣一種演算法：不需要顯式算出1/pi的精確值，而利用它與圓的關係，隱式地獲得1/pi的概率。
這樣的演算法也是存在的！
我們可以採用對半砍矩形的方法。

如圖，藍色正方形邊長為2，綠色的圓是它的內切圓，紅色正方形邊長為1。
綠圓的面積為pi，紅色正方形面積為1。
演算法如下：
while (true) {
取藍色正方形為當前矩形
while (綠色圓或紅色正方形邊界上存在某點在當前矩形內部) {
拋一次硬幣，若為反，砍去當前矩形的左一半；若為正，砍去當前矩形的右一半
再拋一次硬幣，若為反，砍去當前矩形的下一半；若為正，砍去當前矩形的上一半
}
if (當前矩形位於紅色正方形內部) return 「正」
if (當前矩形位於紅色正方形外、綠色圓內) return 「反」
}

這個演算法輸出「正」的概率是紅色正方形與綠色圓的面積之比，即1/pi。
它不需要計算1/pi的值，只需要判斷「當前矩形」的四個角是在圓外還是在圓內。

與對半砍區間法類似，這個演算法拋硬幣的期望次數是有限的。
但拋硬幣的次數沒有上限，原因是兩個while循環在特殊情況下都可能執行很多次。

另外，(3)中的演算法是針對1/pi專門設計的。如果把概率換成別的無理數，比如1/e，就要重新設計演算法。可見其通用性不強。

// 下面這個函數以概率 p 返回true，概率 1 - p 返回false。
// 調用一個以0.5概率返回true的函數 toss。來源：網路。
bool prob(double p) {
if (p &< 0.5) return !prob(1 - p);
return toss() ? true : prob(2 * p - 1);
}

給個非常簡單的辦法：
把1/pi寫成2進小數。
然後不停地投硬幣，產生一個浮點隨機數。第一次得到小數點後第一位的值，第二次得到第二位的值...直到這個不斷變得精確的小數與1/pi有確定的大小關係。
這樣，得到一個小於1/pi的數的概率是1/pi。
你要是運氣好到投這個無窮序列正好是1/pi，踩中了概率為0的事件，我也是服氣。

這麼容易的問題不需要看正反面-_-

不可能，要求2滿足不了。只是普通拋硬幣的話，你能得到的概率被限定在 rational number 裡面，但是 1/pi 本身是一個 real number。Rational number is countable, real number is uncountable。

用對角線方法 [1] 可以證明不可能得到一個一一對應。

Trick: 把硬幣立起來當指針旋轉能得到你要的概率。

[1] http://en.wikipedia.org/wiki/Cantor%27s_diagonal_argument

通過重複隨機模擬試驗來粗暴的獲得一個數值解是典型的的monte carlo模擬。

不是很理解題主的意思，但是我們這裡假設擲硬幣次數足夠大。

首先呢，我們要產生隨機數。因為硬幣只能給出0或者1兩個數值，所以這裡採用二進位。也就是說某個倒霉孩子連續拋硬幣n次，就可以得到一個n位的隨機二進位數。假設這個小孩連續拋了2500次，我們就得到了50個50位的二進位數。我們把這50個二進位數轉換成十進位小數(原值除以2^50)，在畫出直方圖。

呃....看上去不是很隨機哈...我們來拋50萬次看看.....

嗯，這次看上去在每個區間的概率都差不多了，也就是有10000個服從(0,1) uniform distribution的隨機數了。

現在我們假設有個1x1正方形的大靶子，裡面有一個內切的圓形。我們把之前的一萬個隨機數兩兩分組當做坐標(也就是5000組坐標)，如果坐標在圓形內算得1分，在外面算0分。這樣我們就可以通過（得分/總數）估算圓形的面積了。

先來十組坐標看看....

好吧，不是很直觀。現在估算pi是3.602。我們把5000組坐標都扔上去看看....

這下很直觀了，我們根據圓形面積公式得到了Pi的估計值3.1416。不錯哈~

另外附送1/e的估算方法。

首先已知 e=lim(1+1/n)^n 那麼1/e=lim(1-1/n)^n。

假設我們把上面的靶子均勻分成了n個小格子。那麼每個坐標落在一個格子的概率是1/n。那麼某一個格子n組坐標都沒有命中的概率就是(1-1/n)^n。

好了，現在我們設計實驗。之前有個倒霉孩子已經通過500000次擲硬幣給我們製造了5000組隨機坐標。我們只要把靶子割成5000份，然後把這些坐標標註出來。之後在其中數出沒有被擊中過的格子的數量(m)。那麼m/5000就是我們估算的1/e的值。

這些就是monte carlo方法最基本的應用。簡而言之就是在我們對系統啥的不甚了解的情況下，用最簡單最粗暴的隨機實驗來獲得數值解的方法。

為何你們都那麼天真只想著正反面呢，硬幣是圓的，咱們可以標個方向，然後在硬幣上確認一個2弧度的角，這樣不就是1/pi了嗎。←_←

今天看見一個布豐投針問題

設我們有一個以平行且等距木紋鋪成的地板，現在隨意拋一支長度比木紋之間距離小的針，求針和其中一條木紋相交的概率。這就是布豐投針問題（又譯「蒲豐投針問題」）。

布豐投針問題

在這裡面，讓平行線間距是硬幣直徑的兩倍，再在硬幣上沿直徑畫一根線當成針，按照這個實驗做，畫的線與平行線相交的概率就是1/pi

其實很簡單，只要跳出「硬幣一定是用來拋出看正反面」的思維定式就可以了。

因為pi是個超越數，所以你的方法中一定要用到pi，只看正反面是不可能完成的。
所幸硬幣本身就有pi——它是圓的。想到這裡方法就太多了，我隨便說兩個吧。

方法一：
用細繩比出直徑，在硬幣的側面畫出直徑那麼長的線，那麼線長為1，直徑為pi。
隨便拋硬幣，硬幣落下後正上方12點位置落在側面所畫線內的概率就是1/pi。

方法二：
繞著硬幣在紙上畫個圓，再用細繩比出直徑，對摺成半徑，用這段半徑在圓里畫個正方形。
那麼圓的面積為pi，正方形的面積為1。
隨便在硬幣的側面標個點，扔硬幣，該點落在正方形內的概率是1/pi，當然落在圓形外的次數不算。

假設：有一個函數rand()的輸出為True/False，輸出True的概率為1/2，並且各個輸出彼此獨立
結論：對於任意實數c在0到1之間，存在一個演算法輸出True/False，演算法停止的時間的期望有限，並且演算法輸出True的概率為c

可以證明任意在0到1之間的實數可以被展開為如下形式（其實就是寫成二進位）：

其中

定義函數f()：

count = 1 while (true) { x = rand() if (x &< c[count]) { return true } else if (x &> c[count]) { return false } }

現在證明
如果f輸出true，那麼只可能在第1步輸出true（設此事件為C1），第2步輸出true（C2），...，第n步輸出true（Cn)，... （彼此互斥）
如果在第i步輸出true，那麼說明x &< c[i]，並且

所以Ci的概率為（ci = 0的時候第i步輸出true的概率為0，ci=1的時候第i步輸出true的概率為1/2）

那麼求和就可以得到f()輸出true的概率為

接下來證明演算法停止時間期望有限：
在第一步停止的概率為rand()=c1的概率，也就是1/2，同理可以知道在第n步停止的概率是1/2^n，那麼演算法停止時間的期望為

其實有一個定理：用(1/2, 1/2)的Bernoulli可以模擬[0, 1]區間上的均勻分布，用的是相同的技巧，應該很多概率論的書上都有。

這就是最經典的布馮投針實驗呀，用蒙特卡洛方法得到pi相關的概率

都有一個硬幣了，那不是分分鐘造成一個隨機遊動嗎？都有隨機遊動了那也是分分鐘逼近一個布朗運動了。布朗運動都有了，分分鐘用隨機積分造出所有連續鞅。連續鞅都有了，用Feynmankac公式分分鐘都能解微分方程了

-_-|| 區區一個1/pi....

我沒那麼高的智商，不過我覺得這個方法完全不必這麼麻煩
選定硬幣上一個點，如果你的硬幣是1cm直徑這個點你就選1cm，這個點的大小可以等於直徑大小，然後規定一個方向，如果你這個點在拋擲後正好落在這個方向上就算成立，不能半包含。利用外圍周長來計算不是更好么？

如果硬幣是公平的，可以構造一列伯努利隨機變數ξi i.i.d，P(ξ1=1)=P(ξ1=-1)=1/2，記Sn=ξ1+ξ2+...+ξn，S0=0，隨機過程{Sn,n≥0}就是一維隨機遊動。

記P(2k,2n)為質點在[0,2n]時段的正方上度過2k個單位時間的概率，根據反正弦定律

令x=(sin1/2)^2，從而有質點位於正方的時間佔比不大於x的概率趨於1/π。具體操作如下：

數學不好，如有錯誤歡迎指正。
---------------------------------------
後面的內容似乎有一點偏題，如果要模擬1/π的概率，只看那個概率就好了。