從伯努利試驗到最大概似估計


本文主要是接續介紹了中央極限定理和信賴區間的上一篇文(為什麼信心水準不是機率?)的應用。這段推論統計曾是高中選修數學中被用來演練信賴區間的典型模板,只是信賴區間在現行的 108 課綱中已經被徹底刪除(倒是有類似的概念在選修物理不起眼地以不確定度的姿態出現)。不過撰寫本文的初衷其實是為了解決另一個大問題,這個問題和本專欄《見微知著》的首篇文章(為什麼樣本標準差要除以 n-1?)息息相關,那不如在現在就一次搞定,上一篇文末賣的關子在下篇文再繼續。

伯努利試驗與分布

這是一個在統計學開始發展之前就出現的經典機率問題,17 世紀中已有許多數學家(例如法國的費馬和帕斯卡、荷蘭的惠更斯)從賭博遊戲的問題出發開創了機率論,瑞士數學家雅各布・伯努利(Jacob Bernoulli)在 17 世紀末也接著做了許多分析,其中大量討論到結果僅有兩種可能的簡單試驗,例如「成功」或「失敗」、硬幣的「正」或「反」等,我們現在把這類試驗叫做「伯努利試驗(Bernoulli trials)」。

用現在符號表示的話,這相當於隨機變數 \(X\) 只有 \(0\) 或 \(1\) 兩種的統計分布,令試驗成功的機率為 \(p\),期望值與變異數分別為:

\[
E\left(X\right)=1\times p+0\times \left(1-p\right)=p\\
\mathrm{Var}\left(X\right)=E\left(X^2\right)-E^2\left(X\right)=p-p^2=p\left(1-p\right)
\]

它的機率分布被叫做「伯努利分布(Bernoulli distribution)」,記為 \(X\sim\mathrm{Ber}(p)\),其機率質量函數也當然就是:

\[
P\left(X=k\right)=
\begin{cases}
p,\, &k=1\\
1-p,\, &k=0
\end{cases}\\
\]

在前文已經介紹過中央極限定理,我們知道無論抽出的 \(X\) 來自於何種分布的母體,當抽樣的樣本數 \(n\to\infty\),隨機變數 \(\overline{X}\) 形成的「抽樣分布」都會趨近常態分布,即 \(\overline{X}\sim N\left(\mu,\displaystyle{\frac{\sigma^2}{n}}\right)\)。那現在對於伯努利分布的母體當然也適用,「抽 \(n\) 個樣本」在這裡相當於「做 \(n\) 次伯努利試驗」,得到的 \(0\) 或 \(1\) 的平均值 \(\overline{X}\) 可以解讀成這些重複試驗中成功次數(\(1\) 的個數)的比例(proportion),則其分布等價於抽樣分布,期望值與變異數為:

\[
E\left(\overline{X}\right)=\frac{1}{n}\cdot E\left(\sum_{i=1}^{n}X_i\right)=\frac{np}{n}=p\\
\mathrm{Var}\left(\overline{X}\right)=\frac{1}{n^2}\cdot \mathrm{Var}\left(\sum_{i=1}^{n}X_i\right)=\frac{np(1-p)}{n^2}=\frac{p(1-p)}{n}
\]

機率函數也成了連續的機率密度函數,自行代入常態分布即可。所以我們可以說 \(p\) 相當於 \(1\) 的比例,在這裡一般會把樣本比例 \(\overline{X}\) 改用 \(\hat{p}\) 來表示對母體比例 \(p\) 的估計,並再次提醒第二條的成立前提是 \(X_i\) 彼此都獨立,這在之前已經提到過。於是根據中央極限定理,當試驗次數 \(n\to\infty\) 時,\(\hat{p}\sim N\left(p,\displaystyle{\frac{p(1-p)}{n}}\right)\)。

課本在介紹以上內容時往往會用引入「二項分布(binomial distribution)」來推導,這只是把重複 \(n\) 次的伯努利試驗也當成一個機率分布來描述期望值與變異數,令隨機變數 \(Y=\displaystyle{\sum_{i=1}^{n}}X_i\sim B\left(n,p\right)\),則:

