最大似然估計

在統計學中，最大似然估計（maximum likelihood estimation，簡作MLE），也稱極大似然估計，是用來估計一個機率模型的參數的一種方法。

## 預備知識

下方的討論要求讀者熟悉機率論中的基本定義，如機率分布、機率密度函數、隨機變量、數學期望等。讀者還須先熟悉連續實函數的基本性質，比如使用微分來求一個函數的極值（即極大值或極小值）。

同時，讀者須先擁有似然函數的背景知識，以了解最大似然估計的出發點及應用目的。

## 最大似然估計的原理

給定一個機率分布$D$，已知其機率密度函數（連續分布）或機率質量函數（離散分佈）為$f_D$，以及一個分佈參數$\theta$，我們可以從這個分布中抽出一個具有$n$個值的採樣$X_1, X_2,\ldots, X_n$，利用$f_D$計算出其似然函數：

$$\mbox{L}(\theta\mid x_1,\dots,x_n) = f_\theta(x_1,\dots,x_n).$$

若$D$是離散分布，$f_\theta$即是在參數為$\theta$時觀測到這一採樣的機率；若其是連續分布，$f_\theta$則為$X_1, X_2,\ldots, X_n$聯合分布的機率密度函數在觀測值處的取值。一旦我們獲得$X_1, X_2,\ldots, X_n$，我們就能求得一個關於$\theta$的估計。最大似然估計會尋找關於$\theta$的最可能的值（即，在所有可能的$\theta$取值中，尋找一個值使這個採樣的「可能性」最大化）。從數學上來說，我們可以在$\theta$的所有可能取值中尋找一個值使得似然函數取到最大值。這個使可能性最大的$\widehat{\theta}$值即稱為$\theta$的最大似然估計。由定義，最大似然估計是樣本的函數。

### 注意

- 這裡的似然函數是指$x_1,x_2,\ldots,x_n$不變時，關於$\theta$的一個函數。
- 最大似然估計不一定存在，也不一定唯一。

### 推導

最大似然估計可以從相對熵推導而來。相對熵衡量了使用一個給定分布$Q$來近似另一個分布$P$時的資訊損失，對於離散型隨機變量，可以用以下公式：

$$D_{\text{KL}}(P||Q)=\sum_i P(i) \log \frac{P(i)}{Q(i)}$$

其中，$P$是真實分布，$Q$是近似分布。在最大似然估計的情景下，假設分布擁有一系列參數$\theta$，我們希望通過樣本得到參數的估計值$\hat{\theta}$。我們可以利用相對熵來評判估計的好壞：

$$D_{\text{KL}}(p_{\theta} (x)||p_{\hat{\theta}} (x))=\sum_{x \in E} p_\theta(x) \log \frac{p_\theta(x)}{p_{\hat{\theta}}(x)}$$

根據期望的定義，我們可以將上式改寫為：

$$\begin{align} D_{\text{KL}}(p_{\theta} (x)||p_{\hat{\theta}} (x))&=\mathbb{E}_{\theta} \left[\log\left(\frac{p_{\theta}(x)}{p_{\hat{\theta}}(x)} \right)\right]\\ &=\mathbb{E}_{\theta}[\log p_\theta (x)] - \mathbb{E}_{\theta} [\log p_{\hat{\theta}}(x)] \end{align}$$

KL值越大，參數估計越壞，因此，需要通過改變估計參數$\hat{\theta}$的值來獲得最小的值，所對應的參數極為最佳估計參數。即：

$$\hat{\theta}_\text{best}=\arg\min_{\hat{\theta}} D_{\text{KL}}(p_{\theta} (x)||p_{\hat{\theta}} (x))$$

假設有$n$個樣本，根據大數定理，可以進行替換：

$$\mathbb{E}_{\theta} [\log p_\hat{\theta}(x)]\rightsquigarrow\frac{1}{n}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)$$

即，可以通過下式評估：

$$D_{\text{KL}}(p_\theta (x)||p_{\hat{\theta}}(x))=\mathbb{E}_{\theta}[\log p_{\theta}(x)]-\frac{1}{n}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)$$

對於一個已知的分布，其參數$\theta$是確定的。因此，$\mathbb{E}_{\theta}[\log p_\theta (x)]$為常數。因此，我們可以通過最小化KL值獲得最佳估計參數：

$$\begin{align} \hat{\theta} &=\arg\min_{\hat{\theta}}\mathbb{E}_{\theta}[\log p_{\theta}(X)]-\frac{1}{n}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)\\ &\Rightarrow\arg\min_{\hat{\theta}}-\frac{1}{n}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)\\ &\Rightarrow\arg\max_{\hat{\theta}}\frac{1}{n}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)\\ &\Rightarrow\arg\max_{\hat{\theta}}\sum_{i=1}^n\log p_{\hat{\theta}}(x_i)\\ &\Rightarrow\arg\max_{\hat{\theta}}\log\left[\prod_{i=1}^n p_{\hat{\theta}}(x_i) \right]\\ &\Rightarrow\arg\max_{\hat{\theta}}\prod_{i=1}^n p_{\hat{\theta}}(x_i) \\ \end{align}$$

