极大似然估计与后验概率估计

两种参数估计的统计学方法

---

今天主要复习一下两种参数估计的统计学方法，分别是极大似然估计（MLE）和最大后验概率估计（MAP）。

• 问题背景
• MLE
• MAP
• MLE 与 MAP各自的优缺点
• MLE与MAP之间的区别与联系

问题背景

以掷硬币为例。现在我们一共抛掷了10次硬币，其结果为{H,T,T,T,H,H,H,H,H,H}。我们假设硬币朝上的概率为$P(\theta )$,现在问题来了？如何从我们的观测数据当中近似得到$P(\theta )$的估计值呢？
对于上述问题，我们可以从两个角度出发进行思考，一是从频率学派的角度，二是从贝叶斯学派的角度。

MLE

使用极大似然估计的方法得到$P(\theta )$的近似值主要分为以下两个步骤：

• 先根据数据集对数据分布做出假设（或者构建模型）。以掷硬币为例，一种很自然的想法即为：假设掷硬币朝上的概率$P(H)$满足二项分布。对于参数为$(n,p)$的二项分布，我们知道其分布列为$$p(k)=\left(\genfrac{}{}{0ex}{}{n}{k}\right)p^k(1-p{)}^{n-k}$$具体到掷硬币来说，我们可以得到$$P(D\mid \theta )=\left(\begin{array}{c}{n}_H+n_T\\ n_H\end{array}\right){\theta }^{n_H}(1-\theta {)}^{n_T},$$
  其中$n_H$表示硬币朝上的次数，$n_T$表示硬币朝下的次数，$\theta$表示硬币朝上的概率。
• 确定模型参数，使得我们所观测到的数据尽可能具有代表性（make our data more likely to be observed in real world)，也即:$${\hat{\theta }}_{MLE}={\mathrm{argmax}}_{\theta }P(D;\theta )$$，上述问题的求解过程如下：
  $$\begin{gathered} {\hat{\theta }}_{MLE}={\mathrm{argmax}}_{\theta }P(D;\theta ) \\ ={\mathrm{argmax}}_{\theta }\left(\begin{array}{c}{n}_H+n_T\\ n_H\end{array}\right){\theta }^{n_H}(1-\theta {)}^{n_T} \\ ={\mathrm{argmax}}_{\theta }\log \left(\begin{array}{c}{n}_H+n_T\\ n_H\end{array}\right)+n_H\cdot \log (\theta )+n_T\cdot \log (1-\theta ) \\ ={\mathrm{argmax}}_{\theta }{n}_H\cdot \log (\theta )+n_T\cdot \log (1-\theta ) \end{gathered}$$
  对于上述表达式求导：
  $$\frac{n_H}{\theta }=\frac{n_T}{1-\theta }⟹n_H-n_H\theta =n_T\theta ⟹\theta =\frac{n_H}{n_H+n_T}$$这样一来，我们就得到了关于参数$\theta$的估计。
  仔细一看$\hat{\theta }$的值，我们不难发现，其恰好为我们平时说的频率。也就是说，于掷硬币这个例子而言，极大似然估计最终以频率来估计硬币朝上的概率。

MAP

以掷硬币实验中硬币朝上事件发生的频率来估计硬币朝上的概率，这似乎是一种很自然的做法。下面我们来考虑另外一种极端情况：同样是抛10次硬币，结果为{H,H,H,H,H,H,H,H,H,H}。此时我们使用MLE来估计参数$\theta$的值，就会得到$P(H)\approx \hat{\theta }=1$.这样的结果你觉得可靠吗？至少我觉得是不可靠的。因为凭我的直觉，我认为硬币朝上的概率和硬币朝下的概率是五五开的。于是，贝叶斯学派认为，纯粹使用频率来估计概率是不行的。他们认为参数$\theta$是一个随机变量，在估计硬币朝上的概率$P(H)$的时候，还应当引入先验知识，将$P(H)$也即$\theta$本身服从的概率分布也考虑进来。
所以MAP要解决的是一个什么样的问题呢？简而言之，MAP就是要通过我们现有的观测数据去寻找最有可能的$\theta$，即：$${\hat{\theta }}_{MAP}={argmax}_{\theta }P(\theta \mid D)$$
回顾一下贝叶斯公式：$$P(\theta \mid D)=\frac{P(D\mid \theta )\ast P(\theta )}{P(D)}=\frac{P(D\mid \theta )\ast P(\theta )}{P(D)}$$
在上面表达式当中，各个参数代表的含义如下：

• $P(D)$：对于分母这一项，如果是一个连续性随机变量的话，我们可以写成${\int }_{\theta }P(D\mid \theta )\mathrm{d}\theta$，如果是一个离散型随机变量我们可以写成${\sum }_{\theta }P(D\mid \theta )$。不论是何种情况，$P(D)$的取值都不会影响我们求解最佳的参数$\theta$。
• $P(D\mid \theta )$：这一项恰好就是我们先前在MLE中要最大化的量，也即在给定参数$\theta$的前提下，我们所观测到的数据在真实情况下出现的可能性。
• $P(\theta )$：参数$\theta$的先验概率。
• $P(\theta \mid D)$：基于观测数据，我们对于参数$\theta$所满足的先验概率的矫正。