\[
E\left(Y\right)=E\left(\sum_{i=1}^{n}X_i\right)=np\\
\mathrm{Var}(Y)=\mathrm{Var}\left(\sum_{i=1}^{n}X_i\right)=np\left(1-p\right)
\]

其機率質量函數各位也理應相當熟悉:

\[
P\left(Y=k\right)=C^n_k p^k(1-p)^{n-k},\, k=0,1,2,\dots,n
\]

但我為了和前面的論述整齊對應而刻意迴避了這個術語,補在這裡提供讀者比較。話說其實這番推論是逆著歷史脈絡走的,最初的中央極限定理本就出自於深受伯努利影響的另一位法國數學家棣美弗(Abraham de Moivre),依此得到了所謂的 de Moivre–Laplace 定理,在現在看來就只是中央極限定理的伯努利試驗版本。

同場加映:幾何分布(geometric distribution)

高中數學在 108 課綱雖然刪除了信賴區間,但新增了隨機變數 \(X\) 為進行重複伯努利試驗直到成功一次所需要的試驗次數(含成功的那一次)的機率分布,機率質量函數為:
\[
P\left(X=k\right)=\left(1-p\right)^{k-1} p,\, k=1,2,3,\dots
\]

函數值會是公比為 \(1-p\) 的等比數列,所以被叫做「幾何」分布,\(X\sim \mathrm{Geo}(p)\)。至於課本提供的期望值與變異數的證明可能會有點繁瑣,因為呈等比的機率會要再乘上對應的隨機變數,隨機變數的取值本身又是公差為 \(1\) 的等差數列,而形成一個無窮項的差比級數(然後通分再彼此交錯相減)。這裡我就不再呈現一樣的東西了,改用其他方式來說明供讀者參考。

期望值可以用非常直觀的方式理解。以擲硬幣的情景為例,若正面的機率為 \(\displaystyle{\frac{1}{5}}\),表示我們可以預期長期下來平均每丟 \(5\) 次會有 \(1\) 次是正面,換句話說就是丟出 \(1\) 次正面平均需要 \(5\) 次,於是期望值就是這試驗成功機率的倒數:
\[
E\left(X\right)=\frac{1}{p}
\]

變異數我想直接用 \(\mathrm{Var}\left(X\right)=E\left(X^2\right)-E^2\left(X\right)\) 來求,只要求出前面的 \(E\left(X^2\right)\) 就能立刻得證。首先把情況分割成兩部分:第一次試驗就成功(機率為 \(p\))和第一次試驗之後才成功(機率為 \(1-p\));隨機變數 \(X\) 在前者取值即為 \(1\),在後者則令它為 \(1+X’\)。可得:
\[
\begin{align*}
E\left(X^2\right)&=p\times 1+\left(1-p\right)E\left(\left(1+X’\right)^2\right)\\
&=p+\left(1-p\right)\left(1+2E\left(X’\right)+E\left(X’^2\right)\right)
\end{align*}
\]
而你應可以很直觀地理解,在第一次失敗的前提下,接下來到成功需要的次數,跟打從一開始需要的次數,機率分布是完全相同的,這其實就是所謂的「無記憶性(memoryless property)」,且它相當於告訴我們 \(E\left(X’\right)=E\left(X\right)\) 和 \(E\left(X’^2\right)=E\left(X^2\right)\),於是:
\[
\begin{align*}
&E\left(X^2\right)=p+\left(1-p\right)\left(1+2E\left(X\right)+E\left(X^2\right)\right)\\
&\Rightarrow E\left(X^2\right)=\frac{2-p}{p^2}\\
&\Rightarrow \mathrm{Var}\left(X\right)=\frac{2-p}{p^2}-\frac{1}{p^2}=\frac{1-p}{p^2}
\end{align*}
\]

模型導向或設計導向

當初我在讀完以上高中課程內容之後,腦中總有一種難以名狀的感覺揮之不去,如果要先大概描述我所困惑的點的話,或許可以說是:究竟以上的「分布」和現實的連結是什麼?