因此，要得到最佳參數估計值，只需要最大化$\prod_{i=1}^n p_{\hat{\theta}}(x_i)$，這就是最大似然函數。對於連續型隨機變量，有相同的結論。

## 例子

### 離散分布，離散有限參數空間

考慮一個拋硬幣的例子。假設這個硬幣正面跟反面輕重不同。我們把這個硬幣拋80次（即，我們獲取一個採樣$x_1=\mbox{H}, x_2=\mbox{T}, \ldots, x_{80}=\mbox{T}$並把正面的次數記下來，正面記為H，反面記為T）。並把拋出一個正面的機率記為$p$，拋出一個反面的機率記為$1-p$（因此，這裡的$p$即相當於上方的$\theta$）。假設我們拋出了49個正面，31個反面，即49次H，31次T。假設這個硬幣是我們從一個裝了三個硬幣的盒子裡頭取出的。這三個硬幣拋出正面的機率分別為$p=1/3$, $p=1/2$, $p=2/3$，這些硬幣沒有標記，所以我們無法知道哪個是哪個。使用最大似然估計，基於二項分布中的機率質量函數公式，通過這些試驗數據（即採樣數據），我們可以計算出哪個硬幣的可能性最大。這個似然函數取以下三個值中的一個：

$$\begin{matrix} \mathbb{L}(p=1/3 \mid \mbox{H=49, T=31 }) & = & \mathbb{P}(\mbox{H=49, T=31 }\mid p=1/3) & = & {80\choose 49}(1/3)^{49}(1-1/3)^{31} \approx 0.000 \\ &&\\ \mathbb{L}(p=1/2 \mid \mbox{H=49, T=31 }) & = & \mathbb{P}(\mbox{H=49, T=31 }\mid p=1/2) & = & {80\choose 49}(1/2)^{49}(1-1/2)^{31} \approx 0.012 \\ &&\\ \mathbb{L}(p=2/3 \mid \mbox{H=49, T=31 }) & = & \mathbb{P}(\mbox{H=49, T=31 }\mid p=2/3) & = & {80\choose 49}(2/3)^{49}(1-2/3)^{31} \approx 0.054 \\ \end{matrix}$$

我們可以看到當$\widehat{p}=2/3$時，似然函數取得最大值。

顯然地，這硬幣的公平性和那種拋出後正面的機率是2/3的硬幣是最接近的。這就是$p$的最大似然估計。

### 離散分布，連續參數空間

現在假設例子1中的盒子中有無數個硬幣，對於$0\leq p \leq 1$中的任何一個$p$，都有一個拋出正面機率為$p$的硬幣對應，我們來求其似然函數的最大值：

$$\begin{matrix} \mbox{L}(\theta) & = & f_D(\mbox{H=49,T=80-49}\mid p) = {80\choose 49} p^{49}(1-p)^{31} \\ \end{matrix}$$

其中$0\leq p \leq 1$.

我們可以使用微分法來求極值。方程兩邊同時對$p$取微分，並使其為零。

$$\begin{matrix} 0 & = & {80\choose 49}\frac{d}{dp} \left( p^{49}(1-p)^{31} \right) \\ <pre> & & \\ & \propto & 49p^{48}(1-p)^{31} - 31p^{49}(1-p)^{30} \\ & & \\ & = & p^{48}(1-p)^{30}\left[ 49(1-p) - 31p \right] \\ </pre> \end{matrix}$$

<image>

其解為$p=0$, $p=1$，以及$p=49/80$.使可能性最大的解顯然是$p=49/80$（因為$p=0$和$p=1$這兩個解會使可能性為零）。因此我們說最大似然估計值為$\widehat{p}=49/80$.

這個結果很容易一般化。只需要用一個字母$t$代替49用以表達伯努利試驗中的被觀察數據（即樣本）的「成功」次數，用另一個字母$n$代表伯努利試驗的次數即可。使用完全同樣的方法即可以得到最大似然估計值:

$$\widehat{p}=\frac{t}{n}$$

對於任何成功次數為$t$，試驗總數為$n$的伯努利試驗。

### 連續分布，連續參數空間

最常見的連續機率分布是正態分布，其機率密度函數如下：

$$f(x\mid \mu,\sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(x-\mu)^2}{2\sigma^2}}$$

現在有$n$個正態隨機變量的採樣點，要求的是一個這樣的正態分布，這些採樣點分布到這個正態分布可能性最大（也就是機率密度積最大，每個點更靠近中心點），其$n$個正態隨機變量的採樣的對應密度函數（假設其獨立並服從同一分布）為：

$$f(x_1,\ldots,x_n \mid \mu,\sigma^2) = \left( \frac{1}{2\pi\sigma^2} \right)^\frac{n}{2} e^{-\frac{ \sum_{i=1}^{n}(x_i-\mu)^2}{2\sigma^2}}$$

也可以寫為：

$f(x_1,\ldots,x_n \mid \mu,\sigma^2) = \left( \frac{1}{2\pi\sigma^2} \right)^{n/2} \exp\left(-\frac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{2\sigma^2}\right)$,

