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)
• 當我們重複步驟愈多次 (如擲骰的次數)，那麼每面出現的頻率應該是一個可以被定義的量 (即收斂，如果該量會產生波動則不收斂)

但對於貝氏推論而言，它則是基於如下的態度 (stance)：

• 機率是對主觀信念 (subjective beliefs, prior) 的描述，而不受限於資料出現的頻率長遠來看是否是一個定值
• 我們對於模型參數 (例如每個面出現的機率) 的推論給出的是它位於“某個值”的機率分佈，這個分佈代表的是我們對這個參數有多清楚 (不清楚的話該參數就會是一個很大的範圍)

關於古典推論的原則是沒有問題的，非黑即白，但對於貝氏推論的部分，我們應該多做點說明。所謂的主觀信念，就是我們“認為”這個骰是公正 (或不公正) 的。假定是公正的，那麼骰出面是1的機率是$p(1)=1/6\approx 0.167$，基於這樣的prior搭上現有的資料，透過貝氏推論，我們可以給出一個posterior假設是95% C.L.的$p(1)=0.167\pm 0.04$區間。這個意思大概是從資料來看，如果我們找100個人來擲這枚骰子每個人都擲了1000次，其中約有95人會得到骰子骰出面1的次數是介於127和207次之間 (完美的公正骰約是167次)。
　以我們現有的資料，我們只能給出這樣的結果，我沒有辦法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)¹，這個方法在古典推論中被大量的使用。一開始我們先談談什麼是似然函數$\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檢定，這裡就不贅述了。

---

¹ 注意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}$$ 其中每一項的含義為：

1. $p(M,\boldsymbol{\theta}|D)$即posterior，是我們想知道跟$D$有關的分布的機率
2. $p(D|M,\boldsymbol{\theta})$即data likelihood
3. $p(M,\boldsymbol{\theta})$即prior，是在得到或更新觀測資料$D$前我們對$M$的主觀假定
4. $p(D)$是一歸一化常量 (normalization constant)，又稱標準化常量，是所有任何可能得到$D$的$M$的聯合機率分布

注意這裡我們不再把$\boldsymbol{\theta}$寫成$M(\boldsymbol{\theta})$，由於$M$不再是一個已知的分佈，不像古典推論一樣必須先給定比如Gaussian分佈，所以可以知道有什麼參數$\boldsymbol{\theta}$。在這裡我們只能說有一個分佈$M$以及參數$\boldsymbol{\theta}$，在給定$D$的情況下，他們最有可能的分佈區間為何。另外在最一般的情況下，我們根本無法確切知道有哪些$M$會有非0的機率去產生$D$，所以$p(D)$本身的實際意義並不明確，正如同P&S5診斷例子中的$p(痘斑)$，無法確定所有可能產生痘斑的疾病在所有人口中的機率分佈，更遑論求出確切的值！不過我們稍後也會看到是否知道$p(D)$並不影響我們的求解。

1. 如何選擇prior?　

由於除了我們要分析的資料外沒有其他關於prior的資訊了，所以初始prior的選擇可以基於一些規則或經驗，例如說“按照模型參數所給出的變異數必須是正定的”，不過變異數本來就是二次方的矩 (moment)，不管模型參數多麼稀奇古怪，它必然是正定的；而診斷例子中的$p(天花)=10^{-4}$，只能說good doctor認為不失一般性下染上天花的機率很小而已，本質上這些prior好像說了什麼，細看又什麼都沒說，這些prior就稱作無信息先驗機率 (uninformative prior)。這裡的無信息表明了這些prior盡可能是客觀不帶有偏見地去影響posterior的結果。舉例像是統派會說：“長遠來看兩岸合久必分分久必合，統一是必然的”，基於這個prior我們可以推論出一個posterior那就是兩岸必然統一，但什麼時候會統一，這個prior的資訊是完全沒有畫面的！ (它只是客觀地講出他認為最終的可能，但沒法告訴你是明天、下個月或是明年會統一)² 這裡還是要強調一下，即便是最無信息的prior，它仍是我們計算posterior的一部分 (權重)，所以它仍然會造成一定的影響，而貝氏推論最後也不必然等於古典推論的結果，而最大的那個datum likelihood (maximum likelihood estimate) 也不一定影響posterior最多，它可能因為得到比較小的prior而效應變差 (像$p(天花|痘斑)$雖然是MLE，但$p(天花)$小很多，最後也壓低了病人得的是天花的後驗機率)。

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。

---

² 一個例子像是隨便假設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$的支持只能算勉勉強強，基本上兩個模型不分軒輊

不過通常，我們會用Eq. (6.5)改寫Eq. (6.10)進而得到 $$Q_{21}=\frac{p(D|M_2)p(M_2)}{p(D|M_1)p(M_1)}.\tag{6.11}$$ 這裡的$M$假設是由參數$\boldsymbol{\theta}$所張成的 (spanned)，那麼 $$p(D|M)\equiv E(M)=\int p(D|M,\boldsymbol{\theta})d\boldsymbol{\theta}\tag{6.12}$$ 稱為$M$的全域似然度 (global likelihood)，指涉若$M$是正確的則$D$被觀測到的機率。而比$B_{21}=E(M_2)/E(M_1)$則叫作貝氏因子 (Bayes factor)。

2.1 再一個例子：拋硬幣

數學總是很抽象，我們舉的例子來輔助。我手邊有一枚硬幣，拋了$N$次後要拋出幾次正面才能說硬幣是公正的?

表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拆成兩部分，剩下的留到下篇再討論囉！