對於典型的「抽樣(sampling)」的情境,我們可以想像母體是一群真實存在的數值,而抽樣則是在這些數值中選取,其隨機性來自於我們「設計」的抽樣行為本身,但我們不可能真正得知母體的所有數據,所以在合理的前提下以我們發展的理論「模型」來分析,其隨機性則來自於母體本身的隨機變數的生成。可是在伯努利試驗的情境中,我們雖然可以刻意想像母體好比存在無數個由 \(0\) 或 \(1\) 組成的數列,每一次試驗就是從中抽取 \(1\) 個數,但如果這個試驗的情境是擲硬幣卻會變得非常奇怪,我們好像不能把它理解為「抽樣」?試想:這個抽樣是從哪裡「抽」出來的呢?

不過如果把情境換成也可以想像母體有確切存在的對象被抽取的話,比方說公投,可以想像每個人都有同意與不同意的選擇,那在公投之前抽幾位民眾來民調的結果就能符合「抽樣」的想像了。

後來我才知道,這是統計理論很根本的哲學問題,究竟我們要用「模型導向(model-based)」還是「設計導向(design-based)」的思路來理解母體?目前我所討論的主要依循的是前者,先建立統計分布的模型從它來產出隨機數據,分析其性質,再套用在現實中來說明觀察到的結果,如同以上我們從理論(中央極限定理)來說明抽樣分布會如何近似到常態,又如同物理學和數學的各種建模皆是如此;然而現實中很多情景是基於後者,想像的模型並非實體,先是母體的數據本就真實存在(只是仍然未知),再來才藉由我們設計的隨機試驗來判斷母體可能潛在的模式,這在生物醫學領域非常廣泛,很多時候我們其實不清楚根本的機制,但卻可以從經評估夠可靠的統計結果來論事——繞過演繹並訴諸歸納。

那回顧前面的分析你就會發現,擲硬幣或公投的情境都符合伯努利試驗的描述,但在現實層面上,前者似乎更符合「模型導向」的哲學,因為可以想像硬幣的表面形狀、質量、投擲的角度、力道、環境等因素理論上如果被完整計算的話,我們理應可以從每一次的隨機參數來產出各種可預知但不確定的結果(隨機試驗);後者則更符合「設計導向」的哲學,畢竟大家對公投內容的同意與否就算不做民調也本就存在,不是有了試驗才被隨機生成。而有趣的正是,這兩件事的結果是等價的,所以即便我們的理論主要出自「模型導向」的方式,它也可以被用在「設計導向」的情境上。

於是即便隨機性的來源不同,兩種哲學思路最後都同樣來到建立在「抽樣分布」的基礎上的標準誤差(standard error)和信賴區間(confidence interval)等問題,雖然「抽樣」這個詞好像偏重設計導向,但建立它的「分布」卻又是很模型導向的行為,在我看來,這著實是很關鍵的思路匯聚的體現!我想也因此 \(\overline{X}\) 這個很「模型導向」的描述才會同時又常用 \(\hat{p}\) 這種很具有「設計導向」意義的方式來表示,畢竟殊途同歸。

更進一步地,我們還可以思考的是究竟因果關係(causality)是出自於哪一種哲學?或者說,哪種思路才可以合理用來斷言結果所揭示的因果關係?我並沒有打算在這裡選邊站,姑且只是想稍微點出這種曖昧的來源而已,起碼在認清這個問題點之後,我所困惑的源頭已經被清楚抓出來了。

回到本文脈絡,「設計導向」的思路其實一直以比較隱晦的方式潛伏在這些統計學討論中,如前面所說,他的隨機性得來自於抽樣本身的設計,但這個課題屬於比較進階的統計學課程,如系統抽樣(systematic sampling)或群集抽樣(cluster sampling)等,所以我們一般頂多只在介紹基本的統計量時(如標準差和期望值)提供「抽樣調查」的情景來做初步的想像;到開始介紹抽樣分布之後卻不著痕跡的切換成「模型導向」的思路,此時隨機性已經來自模型的隨機試驗,於是現在要處理的問題變成是,我們獲得的數據究竟如何對應這個模型?