下面我们就逐一分析分子中的各项：

- $P(\theta )$

通常情况下，对于掷硬币这个例子而言，我们可以假设$P(H)$服从Beta分布，即$$P(\theta )=\frac{{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1}}{B(\alpha ,\beta )}$$
在上面的表达式当中，分母$B(\alpha ,\beta )=\frac{\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha +\beta )}$。
为什么我们使用Beta分布来模拟硬币朝上的概率$P(H)$所满足的分布呢？原因主要有以下两点：

• 其与二项分布是共轭先验的（Conjugate_prior）。所谓共轭先验就是先验分布是beta分布，而后验分布同样是beta分布。$$P(\theta \mid D)\propto P(D\mid \theta )P(\theta )\propto {\theta }^{n_H+\alpha -1}(1-\theta {)}^{n_T+\beta -1}$$
• Beta分布特指一组定义在(0,1) 区间的连续概率分布，而参数$\theta$的取值范围恰好为(0,1）。

- $P(D\mid \theta )$

这一项我们已经在MLE中分析过来，这里就不再加以阐述。

寻找MAP的最优解

$$\begin{gathered} {\hat{\theta }}_{MAP}={\mathrm{argmax}}_{\theta }P(\theta |Data) \\ ={\mathrm{argmax}}_{\theta }\frac{P(Data|\theta )P(\theta )}{P(Data)}\text{(By Bayes rule)} \\ ={\mathrm{argmax}}_{\theta }\log (P(Data|\theta ))+\log (P(\theta )) \\ ={\mathrm{argmax}}_{\theta }{n}_H\cdot \log (\theta )+n_T\cdot \log (1-\theta )+(\alpha -1)\cdot \log (\theta )+(\beta -1)\cdot \log (1-\theta ) \\ ={\mathrm{argmax}}_{\theta }(n_H+\alpha -1)\cdot \log (\theta )+(n_T+\beta -1)\cdot \log (1-\theta ) \\ ⟹{\hat{\theta }}_{MAP}=\frac{n_H+\alpha -1}{n_H+n_T+\beta +\alpha -2} \end{gathered}$$

MLE与MAP的注意点

• MLE

  • 前面有提到MLE需要我们针对问题事先提出一个假设模型，因此模型的正确与否将很大程度上决定我们能否准确估计参数$\theta$的值。
  • MLE在很大程度上受样本个数n的影响。当$n⟶\infty$并且我们的模型假设也合理的时候，MLE能够对参数$\theta$作出较好的估计，相反当$n⟶0$,也即样本个数较少的时候，使用MLE对参数$\theta$进行估计时，容易发生过拟合的情况，比如先前$\hat{\theta }=1$的情况。

• MAP

  • MAP的结果与我们的先验知识有关。如果我们对于参数的先验假设不合理的话，那么使用MAP很可能会得到一个糟糕的结果。这一点我们可以通过接下来的实验进行论述。
  • MAP在$n⟶0$的时候，相较于MLE而言，会取得更好的结果，当然前提是我们的先验假设要正确，否则结果仍然会很糟糕。而当$n⟶\infty$的时候，MLE和MAP都将收敛到参数$\theta$的真实值附近。

为了能够更直观地理解上述结论，我决定自己手动书写代码进行验证：

import numpy as np
import matplotlib.pyplot as plt
def compute_mle_and_map(thets, gama0, gama1,size = 500):
    MLE_es = []
    MAP_es = []
    for i in range(size):
        m = np.random.binomial(i+1, theta)
        prob1 = (m * 1.0)/ (i + 1)
        prob2 = (m + gama1) * 1.0 / (i + gama0 + gama1)
        MLE_es.append(prob1)
        MAP_es.append(prob2)
    MLE_es = np.array(MLE_es)
    MAP_es = np.array(MAP_es)
    return MLE_es,MAP_es

在上面这段代码中，我假定掷硬币朝上的真实概率为0.3，然后通过二项分布来模拟掷硬币的过程，其中prob1是采用MLE得到的结果，prob2是采用MAP得到的结果，而循环变量i主要控制投掷硬币的次数。

size = 250
theta = 0.3
x = np.arange(size) + 1
y1, y2 = compute_mle_and_map(0.3, 42, 18, size)
plt.figure() #创建图形
plt.subplot(2,2,1)
plt.plot(x,y1, color='blue', label= 'MLE')
plt.plot(x,y2, color='red', label='MAP')
plt.hlines(theta, 0, size, colors='green', label='True theta')
plt.legend()
plt.xlabel("Number of coin flips")
plt.ylabel("Low Confident Priors")
plt.title("true theta=0.3,gama0=42,gama2=18")
plt.subplot(2,2,3)
y1, y2 = compute_mle_and_map(0.3, 84, 36, size)
plt.plot(x,y1, color='blue', label= 'MLE')
plt.plot(x,y2, color='red', label='MAP')
plt.hlines(theta, 0, size, colors='green', label='True theta')
plt.legend()
plt.xlabel("Number of coin flips")
plt.ylabel("High Confident Priors")
plt.title("true theta=0.3,gama0=84,gama2=36")
plt.subplot(2,2,2)
y1, y2 = compute_mle_and_map(0.3, 36, 24, size)
plt.plot(x,y1, color='blue', label= 'MLE')
plt.plot(x,y2, color='red', label='MAP')
plt.hlines(theta, 0, size, colors='green', label='True theta')
plt.legend()
plt.xlabel("Number of coin flips")
plt.ylabel("Low Confident Priors")
plt.title("true theta=0.3,gama0=36,gama2=24")
plt.subplot(2,2,4)
y1, y2 = compute_mle_and_map(0.3, 72, 48, size)
plt.plot(x,y1, color='blue', label= 'MLE')
plt.plot(x,y2, color='red', label='MAP')
plt.hlines(theta, 0, size, colors='green', label='True theta')
plt.legend()
plt.xlabel("Number of coin flips")
plt.ylabel("Low Confident Priors")
plt.title("true theta=0.3,gama0=72,gama2=48")
plt.show() 

上面这段代码主要根据硬币的投掷结果做出图像以方便分析：

• 在左上角的第一幅图像当中，我们的先验知识认为硬币朝上的概率是0.3，这与真实情况相符合，因此在数据量比较小的时候，采用MAP方法能够较为准确地估计参数$\theta$的值。此外我们也可以很明显地看到，采用MLE估计参数$\theta$的值，在数据量较小的时候，很容易出现过拟合的情况，表现出较为剧烈的抖动。但是从最终的变化趋势来看，随着实验次数的不断增加，MLE对应的曲线最只能怪还是收敛到了$\theta$的真实值0.3附近。
• 在左下角的图形当中，我特意将$\alpha$和$\beta$的值变为了左上角所对应图像的两倍，这样做的目的相当于增加了我们对于先验知识的置信度，由于我们对于$\theta$的先验符合真实情况，所以左上角和左下角的两幅图像并无本质区别。
  那么一旦我们的先验出现了偏差，结果会如何呢？

下面我们把视角转向右侧的两幅图像。

• 先看右上角的图像，我们设定$\alpha$和$\beta$的值分别为24、36，也即我们对于硬币朝上概率的先验假设为0.4。此时我们发现MAP在实验次数比较少的时候，并没有得到关于参数$\theta$的准确估计，而是随着实验次数的增加，慢慢从0.4收敛到0.3.
• 再看看右下角的图像，我们进一步增大我们对于先验的置信度，也即设定$\alpha =48$、$\beta =72$。对比右侧的上下两个图像，我们可以发现对于右上角的图像而言，实验次数约为100次时,MAP对应的曲线收敛到0.3,；而对于右下角的图像而言，实验次数则在200次左右时才收敛至0.3.
• 这说明了什么？这说明了，使用MAP方法得到的结果一定程度上依赖于我们的先验知识。如果我们的先验假设是正确的，那么无论实验次数是多是少，MAP都会得到一个可靠的结果。相反，若我们的先验本身就是错误的，那么MAP将会得到一个不可靠的结果，且随我们对先验知识的置信程度的增加，我们需要获取更多的观测数据，进行更多次的实验，才能进一步校正我们的先验知识，进而得到对参数$\theta$的可靠估计。

MLE 与MAP的联系

在我来看，MLE是MAP的一种特殊情况。当我们认为$P(\theta )$满足均匀分布的时候，MAP就转化为了MLE。换句话说，在频率学派来看，参数$\theta$在其各个取值上是等概率的，而贝叶斯学派则认为，参数$\theta$的分布是由观测者的先验知识决定的，不能简单认为其就满足均匀分布。

补充

下面的代码我自己对于Tom Mitchell « Machine Learning » ch2 « Estimating Probabilities» 课后的习题解答，解答未必正确，仅供参考。

from scipy import stats
p = np.linspace(0, 1, 1000)
beta0 = 43 
beta1 = 19
pbeta = stats.beta.pdf(p, beta0, beta1)
plt.plot(p, pbeta, label='prior probability distribution')
plt.legend()
plt.show()
pbeta = stats.beta.pdf(p, beta0+6, beta1+9)
plt.grid()
plt.plot(p, pbeta, label='MAP distribution(prior gama0=42,gama1=18)')
plt.show()

实验结果如下图：

beta0 = 420
beta1 = 180
pbeta = stats.beta.pdf(p, beta0+6+1, beta1+9+1)
plt.grid()
plt.plot(p, pbeta, label='MAP distribution(prior gama0=420,gama1=180)')
plt.legend()
plt.show()
beta0 = 32
beta1 = 28
pbeta = stats.beta.pdf(p, beta0+6+1, beta1+9+1)
plt.grid()
plt.plot(p, pbeta, label='MAP distribution(prior gama0=32,gama1=28)')
plt.legend()
plt.show()

实验结果如下图：