利用馬可夫鏈計算擲幣事件發生的機率

本文的閱讀等級:中級

美國康乃爾大學心理學教授季洛維奇 (Thomas Gilovich) 每年都會在他的統計學課堂中安排一個實驗[1]。他要求每位學生各自寫下一組心中模擬投擲一枚公正硬幣20次所產生的隨機序列,分別以 O 和 X 代表正面和反面。但是,其中一位學生則被指派實際投擲一枚硬幣20次,也寫下他的實驗結果。季洛維奇在實驗進行前走出教室,等他返回教室後,他將接受一項挑戰:檢視所有學生繳交的實驗記錄,然後判斷其中那一張紙記載了實際擲幣產生的序列。季洛維奇總是能令學生們驚訝不已,他無一次例外地揀選出真實的擲幣序列。究竟他是怎麼辦到的?季洛維奇既沒有暗藏機關也不具特異能力,他掌握的技能不過就是「資訊不對稱」。身為心理學教授,他知道絕大多數人──包括教室中的學生──總是低估了出現連續正面或反面的機率。真實的擲幣結果幾乎都是那張記錄著最長的連續正面或反面的序列,例如:

OXXXXXOXOOXOOXOOXOOX,

而學生們想像出來的擲幣序列則經常如下:

XXOXOOOXOOXOXXOOXXOO。

本文的主題即在破解季洛維奇的戲法:釐清投擲一枚公正硬幣 n 次,計算出現至少連續 m 次正面的機率。

 
考慮最簡單的情況:m=2,即投擲一枚硬幣 n 次,問出現至少連續2次正面的機率為何?明顯地,組合數學是最直接的解法,也就是計算 2^n 種可能結果中包含至少連續2次正面的序列個數。不過這個方法需要縝密的分析與推理,有興趣深入了解的讀者請參閱文獻[2,3]。下面我介紹一個線性代數解法,主要的想法是利用馬可夫鏈 (Markov chain)[4]來表示擲幣的隨機過程 (見下圖)。根據問題設定,我們需要三種狀態:S_1 (反面,X),S_2 (正面,O),S_3 (連續兩次正面,OO)。若第一次擲幣結果為反面,則進入 S_1,否則進入 S_2,往後的每次擲幣結果決定如何轉移至下一個狀態,所有可能的情況如下:

\begin{aligned} S_1&\xrightarrow{X}S_1,~S_1\xrightarrow{O}S_2\\ S_2&\xrightarrow{X}S_1,~S_2\xrightarrow{O}S_3\\ S_3&\xrightarrow{X}S_3,~S_3\xrightarrow{O}S_3\end{aligned}

注意,一旦抵達 S_3,即表示已經獲得連續 2 次正面,因此不論之後的投擲結果如何,我們仍然處於此狀態。譬如,對應 XOXXOOX 的狀態序為

S_1\to S_2\to S_1\to S_1\to S_2\to S_3\to S_3

令正面出現的機率為 p,反面出現的機率為 q=1-p。又令 P_i(k) 表示投擲完第 k 次後,k\ge 1,處在狀態 S_i 的機率。因為對於所有 kP_1(k)+P_2(k)+P_3(k)=1,我們的目標是計算 P_3(n)1-P_1(n)-P_2(n)

Markov chain: Two heads in a row

 
根據上圖寫出狀態轉移機率方程組:

\begin{aligned} P_1(k)&=qP_1(k-1)+qP_2(k-1)\\ P_2(k)&=pP_1(n-1)\\  P_3(k)&=pP_2(n-1)+P_3(n-1),\end{aligned}

並有初始條件:P_1(1)=q, P_2(1)=p, P_3(1)=0 (投擲一次硬幣不可能產生 OO)。將線性方程組以矩陣形式表示如下:

\begin{bmatrix} P_1(k)\\ P_2(k)\\ P_3(k) \end{bmatrix}=\begin{bmatrix} q&q&0\\  p&0&0\\  0&p&1  \end{bmatrix}\begin{bmatrix}  P_1(k-1)\\ P_2(k-1)\\ P_3(k-1)  \end{bmatrix}

也就得到一個典型的差分方程 \mathbf{u}_k=A\mathbf{u}_{k-1},其中

\mathbf{u}_k=\begin{bmatrix} P_1(k)\\ P_2(k)\\ P_3(k) \end{bmatrix},~A=\begin{bmatrix}  q&q&0\\  p&0&0\\  0&p&1  \end{bmatrix}

轉移矩陣 (transition matrix) A 也稱為馬可夫矩陣,其中所有元皆不為負,而且每一行的所有元總和為 1。關於馬可夫矩陣的基本性質,請見“馬可夫過程”。

 
再來是解差分方程的計算工作。令 A 的特徵值為 \lambda_1,\lambda_2,\lambda_3,對應的特徵向量為 \mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3。若 A 有相異的特徵值,則 A 可對角化,\mathbf{u}_k=A\mathbf{u}_{k-1} 的解即為

\mathbf{u}_k=c_1\lambda_1^k\mathbf{x}_1+c_2\lambda_2^k\mathbf{x}_2+c_3\lambda_3^k\mathbf{x}_3

寫出矩陣 A 的特徵多項式:

\det(A-\lambda I)=\begin{vmatrix} q-\lambda&q&0\\ p&-\lambda&0\\ 0&p&1-\lambda  \end{vmatrix}=(1-\lambda)(\lambda^2-q\lambda-pq)

解得

\displaystyle\lambda_1=\frac{q+r}{2},~\lambda_2=\frac{q-r}{2},~\lambda_3=1

其中 r=\sqrt{q^2+4pq}。注意,只要 pq 不等於 0\lambda_1,\lambda_2,\lambda_3 必定相異且 r\neq 0。另外,\lambda_1+\lambda_2=q\lambda_1-\lambda_2=r\lambda_1\lambda_2=-pq,這些關係式將於以下化簡運算中使用。下一個步驟是解出零空間 N(A-\lambda_iI) 基底,i=1,2,3,此即對應的特徵向量。請讀者自行計算,結果如下:

\mathbf{x}_1=\begin{bmatrix} \lambda_1\\ p\\  p^2/(\lambda_1-1)  \end{bmatrix},~\mathbf{x}_2=\begin{bmatrix}  \lambda_2\\  p\\  p^2/(\lambda_2-1)  \end{bmatrix},~\mathbf{x}_3=\begin{bmatrix}  0\\  0\\  1  \end{bmatrix}

於是得到

\mathbf{u}_k=c_1\lambda_1^k\begin{bmatrix} \lambda_1\\ p\\  p^2/(\lambda_1-1)  \end{bmatrix}+c_2\lambda_2^k\begin{bmatrix}  \lambda_2\\  p\\  p^2/(\lambda_2-1)  \end{bmatrix}+c_3\begin{bmatrix}  0\\  0\\  1  \end{bmatrix}

最後還要代入初始條件解出係數 c_1,c_2,c_3。考慮 k=1,將 \mathbf{u}_1=\begin{bmatrix}  q\\  p\\  0  \end{bmatrix} 代入上式,若 p\neq 0 (這符合一般情況),就有

\displaystyle\begin{aligned} c_1\lambda_1^2+c_2\lambda_2^2&=q\\ c_1\lambda_1+c_2\lambda_2&=1\\  c_1p^2\frac{\lambda_1}{\lambda_1-1}+c_2p^2\frac{\lambda_2}{\lambda_2-1}+c_3&=0,\end{aligned}

第二式通乘 \lambda_2,與第一式相減,可得

\displaystyle c_1=\frac{q-\lambda_2}{\lambda_1(\lambda_1-\lambda_2)}=\frac{\lambda_1}{\lambda_1r}=\frac{1}{r}

此結果代回第二式,立得

\displaystyle c_2=\frac{1}{\lambda_2}\left(1-\frac{\lambda_1}{r}\right)=\frac{r-\lambda_1}{\lambda_2r}=\frac{-\lambda_2}{\lambda_2r}=-\frac{1}{r}

再將得到的 c_1c_2 代入第三式,

\displaystyle\begin{aligned} c_3&=\frac{p^2}{r}\left(\frac{\lambda_1}{1-\lambda_1}-\frac{\lambda_2}{1-\lambda_2}\right)\\ &=\frac{p^2}{r}\cdot\frac{\lambda_1-\lambda_2}{(1-\lambda_1)(1-\lambda_2)}\\ &=\frac{p^2}{1-\lambda_1-\lambda_2+\lambda_1\lambda_2}\\ &=\frac{p^2}{1-q-pq}=\frac{p^2}{p-pq}\\ &=\frac{p^2}{p(1-q)}=\frac{p}{1-q}=1\end{aligned}