或是更學術一點地說:如何去進行機率分布的擬合(fitting)?

最大概似估計(MLE)

這是最被廣泛使用的機率分布的擬合方法,不過在這之前我得先說明「概似(likelihood)」的概念。

假設有一個封閉的水域裡面有不知多少的 \(N\) 條魚,我隨機撈了 \(10\) 隻出來標記後放回去,過一陣子再回來隨機撈了 \(20\) 隻,發現其中有 \(4\) 隻是有被我標記的魚,於是我們可以很直覺的估計:

\[
\frac{10}{\hat{N}}=\frac{4}{20} \Rightarrow \hat{N}=50
\]

讀過生物的讀者肯定都聽過這就是所謂的「捉放法(capture/recapture method)」,但不知各位是否曾注意到這個推論非常跳,因為即便每隻魚被取到的機率相同讓取到被標記的魚的機率真的是 \(\displaystyle{\frac{10}{N}}\),憑什麼我們就能認為抽出的樣本比例也會滿足這個機率(母體的比例)?仔細想想,樣本並沒有必要契合這件事啊。

事實上,我們就是無法知道母體 \(N\) 到底是多少,也無法回顧當初抽樣時產生現在這個樣本的機率究竟真正是多少,不過我們可以分析對於不同的 \(N\),產生這個樣本的機率分別是多少(即便我們不知哪一個是事實),這些機率被特別稱為「概似」(也有人翻譯為「似然」),而所謂的「最大概似估計(maximum likelihood estimate, MLE)」就是指我們打算選用機率最大的那種情況來作為估計!換言之,母體對應的參數是多少,才最有可能(有最大的概似)會產生這種樣本?接下來將繼續以這情況為例來簡單呈現為何在前例中的 \(N\) 如果是 \(\hat{N}\) 的話,最有可能產生這樣的樣本,於是 \(\hat{N}\) 是 \(N\) 的 MLE。

對於未知的參數 \(N\),他的值可以是大於等於 \(20\)(上述情景起碼有的魚的數量)的任意整數,其中有 \(10\) 隻被標記的魚,而我們現在抽出 \(20\) 隻中有 \(4\) 隻也有標記的機率(概似)是:

\[
\frac{C^{10}_{4}\times C^{N-10}_{20-4}}{C^N_{20}}
\]

它的(概似)函數圖形會長這樣:

最大值發生在N為50的概似函數圖形

可以看到,在 \(N=50\) 的時候會有最大值(其實在 \(N=49\) 時也剛好是),確認了 \(\hat{N}=50\) 的確是 MLE。不過如果要推導出更一般性的通式會稍微有點複雜,往往需要些計算量:假設當初標記的總數量是 \(T\),樣本數 \(n\) 中有 \(t\) 個是被標記的,上述的機率(概似)和「超幾何分布(hypergeometric distribution)」有一樣的數學形式(前後項會是一個固定的有理函數),只是在超幾何分布會用抓到的魚中有被標記的數量 \(t\) 當作隨機變數,母體參數 \(N\) 會是固定的,也因此這樣才會是同一個機率空間中的互斥事件,於是會是個機率分布;然而現在計算的概似使用的母體是變動的,上圖的各個機率總和也就不必然會是 \(1\),所以不能說他是機率質量函數而是「概似函數」請特別注意,統計學在這裡也常很謹慎地一律改用「概似」並迴避「機率」這個可能讓你誤以為是機率分布的術語。

那麼對於不同可能的母體數量 \(N\),概似 \(L\) 的前後項比例會是:

\[
\frac{L_{N}}{L_{N-1}}=\frac{\displaystyle{\frac{C^{T}_{t}\times C^{N-T}_{n-t}}{C^N_{n}}}}{\displaystyle{\frac{C^{T}_{t}\times C^{N-1-T}_{n-t}}{C^{N-1}_{n}}}}=\frac{\left(N-T\right)\left(N-n\right)}{N\left(N-T-n+t\right)}
\]

如果它 \(>1\) 的話表示機率嚴格遞增,可以找到以下範圍:

\[
\begin{align*}
&\left(N-T\right)\left(N-n\right)>N\left(N-T-n+t\right)\\
&\Rightarrow N^2-nN-TN+nT>N^2-TN-nN+tN\\
&\Rightarrow N<\frac{nT}{t}
\end{align*}
\]

所以會有最大概似的 \(N\) 會是不超過 \(\displaystyle{\frac{nT}{t}}\) 的最大整數,在這裡 \(\displaystyle{\frac{20\times 10}{4}=50}\) 剛好是整數,讓 \(49\) 和 \(50\) 都在概似函數的最高處,但當然是後者才貼合最初我們的直覺啦。MLE 的概念是 20 世紀初的英國演化生物學家 R. A. Fisher 提出的,本就稍微有接觸統計學的肯定都聽過這位幾乎可以說是開創整個統計學的天才,可是在他之前就已有捉放法這種估算了(後來還稱之為「Lincoln–Petersen 估計」),只是 Fisher 把這種直覺用很明確的邏輯表達出來,精準描述其數學意義。

希望以上稍微繁複的推演沒有模糊了焦點,我要說明的只是最大概似的這個選擇、這種估計法。同樣的邏輯可以套用在其他機率分布,但在何時有最大值的證明也會更加複雜,對於連續且多參數的機率密度函數還時常要取對數來方便後續對特定參數的偏微,這裡我打算只用前例來交代簡單的定性描述,直接呈現其他結論,證明就留給有興趣的讀者去翻數理統計課本了:

  1. 伯努利分布的期望值 \(p\) 的 MLE 是樣本的 \(\hat{p}\),變異數 \(p\left(1-p\right)\) 的 MLE 是 \(\hat{p}\left(1-\hat{p}\right)\)
  2. 常態分布的期望值 \(\mu\) 的 MLE 是樣本的 \(\overline{X}\),變異數 \(\sigma^2\) 的 MLE 是 \(\displaystyle{\frac{\displaystyle{\sum_{i=1}^{n}}(X_i-\overline{X})^2}{n}}\)

簡單來說,就是直接把樣本當母體算就對了,而且要特別注意常態分布的變異數(或標準差)不是經過 Bessel 校正的 \(s^2\) 喔,為加以區分可以用 \(\hat{\sigma}^2\) 表示,於是你肯定會有疑問,這樣豈不是就沒有「不偏(unbiased)」了?沒錯,而且這本就和 MLE 是兩碼子事,關於這點請繼續看下去。

再談「不偏性」

了解到 MLE 作為估計量的選用的判斷基準後,最後再來回顧一下之前證明「樣本標準差」要除以 \(n-1\) 所用到的「不偏性(unbiasedness)」究竟是什麼樣的概念。

在許多統計教材中會在介紹點估計時提及估計量的眾多性質,來綜合描述該估計量是否是一個「好的」估計,不偏性就是其一,除此之外還有「一致性(consistency)」、「有效性(efficiency)」等等。這些和 MLE 並不是同一個層級的問題,而是在那之後所做的評估和調整,於是例如我們可以針對小樣本的抽樣來說因為 \(\hat{\sigma}\) 不滿足不偏性,最好再用「Bessel 校正」來微調成不偏的 \(s\) 來做點估計,甚至區間估計的標準誤也比照辦理;不過因果關係更正確的說法應該是從自由度的角度來詮釋的那部分,這直接影響後續衍生出的對小樣本非常實用的「t-分布」,不偏性則只是如此設計的漂亮結果之一(關於 t-分布的介紹下篇文再繼續)。

如果樣本夠大的話,到底要不要把 \(n\) 改用 \(n-1\) 其實也相對不重要了,減 \(1\) 對數據的影響已經相當不顯著(t-分布也趨近常態分布),MLE 的結果也隨著樣本數越大而越不偏,這被稱為「漸進不偏性(asymptotic unbiasedness)」,所以課本在介紹上述以伯努利試驗出發的統計推論時往往預設樣本取得夠大的前提(我記得以前高中課本的說法是超過 \(30\) 個),並告訴你直接將 \(\hat{p}\) 當 \(p\) 來用即可,依此來計算信賴區間。

