PS6 - 公正不公正? 不確定性的度量: 貝氏推論 (上)
 |
| PASIEKA / SCIENCE PHOTO LIBRARY |
Data analysis is simply a dialogue with the data
—S. F. Gull
一、機率到底告訴我們了什麼?
假定我們有兩枚六面骰A和B每面為標記為1~6的整數，我們如何確定兩枚骰子都是公正的? 即每個面骰出的機率皆為1/6。古典機率的概念告訴我們若擲兩枚骰子多次，並記錄每面出現的次數，如果一枚骰子是公正的，那麼理論上每面出現的次數相同或接近相同。
 |
| 圖1. 骰A (粉藍) 與骰B (粉紅) 各擲50次每面出現的次數 |
現在我擲50次並記錄兩枚骰子每個面出現的次數，並用將其畫成直方圖如圖1所示，從記錄的結果來看，我可以斷定哪個骰子是公正的嗎? 看起來好像不行，那如果我改成擲100次、1000次就可以分別公不公正了嗎? 換個方向來看，也許我們不在乎骰子是否公正，我們在乎的是用我們觀測的資料集\(D\)，它是否可以告訴我們骰子每面出現的機率分佈是什麼? 都是同一組資料，一個是基於古典推論 (classical inference)，另一個則是貝氏推論 (Bayesian inference)。對於古典推論，即頻率論而言，我們有以下原則 (tenets)：
- 機率對應到的是事件出現的頻率，它們是真實世界的客觀呈現 (objective properties)
- 參數是固定的 (未測量前是未知的)，比如六面骰每面擲出的機率可以是這裡提到的參數。這裡的參數對應的是柏拉圖世界的量，如果是一枚六面公正骰，那麼每面擲出的機率一定是1/6，不可能有波動 (fluctuation)
- 當我們重複步驟愈多次 (如擲骰的次數)，那麼每面出現的頻率應該是一個可以被定義的量 (即收斂，如果該量會產生波動則不收斂)
- 機率是對主觀信念 (subjective beliefs, prior) 的描述，而不受限於資料出現的頻率長遠來看是否是一個定值
- 我們對於模型參數 (例如每個面出現的機率) 的推論給出的是它位於“某個值”的機率分佈，這個分佈代表的是我們對這個參數有多清楚 (不清楚的話該參數就會是一個很大的範圍)
以我們現有的資料，我們只能給出這樣的結果，我沒有辦法100%肯定他一定是公正的，但我可以給出它會偏離公正骰有多遠，一個信心區間 (confidence level, C.L.)。當然透過累積更多的觀測資料，我們可以不斷更新我們的posterior而縮小這個不確定性。所以貝氏推論描述的是基於我們所擁有的知識 (即資料\(D\)) 能給出推定的範圍，是知識論 (epistemology) 的上限。
在我們往下以前，先揭露骰A確實是公正骰，但骰B每個面的權重是\(p=\{\frac{1}{12},\frac{1}{6},\frac{1}{4},\frac{1}{6},\frac{1}{4},\frac{1}{12}\}\)，只是骰的次數不夠多次無法直觀地看出這個差異。
二、最大似然估計法
對於參數的估計，我們可以採用最大似然估計法 (maximum likelihood estimation, MLE)1，這個方法在古典推論中被大量的使用。一開始我們先談談什麼是似然函數\(\mathcal{L}\) (likelihood function)。假設我們有一組資料集\(D=\{x_i\}\)其中集內的每個樣本\(x_i\)都是對某個未知分佈\(M(\boldsymbol{\theta})\)的取樣，而\(\boldsymbol{\theta}\)是\(M\)的參數。所以當分佈\(M\)給定後，它產生每個樣本\(x_i\)的機率分佈 (pdf) 就是\(p(x_i|M(\boldsymbol{\theta}))\)，或稱作樣本\(x_i\)的datum likelihood。若有\(n\)個樣本，那麼就有\(n\)個data likelihood：\(p(x_1|M(\boldsymbol{\theta}))\)、\(p(x_2|M(\boldsymbol{\theta}))\)⋯\(p(x_n|M(\boldsymbol{\theta}))\)，而\(\mathcal{L}\)就是所有這些\(p\)的積 (product)，即\[\mathcal{L}\equiv p(\{x_i\}|M(\boldsymbol{\theta}))=\prod_{i=1}^{n} p(x_i|M(\boldsymbol{\theta})),\tag{6.1} \]其中一個\(M\)可能有很多個參數\(\boldsymbol{\theta}=\{\theta_1,\theta_2,...,\theta_k\}\)，若\(M\)是Gaussian分布，則稱Eq. (6.1)為Gaussian likelihood。而MLE即找出一組\(\boldsymbol{\theta}\)讓\(\mathcal{L}\)有最大值，該組\(\boldsymbol{\theta}\)就是我們要的參數。它的想法很簡單，每個datum likelihood都表示模型\(M(\boldsymbol{\theta})\)產生隨機變數\(x_i\)的pdf，最有可能的參數\(\boldsymbol{\theta}\)表示可以讓整組資料\(D=\{x_i\}\)被\(M\)生出來的機率最大，故所有data likelihood的乘積\(\Pi_{i=1}^n p(x_i|M)\)也會是最大，這就是MLE的概念。找出\(\mathcal{L}\)極值的做法便是其對\(\boldsymbol{\theta}\)的微分等於0處，由於有k個\(\theta\)，所以\[\left.\frac{\partial\mathcal{L}}{\partial\theta_{k}}\right|_{\theta_{k}=\theta_{0}}=0.\tag{6.2}\]由於\(\mathcal{L}\)是由許多組pdf相乘所構成，為了計算方便可以將其取ln，所有相乘的項就變相加了，計算因而簡單許多，對應的對數似然函數\(\ln \mathcal{L}\) (log-likelihood function)\[\left.\frac{\partial\ln \mathcal{L}}{\partial\theta_{k}}\right|_{\theta_{k}=\theta_{0}}=0.\tag{6.2}\]假定\(M\)是高斯分佈的情形，參數\(\boldsymbol{\theta}=\{\mu\}\)。對資料集\(D=\{x_i\}\)，透過Eq. (6.2)可得解析解在\[\mu=\mu_0=\frac{\sum_i^n w_i x_i}{\sum_i^n wi}\tag{6.3}\]有MLE極大值。其中\(w_i=\sigma_i^{-2}=(\sigma_o^2+e_i^2)^{-1}\)，而\(\sigma_o\)是樣本分佈的真實標準差、\(e_i\)是每次量測時的誤差。另外\(\mu_0\)的誤差可以表示成\[\sigma=\sigma_\mu=\left(\sum_{i=1}^n w_i \right)^{-1/2}.\tag{6.4}\] 另一個問題就是對同一組\(D\)而言，我們當然可以假設M是高斯分佈，但當然也可能不是，這時候我們就需要做擬合度檢驗 (goodness of fit)。若是高斯，則我們有Gaussian likelihood\[\ln \mathcal{L}={\rm const.}-\frac{1}{2}\sum_{i=1}^n z_i^2={\rm const.}-\frac{1}{2}\chi^2\]其中\(z_i=(x_i-\mu)/\sigma\)。所以\(\mathcal{L}\)可以視為是一個擁有\(N-k\)自由度的卡方分佈，其中\(k\)是模型參數的數量。對於一個擬合度良好的模型，我們期望\[\chi^2_{\rm dof}=\frac{1}{N-k}\sum_{i=1}^n z_i^2\approx 1,\]但若\(\chi^2_{\rm dof}-1\)遠大於\(\sqrt{2/(N-k)}\)，則不太可能\(D\)是來自該高斯分佈。除了卡方檢定外，還有其他檢驗也可以用來測試\(D\)是否是來自高斯分佈，例如Anderson-Darlin檢定和Kolmogorov-Smirnov檢定，這裡就不贅述了。
1 注意maximum likelihood estimation和maximum likelihood estimate的縮寫都是MLE，前者指的是估計模型參數的方法，後者指的是所有的data likelihood裡面最大的那個datum likelihood: \(\max_{x_{i}\in D} p(x_i|M)\)
三、貝氏推論
接下來會很常舉P&S5中疾病診斷判定是否病人得的是天花的例子，忘記了可以再去翻翻看複習一下。在上一個小節中我們簡述了古典推論透過MLE來給定參數\(\boldsymbol{\theta}\)的值，另外可以注意到的是古典推論討論的是\(p(D|M)\)，即給定\(M\)下，它產生\(D=\{x_i\}\)的機率。所以不管怎樣我一定先假定了分佈\(M\)了，可是分佈這麼多我能確定一定是這個\(M\)嗎? 比如說是高斯分佈，但如果檢驗告訴我他不是高斯，那又是什麼? 並不是每種狀況都可以單純從\(D\)中便知道他是什麼分佈，我們更希望的是\(D\)是否可以自己說出多有可能它是來自某個分布\(M\)！所以我們想知道的便是\(p(M|D)\)，也許一開始我們並不知道真正的\(M\)是什麼，可能肇因於資料不夠多，但隨著愈來愈多的觀測資料被加進\(D\)裡，我們同時也不斷地更新\(p(M|D)\)。我們可以把這類對貝氏定理的理解稱為歷時詮釋 (diachronic interpretation)，也是P&S5診斷天花例子的解釋，為什麼醫生最後判定帶痘斑來的病人是天花的機率變高了。利用貝氏定理，我們可以得到\[p(M,\boldsymbol{\theta}|D)=\frac{p(D|M,\boldsymbol{\theta})p(M,\boldsymbol{\theta})}{p(D)}\tag{6.5}\]其中每一項的含義為：
- \(p(M,\boldsymbol{\theta}|D)\)即posterior，是我們想知道跟\(D\)有關的分布的機率
- \(p(D|M,\boldsymbol{\theta})\)即data likelihood
- \(p(M,\boldsymbol{\theta})\)即prior，是在得到或更新觀測資料\(D\)前我們對\(M\)的主觀假定
- \(p(D)\)是一歸一化常量 (normalization constant)，又稱標準化常量，是所有任何可能得到\(D\)的\(M\)的聯合機率分布
1. 如何選擇prior?
由於除了我們要分析的資料外沒有其他關於prior的資訊了，所以初始prior的選擇可以基於一些規則或經驗，例如說“按照模型參數所給出的變異數必須是正定的”，不過變異數本來就是二次方的矩 (moment)，不管模型參數多麼稀奇古怪，它必然是正定的；而診斷例子中的\(p(天花)=10^{-4}\)，只能說good doctor認為不失一般性下染上天花的機率很小而已，本質上這些prior好像說了什麼，細看又什麼都沒說，這些prior就稱作無信息先驗機率 (uninformative prior)。這裡的無信息表明了這些prior盡可能是客觀不帶有偏見地去影響posterior的結果。舉例像是統派會說：“長遠來看兩岸合久必分分久必合，統一是必然的”，基於這個prior我們可以推論出一個posterior那就是兩岸必然統一，但什麼時候會統一，這個prior的資訊是完全沒有畫面的！ (它只是客觀地講出他認為最終的可能，但沒法告訴你是明天、下個月或是明年會統一)2 這裡還是要強調一下，即便是最無信息的prior，它仍是我們計算posterior的一部分 (權重)，所以它仍然會造成一定的影響，而貝氏推論最後也不必然等於古典推論的結果，而最大的那個datum likelihood (maximum likelihood estimate) 也不一定影響posterior最多，它可能因為得到比較小的prior而效應變差 (像\(p(天花|痘斑)\)雖然是MLE，但\(p(天花)\)小很多，最後也壓低了病人得的是天花的後驗機率)。
再來看回六面骰的例子了，prior是每個面被骰出的機率\(p_i\)的離散集合，所以我們可以導出兩個約束條件 (constraints)，第一個是歸一化條件\[\sum_{i=1}^6 p_i=1\tag{6.7}\]而第二個是點數的期望值\(\mu\)\[\sum_{i=1}^6 ip_i=\mu.\tag{6.8}\]為了找出\(S\)在上面這兩個條件成立下的極值，可以直接套用拉格朗日乘子法 (Lagrange's multipliers)\[Q=S+\lambda_0\left( 1-\sum_{i=1}^6 p_i\right)+\lambda_1\left(\mu-\sum_{i=1}^6 ip_i\right)\]其中\(S\)為Eq. (6.6)。所以極值會發生在\(\partial Q/\partial p_i=0\)的地方 (多變數就是找梯度為0的點)，對這個六面骰的例子，我們有\(p_i=\exp(-1-\lambda_0)\exp(i \lambda_1)\)。雖然有兩個未知數\(\lambda_0\)和\(\lambda_1\)，但也有兩個約束Eq. (6.7)和Eq. (6.8)，就可以得到\[p(\theta)=\frac{1}{\mu}\exp \left(-\frac{\theta}{\mu}\right),\tag{6.9}\]而\(\theta\)就是該面的點數。我們便可以將Eq. (6.9)當作求六面骰每面被骰出的機率的prior。
2 一個例子像是隨便假設flat prior: \(p(\theta)\propto C\)且\(C>0\)是一常量，因為在\(\theta\)空間這個prior是發散的\(\int_{-\infty}^\infty p(\theta) d\theta \to \infty\)，這種prior也稱作improper prior。不過只要posterior是well-defined，prior是不是improper沒有關係，另外我們也可以設定\(\theta\)必須處在一個區間，這樣他就不會發散了
給定了\(D\)，我們心中有兩種可能描述\(D\)的說法，即\(M_1\)和\(M_2\)，至於誰是最有可能的說法呢? 則可以透過勝算比 (odds ratio，一說幾率，不是機率) 來估算，定義為\[Q_{21}=\frac{p(M_2|D)}{p(M_1|D)}.\tag{6.10}\]通常來說：
1.1 熱力學觀點：夏儂熵 (Shannon's entropy)
為了盡量減少prior對計算posterior產生的bias，我們當然希望prior愈uninformative愈好，我們來看看如何減低prior所帶的資訊。以下參考自ref. 5。假設有一組pdf含\(n\)個離散值\(p_i\) (像六面骰\(n=6\)，每一面骰到的機率就是\(p({\rm face})\))，歸一化保證了\(\sum_{i=1}^n p_i=1\)，我們便可以定義\(S\)為該分佈的夏儂熵 (Shannon's entropy)\[S=-\sum_{i=1}^n p_i \ln p_i.\tag{6.6}\]另從L'Hôpital's rule我們有\(\lim_{p\to 0} p\ln p=0\)，順便補充當pdf是連續而不是離散的話，Shannon's entropy可以用積分表示\[S=-\int_{-\infty}^\infty p(x)\left(\frac{p(x)}{m(x)}\right)dx,\]其中\(m(x)\)稱為測度 (measure)，當我們對\(x\)做變數變換時，測度保證對應機率的不變性，可以視為是Jacobian。根據熱力學定律，當熵愈大這個系統就愈無序，亦即有意義的資訊愈少，所以我們想要的愈uninformative prior就是可以讓\(S\)有最大的情況，其隱含了最少量的資訊。再來看回六面骰的例子了，prior是每個面被骰出的機率\(p_i\)的離散集合，所以我們可以導出兩個約束條件 (constraints)，第一個是歸一化條件\[\sum_{i=1}^6 p_i=1\tag{6.7}\]而第二個是點數的期望值\(\mu\)\[\sum_{i=1}^6 ip_i=\mu.\tag{6.8}\]為了找出\(S\)在上面這兩個條件成立下的極值，可以直接套用拉格朗日乘子法 (Lagrange's multipliers)\[Q=S+\lambda_0\left( 1-\sum_{i=1}^6 p_i\right)+\lambda_1\left(\mu-\sum_{i=1}^6 ip_i\right)\]其中\(S\)為Eq. (6.6)。所以極值會發生在\(\partial Q/\partial p_i=0\)的地方 (多變數就是找梯度為0的點)，對這個六面骰的例子，我們有\(p_i=\exp(-1-\lambda_0)\exp(i \lambda_1)\)。雖然有兩個未知數\(\lambda_0\)和\(\lambda_1\)，但也有兩個約束Eq. (6.7)和Eq. (6.8)，就可以得到\[p(\theta)=\frac{1}{\mu}\exp \left(-\frac{\theta}{\mu}\right),\tag{6.9}\]而\(\theta\)就是該面的點數。我們便可以將Eq. (6.9)當作求六面骰每面被骰出的機率的prior。
2 一個例子像是隨便假設flat prior: \(p(\theta)\propto C\)且\(C>0\)是一常量，因為在\(\theta\)空間這個prior是發散的\(\int_{-\infty}^\infty p(\theta) d\theta \to \infty\)，這種prior也稱作improper prior。不過只要posterior是well-defined，prior是不是improper沒有關係，另外我們也可以設定\(\theta\)必須處在一個區間，這樣他就不會發散了
2. 模型M的選擇
就像天花的例子一樣，如果只有一種模型 (沒有水痘可比)，就沒有什麼必要使用貝氏推論，長斑痘無疑是染上天花了。需要貝氏推論的地方在於對同一種狀況 (觀測資料\(D\))，我們可能有多個模型都可以符合，但基於每個模型的主觀先驗 (prior)，我們希望可以得出哪個模型符合觀測的勝算最大。給定了\(D\)，我們心中有兩種可能描述\(D\)的說法，即\(M_1\)和\(M_2\)，至於誰是最有可能的說法呢? 則可以透過勝算比 (odds ratio，一說幾率，不是機率) 來估算，定義為\[Q_{21}=\frac{p(M_2|D)}{p(M_1|D)}.\tag{6.10}\]通常來說：
- \(Q_{21}>100\)：資料支持\(M_2\)的證據是決定性的
- \(10\lesssim Q_{21}<100\)：資料強烈支持\(M_2\)
- \(Q_{21}<3\)：資料對\(M_2\)的支持只能算勉勉強強，基本上兩個模型不分軒輊
2.1 再一個例子：拋硬幣
數學總是很抽象，我們舉的例子來輔助。我手邊有一枚硬幣，拋了\(N\)次後要拋出幾次正面才能說硬幣是公正的?
\(N=10\)
|
\(N=20\)
|
|
正面
|
3
|
11
|
反面
|
7
|
9
|
從表1中我給出兩個假設，一個是硬幣是公正的\(M_1\)，拋出正面的機率是\(b_1=0.5\)，而從拋10次的結果來看，我覺得硬幣比較偏向反面，對\(M_2\)而言我可以認為\(b_2=0.3\)，另外從資料來看我已經假定了硬幣是不公正的，我可以選擇讓\(p(M_2)\)這個prior比\(p(M_1)\)大一點，假定跟拋10次的結果一樣是7-3比。另外拋硬幣是個二項式檢定，它的data likelihood可以表示成\[p(k|M_i)=C^N_k b_i^k (1-b_i)^{N-k}\]其中\(k\)是拋出正面的次數。對拋10次而言，我有勝算比 (二項式係數在這邊都相互消去了)\[Q_{21}=\frac{0.3^3(1-0.3)^7\times 7}{0.5^3(1-0.5)^7\times 3}\approx 5.31.\]所以從拋10次的資料來看，加上我的偏見，有一定的勝算顯示硬幣是不公正的，不過由於\(Q_{21}\)沒有大於10，我還不敢斷定。為了更加確定，我又多拋了幾次，現在我的資料\(D\)來到了拋出20次，結果也如表1所示。我可以依照新的資料更新勝算比得到\[Q_{21}=\frac{0.3^{11}(1-0.3)^{9}\times 7}{0.5^{11}(1-0.5)^9\times 3}\approx 0.17.\]哇，翻盤！看起來新的資料favor硬幣是公正的，即使我讓不公正的權重一開始比較重 (7-3比)。所以從這個例子也可以看到posterior隨著資料的更新而有變化，這是歷時詮釋的另一個驗證。
拋硬幣的例子算是單純，真實世界的問題可能無法這麼快就猜出data likelihood，另外不同模型的參數數量也可能不一樣，如果勝算比顯示兩個模型解釋\(D\)的能力不分軒輊，一個簡單判定應該選擇哪個模型的方法是利用奧坎剃刀 (Occam's razor) 的原則：較少即是較好，我們會選擇參數數量比較少的那個模型作為答案。接下來我想要談的是用貝氏推論做參數的估計，我們會使用到MLE的方法，不過由於再寫下去篇幅就太大了，所以我們把P&S6拆成兩部分，剩下的留到下篇再討論囉！
留言
張貼留言