最後我們將三個狀態的發生機率整理於下:

\displaystyle\begin{aligned} P_1(k)&=\frac{1}{r}\left(\lambda_1^{k+1}-\lambda_2^{k+1}\right)\\ P_2(k)&=\frac{p}{r}\left(\lambda_1^k-\lambda_2^k\right)\\  P_3(k)&=\frac{p^2}{r}\left(\frac{\lambda_1^k}{\lambda_1-1}-\frac{\lambda_2^k}{\lambda_2-1}\right)+1\end{aligned}

 
觀察發現 P_1(k)P_2(k) 的形式雷同並且比 P_3(k) 單純,因此我們選擇計算 1-(P_1(n)+P_2(n))。考慮公正硬幣,p=q=1/2,則 r=\sqrt{5}/2,且

\displaystyle \lambda_1=\frac{1}{2}\left(\frac{1+\sqrt{5}}{2}\right),~\lambda_2=\frac{1}{2}\left(\frac{1-\sqrt{5}}{2}\right)

代入數值化簡整理後可導出

\displaystyle\begin{aligned} P_1(n)&=\frac{1}{2^n}\cdot\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^{n+1}-\left(\frac{1-\sqrt{5}}{2}\right)^{n+1}\right]=\frac{F_{n+1}}{2^n}\\ P_2(n)&=\frac{1}{2^n}\cdot\frac{1}{\sqrt{5}}\left[\left(\frac{1+\sqrt{5}}{2}\right)^{n}-\left(\frac{1-\sqrt{5}}{2}\right)^{n}\right]=\frac{F_n}{2^n},\end{aligned}

其中 F_n 表示費布納西數列的第 n 項 (見“費布納西數列的表達式”)。利用費布納西數列的遞歸生成方程式:

F_{n+2}=F_{n+1}+F_{n}

即得到一個非常簡潔的公式:

\displaystyle P_3(n)=1-(P_1(n)+P_2(n))=1-\frac{F_{n+2}}{2^n}

換句話說,投擲一枚公正硬幣 n 次,不出現連續2次以上正面的序列總數就是費布納西數 F_{n+2}

 
下一步要如何將 m=2 的結果推廣至一般情況?對於任意 2<m\le n,我們可以重複以上步驟。默想 m=3 的情況,前述馬可夫鏈模型將增加一個新狀態 S_4 以表示 OOO,狀態轉移矩陣 A 則擴大為 4\times 4 階,再經過一番努力算出 P_1(n)+P_2(n)+P_3(n)。但如果我們相信數學的基本形式無所不在,這時應該可以推測投擲一枚公正硬幣 n 次,出現至少連續 m>2 次正面的機率也有和 m=2 相仿的計算公式。確實如此,下為一般公式:

\displaystyle 1-\frac{F^{(m)}_{n+2}}{2^n}

其中廣義費布納西數 F^{(m)}_{n} 等於之前的 m 個 (廣義) 費布納西數之和,並定義 F^{(m)}_{n}\equiv 0,若 n\le 0。譬如,當 m=3,廣義費布納西數列 F^{(3)}_{n}n=0,1,2,\ldots, 的最初幾項是

0, 1, 1, 2, 4, 7, 13, 24, 44, 81, 149,…。

讀者如果懷疑上式的正確性,不妨拿些例子來驗證。若擲幣5次,則出現至少連續3次正面的機率為 1-{F_7^{(3)}}/{2^5}=1-24/32=8/32,8種可能情況如下:

OOOXO, OOOXX, XOOOX, OXOOO, XXOOO, OOOOX, XOOOO, OOOOO。

 
回到本文開頭講述的故事,投擲一枚公正硬幣20次,我們好奇:出現至少連續4次,甚至5次正面的機率有多大?當 n=20m=4,結果是

\displaystyle 1-\frac{F^{(4)}_{22}}{2^{20}}=1-\frac{547,337}{1,048,576}=1-0.5220=0.4780

n=20m=5,結果是

\displaystyle 1-\frac{F^{(5)}_{22}}{2^{20}}=1-\frac{786,568}{1,048,576}=1-0.7501=0.2499

擲幣20次中出現至少連續5次正面的機率與至少連續5次反面的機率相同,大約是1/4,扣除至少連續5次正面與連續5次反面同時發生的機率,真實實驗中出現至少連續5次正面或反面的機率小於 (但接近) 1/2。如此說來,Gilovich 的勝算似乎不太高,但為什麼每年他總能精準命中呢?可能的解釋有二:第一,擲幣20次中出現至少連續4次正面或反面的機率介於 0.4780 與 0.9560 之間 (扣除兩者同時發生的機率);第二,參與模擬實驗的學生們對於出現連續正面或反面的機率估計有顯著的偏差。或許很多人心想擲幣4次出現 OOOO 的機率僅為1/16,致使他們在寫下 OOO 之後,「忍不住」又補上了 X。所以可以這麼說,擲幣20次中出現連續4次以上相同結果的發生率大於50% (這還是保守的數字),以及全體學生有志一同的反向操作,這兩個因素相輔相成才得以成就季洛維奇博得滿堂彩的精湛演出。

 
附註:利用遞歸法也可以計算投擲一枚硬幣 n 次出現至少連續 m 次正面的機率。以 m=2 為例,令 P_k 表示擲幣 k 次出現至少連續2次正面的機率,Q_k 表示擲幣 k 次不發生連續2次正面的機率,故 P_k=1-Q_k。下面我們以遞歸法計算 Q_k。明顯地,Q_1=1Q_2=3/4。當 k>2 時,考慮二種可能情況:(1) 若第一次擲出反面,爾後的 k-1 次投擲不出現連續2次正面的機率是 Q_{k-1};(2) 若第一次擲出正面,則第二次必須擲出反面 (否則將出現連續2次正面),剩下的 k-2 次投擲不出現連續2次正面的機率是 Q_{k-2}。因為第一次擲出反面的機率是 1/2,第一次擲出正面且第二次擲出反面的機率是 1/4,根據上述分析 Q_k 可表示為

\displaystyle Q_k=\frac{1}{2}Q_{k-1}+\frac{1}{4}Q_{k-2},~~k>2

上式等號兩邊同乘以 2^k,可得

\displaystyle 2^kQ_k=2^{k-1}Q_{k-1}+2^{k-2}Q_{k-2},~~k>2

R_k=2^kQ_k,就有

R_k=R_{k-1}+R_{k-2}

這剛好就是費布納西數列的遞歸關係式。因為 R_k 的初始狀態是 R_1=2R_2=3,即知 R_k=F_{k+2},則 Q_k=F_{k+2}/2^k,故投擲一枚硬幣 n 次出現至少連續2次正面的機率是

\displaystyle P_n=1-\frac{F_{n+2}}{2^n}

使用相同方式不難推論擲幣 n 次出現至少連續 m 次正面的機率是

\displaystyle P_n=1-\frac{F^{(m)}_{n+2}}{2^n}

其中 F^{(m)}_n 是廣義費布納西數,如上文所述。

 
參考來源:
[1] Gary Belsky and Thomas Gilovich, Why smart people make big money mistakes and how to correct them, Fireside, pp 119, 1999.
[2] http://elec424.wikia.com/wiki/Two_Heads_in_a_Row_in_n_Tosses
[3] http://marknelson.us/2011/01/17/20-heads-in-a-row-what-are-the-odds/
[4] http://en.wikipedia.org/wiki/Markov_chain

廣告
本篇發表於 機率統計 並標籤為 , , , , , , 。將永久鏈結加入書籤。

3 Responses to 利用馬可夫鏈計算擲幣事件發生的機率

  1. pgcci7339 說道:

    老師您好~對於廣義費布納西數的定義,我可能還不是看的非常清楚~想請問您一下:)

    當m=3時,F^{(3)}_{n} 的前三項 0,1,1 是如何來的呢,第一個0我知道根據定義,那後面兩個1呢?
    若m=4,前4個是否是0,1,1,2?
    若m=5,前5個是否為0,1,1,2,3?

發表迴響

在下方填入你的資料或按右方圖示以社群網站登入:

WordPress.com Logo

您的留言將使用 WordPress.com 帳號。 登出 / 變更 )

Twitter picture

您的留言將使用 Twitter 帳號。 登出 / 變更 )

Facebook照片

您的留言將使用 Facebook 帳號。 登出 / 變更 )

Google+ photo

您的留言將使用 Google+ 帳號。 登出 / 變更 )

連結到 %s