當然,如果你還是想微調成「不偏」也可以,\(E\left(\hat{p}\right)=p\) 在前面呈現過了(注意 \(\hat{p}\) 和 \(\overline{X}\) 是同一件事),接著要進一步確認的是把 \(\hat{p}\) 直接拿來套入伯努利分布的母體變異數 \(\mathrm{Var}\left(X\right)\) 是否不偏,正如同之前我們想直接把 \(\overline{X}\) 直接套入變異數(或標準差)公式來求母體的 \(\sigma^2\):

\[
\begin{align*}
E\left(\hat{p}\left(1-\hat{p}\right)\right)&=E\left(\hat{p}-\hat{p}^2\right)\\
&=p-\left(p^2+\frac{p(1-p)}{n}\right)\\
&=\frac{p}{n}\cdot \left(n-np-1+p\right)\\
&=\frac{p}{n}\cdot \left(n-1\right)\left(1-p\right)\\
&=\frac{n-1}{n}\cdot p(1-p)
\end{align*}
\]

很明顯地,如果讓樣本變異數改使用 \(\displaystyle{\frac{n}{n-1}}\cdot\hat{p}(1-\hat{p})\) 來計算,它的期望值就會是 \(p(1-p)=\mathrm{Var}\left(X\right)\),和我們熟悉的樣本變異數公式一樣,都是乘上 \(\displaystyle{\frac{n}{n-1}}\) 這個因子來修正成「不偏的」估計量。然而課本通常不會在這裡討論到這地步,我個人覺得可能的理由是伯努利試驗本身的樣本很單純(如硬幣正或反、公投同意或不同意皆一拍兩瞪眼),於是實務上獲取大樣本的成本本就不高,不過以上的推演對熟悉數理統計的人來說應也是小菜一碟就是了。

況且即便這樣做的意義又在哪呢?還記得在之前的文章有稍微提過,不偏性畢竟也是描述估計量的期望值,它確保的「不偏」是長期下來的結果,但我們的統計量的計算只是那「一組」\(n\) 個樣本(或 \(n\) 次伯努利試驗),他究竟是否是恰好體現母體的那組樣本(或伯努利試驗),我們始終無法得知,只是在用如此數學性質說明我們對他抱有何種期待,所以其實大可不必死守不偏性,把除以 \(n-1\) 的結果當成是「讓小樣本的誤差範圍更大來更保守估計」來理解反而是更務實的理由,雖然單只是這樣的話就無法解釋為何減 \(1\) 了。

「說了這麼多,結果只是個冠冕堂皇的設計?」或許有人會這樣想吧,統計學好像繞了一大圈又沒什麼進展一樣,但我對此的理解是:「是這樣,但不是這樣。」

這些設計本就不是在讓各位可以發掘出更多新資訊,如果你以為統計學可以透過某種神奇的運算來從樣本中挖出母體的真實情況,那我覺得你誤會大了。雖然說是在用樣本來「估計」母體,但這些本質上只是在「描述」,用定義明確的方式來說明我們以什麼準則出發(MLE)、用什麼方式整理數據(統計量)、又是用什麼策略來下結論(統計推論),否則每個人都可以用自己慣用或偏好的邏輯來呈現各自的調查結果,那即便初始數據本身客觀,後續的計算會隱含了許多主觀因素,最後淪為無法相互比較的雞同鴨講。

帶給我們新資訊的從頭到尾都是「統計」的這個行為本身,不是統計學的設計;這些設計只是一套明確的 SOP,透過經統計學家討論出的最不失公允的脈絡讓你得以用沒有歧義的方式與其他人溝通,同時也為嘗試再現你所做的研究的任何人提供一致的分析流程。毫無疑問地,統計學著實屬於科學體系中相當重要的一環,接下來本專欄的文章預計也會開始朝更多自然科學領域所需的統計推論技術邁進。

靖之
靖之

一個記性不太好卻又什麼都想學的大學生,相信著透過不斷地思考、整理、撰文,終能將無涯學海,塞入兩耳間不到二十公分寬的頭顱中。

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *