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所示，其中实线是基于估计参数绘制的数据分布曲线，虚线是系统自动生成的密度曲线，可见拟合效果还是比较理想的。