2.3　極大似然估計

正如本章最初所講的，統計推斷的基本問題可以分為兩大類：一類是參數估計；另一類是假設檢驗。其中假設檢驗又分為參數假設檢驗和非參數假設檢驗兩大類。本章所講的假設檢驗都屬于是參數假設檢驗的范疇。參數估計也分為兩大類，即參數的點估計和區間估計。用于點估計的方法一般有矩方法和最大似然估計法（Maximum Likelihood Estimate，MLE）兩種。

2.3.1　極大似然法的基本原理

最大似然這個思想最初是由德國著名數學家卡爾·高斯（Carl Gauss）提出的，但真正將其發揚光大的則是英國的統計學家羅納德·費希爾（Ronald Fisher）。費希爾在其1922年發表的一篇論文中再次提出了最大似然估計這個思想，并且首先探討了這種方法的一些性質。而且，費希爾當年正是憑借這一方法徹底撼動了皮爾遜在統計學界的統治地位。從此開始，統計學研究正式進入了費希爾時代。

為了引入最大似然估計法的思想，來看一個例子。假設一個口袋中有黑白兩種顏色的小球，并且知道這兩種球的數量比為3:1，但不知道具體哪種球占3/4，哪種球占1/4。現在從袋子中有返回地任取三個球，其中有一個是黑球，那么試問袋子中哪種球占3/4，哪種球占1/4。

設X是抽取三個球中黑球的個數，又設p是袋子中黑球所占的比例，則有X~B（3，p），即

當X=1時，不同的p值對應的概率分別為

由于第一個概率小于第二個概率，所以我們判斷黑球的占比應該是1/4。

在上面的例子中，p是分布中的參數，它只能取3/4或者1/4。需要通過抽樣結果來決定分布中參數究竟是多少。在給定了樣本觀察值以后再去計算該樣本的出現概率，而這一概率依賴于p的值。所以就需要用p的可能取值分別去計算最終的概率，在相對比較之下，最終所取之p值應該是使得最終概率最大的那個p值。

極大似然估計的基本思想就是根據上述想法引伸出來的。設總體含有待估參數θ，它可以取很多值，所以就要在θ的一切可能取值之中選出一個使樣本觀測值出現的概率為最大的θ值，記為，并將此作為θ的估計，并稱為θ的極大似然估計。

首先來考慮X屬于離散型概率分布的情況。假設在X的分布中含有未知參數θ，記為

P（X=a_i）=p（a_i；θ），　i=1，2，…，θ∈Θ

現從總體中抽取容量為n的樣本，其觀測值為x₁，x₂，…，x_n，這里每個x_i為a₁，a₂，…中的某個值，該樣本的聯合分布為

由于這一概率依賴于未知參數θ，故可將它看成是θ的函數，并稱其為似然函數，記為

對不同的θ，同一組樣本觀察值x₁，x₂，…，x_n出現的概率也不一樣。當P（A）＞P（B）時，事件A出現的可能性比事件B出現的可能性大，如果樣本觀察值x₁，x₂，…，x_n出現了，當然就要求對應的似然函數的值達到最大，所以應該選取這樣的作為θ的估計，使得

如果存在的話，則稱為θ的極大似然估計。

此外，當X是連續分布時，其概率密度函數為p（x；θ），θ為未知參數，且θ∈Θ，這里的Θ表示一個參數空間。現從該總體中獲得容量為n的樣本觀測值x₁，x₂，…，x_n，那么在X₁=x₁，X₂=x₂，…，X_n=x_n時聯合密度函數值為

它也是θ的函數，也稱為似然函數，記為

對不同的θ，同一組樣本觀察值x₁，x₂，…，x_n的聯合密度函數值也是不同的，因此應該選擇θ的極大似然估計，從而使下式得到滿足

2.3.2　求極大似然估計的方法

當函數關于參數可導時，可以通過求導方法來獲得似然函數極大值對應的參數值。在求極大似然估計時，為求導方便，常對似然函數取對數，稱為對數似然函數，它與在同一點上達到最大。根據微積分中的費馬定理，當l（θ）對θ的每一分量可微時，可通過l（θ）對θ的每一分量求偏導并令其為0求得，稱

為似然方程，其中k是θ的維數。

下面就結合一個例子來演示這個過程。假設隨機變量X~B（n，p），又知x₁，x₂，…，x_n是來自X的一組樣本觀察值，現在求P（X=T）時，參數p的極大似然估計。首先寫出似然函數

然后對上式左右兩邊取對數，可得

將l（p）對p求導，并令其導數等于0，得似然方程

解似然方程得

可以驗證，當時，?²l（p）/?p²＜0，這就表明可以使函數取得極大值。最后將題目中已知的條件代入，可得p的極大似然估計為。

再來看一個連續分布的例子。假設有隨機變量X~N（μ，σ²），μ和σ²都是未知參數，x₁，x₂，…，x_n是來自X的一組樣本觀察值，試求μ和σ²的極大似然估計值。首先寫出似然函數

然后對上式左右兩邊取對數，可得

將l（μ，σ²）分別對μ和σ²求偏導數，并令它們的導數等于0，于是可得似然方程

求解似然方程可得

而且還可以驗證和可以使得l（μ，σ²）達到最大。用樣本觀察值替代后便得出μ和σ²的極大似然估計分別為

因為是μ的無偏估計，但并不是σ²的無偏估計，可見參數的極大似然估計并不能確保無偏性。

最后給出一個被稱為“不變原則”的定理：設是θ的極大似然估計，g（θ）是θ的連續函數，則g（θ）的極大似然估計為。

這里并不打算對該定理進行詳細證明。下面將通過一個例子來說明它的應用。假設隨機變量X服從參數為λ的指數分布，x₁，x₂，…，x_n是來自X的一組樣本觀察值，試求λ和E（X）的極大似然估計值。首先寫出似然函數

然后對上式左右兩邊取對數，可得

將l（λ）對λ求導得似然方程為

解似然方程得

可以驗證它使l（λ）達到最大，而且上述過程對一切樣本觀察值都成立，所以λ的極大似然估計值為。此外，E（X）=1/λ，它是λ的函數，其極大似然估計可用不變原則進行求解，即用代入E（X），可得E（X）的最大似然估計為，這與矩法估計的結果一致。

2.3.3　極大似然估計應用舉例

上一小節演示了通過解方程?l（θ）/?θ_j=0從而求得參數θ的極大似然估計值的基本方法。但顯而易見的是，這個求解過程非常復雜，本節將通過幾個實例來演示在R中進行極大似然估計的方法。

對于不同的分布形式而言，其似然函數的形式也是各式各樣的，所以最后得到的似然方程解（也即是參數的極大似然估計值）的表達式也很難統一。很難找到一種通用的方法來對所有情況下的參數做極大似然估計。因此使用R語言進行極大似然估計，往往是先要確定似然函數的表達式，然后再借助于R中的極值求解函數來完成。

在單參數情況下，可以使用R中的函數optimize（）求極大似然估計值，它的調用格式如下。

函數optimize（）的作用是在由參數interval指定的區間內搜索函數f的極值。這個區間也可以由參數lower（即區間的下界）和upper（即區間的上界）來控制。默認情況下，參數maximum=FALSE表示求極小值，如果將其置為TRUE則表示求極大值。

例如，現在已知某批電子元件的使用壽命服從參數為λ的指數分布，λ未知且有λ＞0。現在隨機抽取一組樣本并測得其使用壽命如下（單位：小時）

518　612　713　388　434

請嘗試用極大似然估計其這批產品的平均壽命。

上一節的最后已經求出了指數函數的對數似然函數形式，可以用R語言代碼將似然函數如下。

然后用optimize（）求使得似然函數取得極值時的參數λ的估計值，結果如下。

由此便求出了參數λ的估計值為0.001 878 689，再根據上一小節最后得出的結論，可知這批電子元件的平均使用壽命E（X）=1/λ。

而且這個結果與之前推導的結論，當X服從指數分布時，E（X）的最大似然估計為，并由此算得的結果是一致的。

再來看一個稍微復雜的例子，這次要估計的參數將有多個。首先在R中導入程序包MASS中的數據geyser，示例代碼如下。

該數據集是地質學家記錄的美國黃石公園內忠實泉（Old Faithful），如圖2-6所示，一年內的噴發數據，數據有兩個變量，分別是泉水持續涌出的時間（eruptions）和噴發相隔的時間（waiting），在這個例子中我們將僅會用到后者。

現在我們打算對變量waiting的分布進行擬合，于是首先通過直方圖來大致了解一下數據的分布形態。執行下面的代碼，其結果如圖2-7所示。

從繪制的結果來看，圖形中有兩個峰，很像是兩個分布疊加在一起而成的結果，于是可以推斷分布是兩個正態分布的混合，故用下面的函數來描述

p（x）=αN（x；μ₁，σ₁）+（1-α）N（x；μ₂，σ₂）

所以在構建的模型中，需要估計的參數有5個，即α、μ₁、σ₁、μ₂和σ₂。上述分布函數的對數極大似然函數為

接下來，在R中定義對數似然函數，示例代碼如下。由于在后面將要使用的極值求解法會在迭代過程中產生一些似然函數不能處理的無效值，盡管這并不會影響到最終的求解結果，但是為了避免出現不必要的警告信息，此處使用了suppressWarnings（）函數來忽略那些警告信息。

為了進行極大似然估計，下面將調用R語言中的程序包maxLik，該包為進行極大似然估計提供了諸多便利，在多參數估計時可以考慮使用它。在進行極值求解時，可以通過修改maxLik（）函數中的參數method來選擇不同的數值求解方法。可選的值有“NR”“BHHH”“BFGS”“NM”和“SANN”5種，默認情況下函數將使用默認值“NR”，即采用Newton-Raphson算法。

最后通過圖形來評估一下采用最大似然法所估計出來的參數擬合效果，示例代碼如下。

執行以上代碼，結果如圖2-8所示，其中實線是基于估計參數繪制的數據分布曲線，虛線是系統自動生成的密度曲線，可見擬合效果還是比較理想的。