這個分布有兩個參數：$\mu,\sigma^2$.有人可能會擔心兩個參數與上方的討論的例子不同，上方的例子都只是在一個參數上對可能性進行最大化。實際上，在兩個參數上的求最大值的方法也差不多：只需要分別把可能性$\mbox{L}(\mu,\sigma) = f(x_1,,\ldots,x_n \mid \mu, \sigma^2)$在兩個參數上最大化即可。當然這比一個參數麻煩一些，但是一點也不複雜。使用上方例子同樣的符號，我們有$\theta=(\mu,\sigma^2)$.

最大化一個似然函數同最大化它的自然對數是等價的。因為自然對數log是一個連續且在似然函數的值域內嚴格遞增的上凹函數。注意：可能性函數（似然函數）的自然對數跟資訊熵以及費雪訊息聯繫緊密。求對數通常能夠一定程度上簡化運算，比如在這個例子中可以看到：

$$\begin{align} 0 & = \cfrac{\partial}{\partial \mu} \log \left( \left( \cfrac{1}{2\pi\sigma^2} \right)^\frac{n}{2} e^{-\frac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{2\sigma^2}} \right) \\ <pre> & = \cfrac{\partial}{\partial \mu} \left( \log\left( \cfrac{1}{2\pi\sigma^2} \right)^\frac{n}{2} - \cfrac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{2\sigma^2} \right) \\ & = 0 - \cfrac{-2n(\bar{x}-\mu)}{2\sigma^2} </pre> \end{align}$$

這個方程的解是$\widehat{\mu} = \bar{x} = \sum^{n}_{i=1}x_i/n$.這的確是這個函數的最大值，因為它是$\mu$裡頭惟一的一階導數等於零的點並且二階導數嚴格小於零。

同理，我們對$\sigma$求導，並使其為零。

$$\begin{align} 0 & = \cfrac{\partial}{\partial \sigma} \log \left( \left( \cfrac{1}{2\pi\sigma^2} \right)^\frac{n}{2} e^{-\frac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{2\sigma^2}} \right) \\ <pre> & = \cfrac{\partial}{\partial \sigma} \left( \cfrac{n}{2}\log\left( \cfrac{1}{2\pi\sigma^2} \right) - \cfrac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{2\sigma^2} \right) \\ & = -\cfrac{n}{\sigma} + \cfrac{ \sum_{i=1}^{n}(x_i-\bar{x})^2+n(\bar{x}-\mu)^2}{\sigma^3} </pre> \end{align}$$

這個方程的解是$\widehat{\sigma}^2 = \sum_{i=1}^n(x_i-\widehat{\mu})^2/n$.

因此，其關於$\theta=(\mu,\sigma^2)$的最大似然估計為：

$\widehat{\theta}=(\widehat{\mu},\widehat{\sigma}^2) = (\bar{x},\sum_{i=1}^n(x_i-\bar{x})^2/n)$.

## 性質

### 泛函不變性（Functional invariance）

如果$\hat{\theta}$是$\theta$的一個最大似然估計，那麼$\alpha = g(\theta)$的最大似然估計是$\hat{\alpha} = g(\hat{\theta})$。函數g無需是一個雙射。

### 漸近線行為

最大似然估計函數在採樣樣本總數趨於無窮的時候達到最小方差，其證明可見於克拉馬－羅下限。當最大似然估計非偏時，等價的，在極限的情況下我們可以稱其有最小的均方差。

對於獨立的觀察來說，最大似然估計函數經常趨於正態分布。

### 偏差

最大似然估計的偏差是非常重要的。考慮這樣一個例子，標有$1$到$n$的$n$張票放在一個盒子中。從盒子中隨機抽取票。如果$n$是未知的話，那麼$n$的最大似然估計值就是抽出的票上標有的$n$，儘管其期望值的只有$(n+1)/2$.為了估計出最高的$n$值，我們能確定的只能是$n$值不小於抽出來的票上的值。

## 歷史

最大似然估計最早是由羅納德·費雪在1912年至1922年間推薦、分析並大範圍推廣的。（雖然以前高斯、拉普拉斯、托瓦爾·尼古拉·蒂勒和F. Y. 埃奇沃思也使用過）。許多作者都提供了最大似然估計發展的回顧。

大部分的最大似然估計理論都在貝葉斯統計中第一次得到發展，並被後來的作者簡化。

## 參見

- 均方差是衡量一個估計函數的好壞的一個量。
- 關於拉奧-布萊克韋爾定理（Rao-Blackwell theorem）的文章中討論到如何利用Rao-Blackwellisation過程尋找最佳不偏估計（即使均方差最小）的方法。而最大似然估計通常是一個好的起點。
- 讀者可能會對最大似然估計（如果存在）總是一個關於參數的充分統計量（sufficient statistic）的函數感興趣。
- 最大似然估計跟廣義動差估計（generalized method of moments）有關。

## 外部連結

- 關於最大似然估計的歷史的一篇論文，作者John Aldrich