极大似然估计的定义

相关定义

极大似然估计（MLE）是建立在极大似然原理的基础上的一种参数估计的方法。所谓极大似然原理，通俗理解就是，一个事件发生的原因就是该事件发生的概率较大。假设在一个随机试验中有若干结果事件$A、B、\dots$，若$A$出现，则可以认为该试验对$A$的出现有利，即该试验中$A$出现的概率较大。例如，设甲箱中有99个白球，1个黑球；乙箱中有1个白球。99个黑球。现随机取出一箱，再从抽取的一箱中随机取出一球，结果是黑球，则这一黑球从乙箱抽取的概率比从甲箱抽取的概率大得多，这时我们自然更多地相信这个黑球是取自乙箱的。

一般来说，事件$A$发生的概率与未知参数$\theta$有关，其发生的概率为$P(A|\theta)$。若在一次试验中$A$发生了，MLE则要求我们应该选择使$P(A|\theta)$达到最大的参数$\theta$。在详细介绍MLE之前，先让我们来了解一下何为参数化模型。参数化模型就是用有限个参数来表示概率密度函数集$\mathrm{q}(\boldsymbol{x} ; \boldsymbol{\theta})$。其中$\boldsymbol{x}$为随机变量，$\boldsymbol{\theta}=\left(\theta^{(1)}, \cdots, \theta^{(b)}\right)^{\mathrm{T}}$表示参数向量，设维度为$b$。根据MLE，我们==确定参数$\boldsymbol{\theta}$时应该尽可能地最大化训练样本${\boldsymbol{x}}^n_{i=1}$的概率==。因此，我们可以将在参数$\boldsymbol{\theta}$下的训练样本${\boldsymbol{x}}^n_{i=1}$的概率定义为似然函数，记为$L(\boldsymbol{\theta})$，并假设$L(\boldsymbol{\theta})$独立同分布：
$$
L(\theta) = \prod^n_{i=1}\mathrm{q}(\boldsymbol{x_i} ; \boldsymbol{\theta})
$$
则MLE的优化问题为最大化该似然函数：
$$
\hat{\boldsymbol{\theta}}{\mathrm{ML}}=\underset{\boldsymbol{\theta} \in \Theta}{\arg \max } L(\boldsymbol{\theta}),\\
密度估计：\hat{p}(\boldsymbol{x})=q\left(\boldsymbol{x} ; \hat{\boldsymbol{\theta}}{\mathrm{ML}}\right)
$$
若$\mathrm{q}(\boldsymbol{x} ; \boldsymbol{\theta})$关于$\boldsymbol{\theta}$可微，则$\hat{\boldsymbol{\theta}}{\mathrm{ML}}$应满足如下：
$$
\left.\frac{\partial}{\partial \boldsymbol{\theta}} L(\boldsymbol{\theta})\right|{\boldsymbol{\theta}=\hat{\boldsymbol{\theta}}{\mathrm{ML}}}=\mathbf{0}{b}
$$
此即为似然方程。应注意的是，该方程的解不一定就是极大似然的解，而极大似然解通常满足该方程。原因如图：

由于$L(\theta)$是连乘的，且对数函数为单调递增的，所以我们通常将其通过对数函数转换为连加的形式，这种形式称为$\text{log}-似然函数$：
$$
\hat{\boldsymbol{\theta}}{\mathrm{ML}}=\underset{\boldsymbol{\theta} \in \Theta}{\arg \max } \log L(\boldsymbol{\theta})=\underset{\boldsymbol{\theta} \in \Theta}{\arg \max }\left[\sum{i=1}^{n} \log q\left(\boldsymbol{x}{i} ; \boldsymbol{\theta}\right)\right]
$$
同样求偏微分令为零可以得到如下似然方程：
$$
\left.\frac{\partial}{\partial \boldsymbol{\theta}} \log L(\boldsymbol{\theta})\right|{\boldsymbol{\theta}=\boldsymbol{\theta}{\mathrm{M} .}}=\mathbf{0}{b}
$$

高斯模型

这里的高斯模型为服从高斯分布的$d$维参数化模型：
$$
q(\boldsymbol{x} ; \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{\frac{d}{2}} \operatorname{det}(\mathbf{\Sigma})^{\frac{1}{2}}} \exp \left(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right)
$$

式中的向量$\boldsymbol{\mu}$和矩阵$\boldsymbol{\Sigma}$都是高斯模型的参数，分别对应期望和方差-协方差矩阵：
$$
\begin{aligned}
&\boldsymbol{\mu}=E[\boldsymbol{x}]=\int x q(\boldsymbol{x} ; \boldsymbol{\mu}, \mathbf{\Sigma}) \mathrm{d} x \
&\boldsymbol{\Sigma}=V[\boldsymbol{x}]=\int(\boldsymbol{x}-\boldsymbol{\mu})(\boldsymbol{x}-\boldsymbol{\mu})^{\mathrm{T}} q(\boldsymbol{x} ; \boldsymbol{\mu}, \mathbf{\Sigma}) \mathrm{d} x
\end{aligned}
$$
`　　对训练样本${\boldsymbol{x}i}^n{i=1}$而言，其对应的高斯模型为$q(\boldsymbol{x} ; \boldsymbol{\mu}, \boldsymbol{\Sigma})$，$\text{log}-$似然估计为
$$
\log L(\boldsymbol{\mu}, \mathbf{\Sigma})=-\frac{n d \log 2 \pi}{2}-\frac{n \log (\operatorname{det}(\boldsymbol{\Sigma}))}{2}-\frac{1}{2} \sum_{i=1}^{n}\left(\boldsymbol{x}{i}-\boldsymbol{\mu}\right)^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}\left(\boldsymbol{x}{i}-\boldsymbol{\mu}\right)
$$
则高斯模型的似然方程为：
$$
\left{\begin{array}{l}
\left.\frac{\partial}{\partial \mu} \log L(\boldsymbol{\mu}, \mathbf{\Sigma})\right|{\mu=\hat{\mu}{\mathrm{ML}}}=\mathbf{0}{d} \
\left.\frac{\partial}{\partial \mathbf{\Sigma}} \log L(\boldsymbol{\mu}, \mathbf{\Sigma})\right|{\Sigma=\hat{\Sigma}{\mathrm{ML}}}=\boldsymbol{O}{d \times d}
\end{array}\right.
$$
其中：
$$
\begin{aligned}
&\frac{\partial \log L}{\partial \boldsymbol{\mu}}=-n \boldsymbol{\Sigma}^{-1} \boldsymbol{\mu}+\boldsymbol{\Sigma}^{-1} \sum_{i=1}^{n} \boldsymbol{x}{i} \
&\frac{\partial \log L}{\partial \boldsymbol{\Sigma}}=-\frac{n}{2} \mathbf{\Sigma}^{-1}+\frac{1}{2} \mathbf{\Sigma}^{-1}\left(\sum{i=1}^{n}\left(\boldsymbol{x}{i}-\boldsymbol{\mu}\right)\left(\boldsymbol{x}{i}-\boldsymbol{\mu}\right)^{\mathrm{T}}\right) \mathbf{\Sigma}^{-1}
\end{aligned}
$$

其中向量和矩阵的偏导公式如下：
$$
\begin{gathered}
\frac{\partial \boldsymbol{\mu}^{\text{T}} \mathbf{\Sigma}^{-1} \mu}{\partial \boldsymbol{\mu}}=2 \boldsymbol{\Sigma}^{-1} \mu, \frac{\partial \boldsymbol{x}^{\text{T}} \boldsymbol{\Sigma}^{-1} \mu}{\partial \boldsymbol{\mu}}=\boldsymbol{\Sigma}^{-1} \boldsymbol{x} \
\frac{\partial \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1} \boldsymbol{x}}{\partial \boldsymbol{\Sigma}}=-\boldsymbol{\Sigma}^{-1} \boldsymbol{x} \boldsymbol{x}^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}, \frac{\partial \log (\operatorname{det}(\boldsymbol{\Sigma}))}{\partial \boldsymbol{\Sigma}}=\boldsymbol{\Sigma}^{-1}, \frac{\partial \operatorname{tr}\left(\tilde{\boldsymbol{\Sigma}}^{-1} \boldsymbol{\Sigma}\right)}{\partial \boldsymbol{\Sigma}}=\overline{\mathbf{\Sigma}}^{-1}
\end{gathered}
$$

求解以上似然方程可以得到如下极大似然估计值：
$$
\begin{aligned}
&\hat{\boldsymbol{\mu}}{\mathrm{ML}}=\frac{1}{n} \sum{i=1}^{n} x_{i} \
&\hat{\mathbf{\Sigma}}{\mathrm{ML}}=\frac{1}{n} \sum{i=1}^{n}\left(\boldsymbol{x}{i}-\hat{\boldsymbol{\mu}}{\mathrm{ML}}\right)\left(\boldsymbol{x}{i}-\hat{\boldsymbol{\mu}}{\mathrm{ML}}\right)^{\mathrm{T}}
\end{aligned}
$$
此时我们可以发现$\hat{\boldsymbol{\mu}}{\mathrm{ML}}$和$\hat{\mathbf{\Sigma}}{\mathrm{ML}}$分别对应着样本均值和样本协方差。而当我们有足够的训练样本时，$\hat{\mathbf{\Sigma}}_{\mathrm{ML}}$是可逆的。

一维高斯模型MLE的Python代码实现

def d1_gauss_mle(n, m, s):
    """
    n: 样本数量
    m: 参数mu
    s: 参数sigma
    """
    x = s*np.random.randn(n, 1)+m  # 生成数据
    mh = np.mean(x)  # 计算mu的估计值(样本均值)
    sh = np.std(x, axis=0)  # 计算sigma的估计值(样本标准差)
    X = np.linspace(-4, 4, 100)
    Y = np.exp(-np.dot(np.mat((X-m)).T,np.mat((X-m)))/(2*s**2))/(2*np.pi*s)  # 真实值
    Yh = np.exp(-np.dot(np.mat((X-mh)).T,np.mat((X-mh)))/(2*sh**2))/(2*np.pi*sh)  # 估计值
    return x, X, Y, Yh

sample, x, y, y_hat = d1_gauss_mle(100, 0, 1)

plt.figure(figsize=(8,4), dpi=100)
plt.plot(x, np.diag(y), 'b-',label='真实值')
plt.plot(x, np.diag(y_hat), 'y--', label='估计值')
plt.plot(sample, np.zeros_like(sample), 'o')
plt.legend()
plt.xlim([-4,4])
plt.ylim([0,0.2])
plt.show()

mle

## 高斯混合模型

现实生活中，若一种模式分布在不同的聚类中，用单一的高斯模型去近似其条件分布是不合适的。如图所示：

当我们用单峰模型来拟合多峰分布时，即使当样本数量趋于无穷，模型的最终表现依旧很差。在此，定义高斯混合模型如下：
$$
q(\boldsymbol{x} ; \boldsymbol{\theta})=\sum_{\ell=1}^{m} w_{\ell} N\left(\boldsymbol{x} ; \boldsymbol{\mu}{\ell}, \mathbf{\Sigma}{\ell}\right)
$$
　　高斯混合模型适用于近似如上的多峰模型分布。我们可以看到，高斯混合模型其实就是多个单峰模型的加权线性组合，其中：
$$
N(\boldsymbol{x} ; \boldsymbol{\mu}, \mathbf{\Sigma})=\frac{1}{(2 \pi)^{d / 2} \operatorname{det}(\boldsymbol{\Sigma})^{1 / 2}} \exp \left(-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{\mu})^{\mathrm{T}} \boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu})\right)
$$
高斯混合模型的参数$\boldsymbol{\theta}$为：
$$
\boldsymbol{\theta}=\left(w_{1}, \cdots, w_{m}, \boldsymbol{\mu}{1}, \cdots, \mu{m}, \boldsymbol{\Sigma}{1}, \cdots, \boldsymbol{\Sigma}{m}\right)
$$
而$q(\boldsymbol{x} ; \boldsymbol{\theta})$要想成为一个概率密度函数应满足：
$$
\forall x \in X, q(x ; \boldsymbol{\theta}) \geqslant 0 \text { 且 } \int_{X} q(x ; \boldsymbol{\theta}) \mathrm{d} x=1
$$
且${\omega_\ell}^m_{\ell=1}$应该满足：
$$
w_{1}, \cdots, w_{m} \geqslant 0 \text { 且 } \sum_{\ell=1}^{m} w_{\ell}=1
$$
通过python生成一个带有$3$个高通分量的高斯混合模型：

n = 30
n_sample = 1000
w = np.array([0.5, 0.3, 0.2])
m = np.array([-4, -1, 3])
s = np.array([1, 1.5, 2])
x_min = -10
x_max = 10
x = np.zeros((len(m),n))
for i in range(len(m)):
    x[i, :] = (s[i]*np.random.randn(n, 1)+m[i]).flatten()
X = np.zeros((len(m),n_sample))
for k in range(len(m)):
    X[k, :] = (np.linspace(x_min, x_max, n_sample))
Y = np.zeros_like(X)
for j in range(len(m)):
    Y[j, :] = np.diag(np.exp(-np.dot(np.mat((X[j]-m[j])).T,np.mat((X[j]-m[j])))/(2*s[j]**2))/(2*np.pi*s[j]))
g = np.dot(w, Y)
plt.figure(figsize=(8, 4), dpi=200)
plt.xlim([-10, 10])
plt.ylim([0,0.2])
for num in range(len(m)):
    plt.plot(X[num, :], Y[num, :], '-.', label='高斯分量g'+str(num+1))
    plt.plot(x[num, :], np.ones_like(x[num, :])*0.025*(num+1), 'o')  
plt.title("g1:m=%.1f, s=%.1f; g2:m=%.1f, s=%.1f; g3:m=%.1f, s=%.1f"%(m[0],s[0],m[1],s[1],m[2],s[2]))
plt.plot(X[num, :], g, 'black', label='高斯混合模型')
plt.legend()
plt.show()

该模型为：
$$
q(x)=0.5 N\left(x ;-4,1^{2}\right)+0.3 N\left(x ; -1,1.5^{2}\right)+0.2 N\left(x ; 3,2^{2}\right)
$$

高斯混合模型的极大似然估计

高斯混合模型中的参数$\boldsymbol{\theta}$可参考上节。模型的似然函数为：
$$
L(\boldsymbol{\theta})=\prod_{i=1}^{n} q\left(\boldsymbol{x}{i} ; \boldsymbol{\theta}\right)=\prod{i=1}^{n} \sum_{\ell=1}^{m} w_{\ell} N\left(\boldsymbol{x}i ; \boldsymbol{\mu}{\ell}, \mathbf{\Sigma}{\ell}\right)
$$
通过MLE来找到使该似然函数最大化的参数。该优化问题为约束优化问题：
$$
\hat{\boldsymbol{\theta}}=\underset{\boldsymbol{\theta}}{\operatorname{argmax} L}(\boldsymbol{\theta})\\
\text{s.t.}\ \ \ \ w{1}, \cdots, w_{m} \geqslant 0 \text { 且 } \sum^m_{\ell=1} w_{\ell}=1
$$
由于问题带约束，因此不能通过简单求导来求解。将其中$w_{1}, \cdots, w_{m}$参数化为：
$$
w_{\ell}=\frac{\exp \left(\gamma_{\ell}\right)}{\sum_{\ell^{\prime}=1}^{m} \exp \left(\gamma_{\ell’}\right)}
$$
式中$\left{\gamma_{\ell}\right}{\ell=1}^{m}$满足约束条件$\gamma{\ell}\in(0,1]$，可以通过算法来学习到。将上式代入到似然函数中去并转化为对数似然函数：
$$
\log L(\boldsymbol{\theta})=\sum_{i=1}^{n} \log \sum_{\ell=1}^{m} \exp \left(\gamma_{\ell}\right) N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell} ; \mathbf{\Sigma}{\ell}\right)-n \log \sum{\ell=1}^{m} \exp \left(\gamma_{\ell}\right)
$$

将上式分别对参数$\gamma_{\ell}、\boldsymbol{\mu}\ell$和$\boldsymbol{\Sigma}\ell$求偏导，得
$$
\begin{aligned}
\frac{\partial}{\partial \gamma_{\ell}} \log L(\boldsymbol{\theta})=& \sum_{i=1}^{n} \frac{\exp \left(\gamma_{\ell}\right) N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell}, \boldsymbol{\Sigma}{\ell}\right)}{\sum{\ell^{\prime}=1}^{m} \exp \left(\gamma_{\ell’}\right) N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell’}, \boldsymbol{\Sigma}{\ell’}\right)}-\frac{n \gamma{\ell}}{\sum_{\ell^{\prime}=1}^{m} \exp \left(\gamma_{\ell^{\prime}}\right)} \
&=\sum_{i=1}^{n} \frac{w_{\ell} N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell}, \mathbf{\Sigma}{\ell}\right)}{\sum{\ell^{\prime}=1}^{m} w_{\ell’} N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell’}, \boldsymbol{\Sigma}{\ell’}\right)}-n w{\ell}\
&=\sum_{i=1}^{n} \eta_{i, \ell}-m w_{\ell}
\end{aligned}
$$
$$
\begin{aligned}
&\frac{\partial}{\partial \boldsymbol{\mu}{\ell}} \log L(\boldsymbol{\theta})=\sum{i=1}^{n} \eta_{i, \ell} \mathbf{\Sigma}{\ell}^{-1}\left(\boldsymbol{x}{i}-\boldsymbol{\mu}{\ell}\right) \
&\frac{\partial}{\partial \boldsymbol{\Sigma}{\ell}} \log L(\boldsymbol{\theta})=\frac{1}{2} \sum_{i=1}^{n} \eta_{i, \ell}\left(\boldsymbol{\Sigma}{\ell}^{-1}\left(\boldsymbol{x}{i}-\boldsymbol{\mu}{\ell}\right)\left(\boldsymbol{x}{i}-\boldsymbol{\mu}{\ell}\right)^{\top} \boldsymbol{\Sigma}{\ell}^{-1}-\boldsymbol{\Sigma}_{\ell}^{-1}\right)
\end{aligned}
$$

其中定义$\eta_{i, \ell}$如下：
$$
\eta_{i, \ell}=\frac{w_{\ell} N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell}, \mathbf{\Sigma}{\ell}\right)}{\sum{\ell^{\prime}=1}^{m} w_{\ell’} N\left(\boldsymbol{x}{i} ; \boldsymbol{\mu}{\ell’}, \boldsymbol{\Sigma}{\ell’}\right)}
$$
令偏导式为零，即可得到如下参数估计：
$$
\left{\begin{array}{l}
\hat{w}{\ell}=\frac{1}{n} \sum_{i=1}^{n} \hat{\eta}{i, \ell} \
\hat{\boldsymbol{\mu}}{\ell}=\frac{\sum_{i=1}^{n} \hat{\eta}{i, \ell} \boldsymbol{x}{i}}{\sum_{i^{\prime}=1}^{n} \hat{\eta}{i^{i}, \ell}} \
\hat{\boldsymbol{\Sigma}}{\ell}=\frac{\sum_{i=1}^{n} \hat{\eta}{i, \ell}\left(\boldsymbol{x}{i}-\hat{\boldsymbol{\mu}}{\ell}\right)\left(\boldsymbol{x}{i}-\hat{\boldsymbol{\mu}}{\ell}\right)^{\mathrm{T}}}{\sum{i^{\prime}=1}^{n} \hat{\eta}_{i \cdot \ell}}
\end{array}\right.
$$

其中参数$\hat{\eta}{i^{\prime}, \ell}$是样本$x_i$的第$\ell$个组成部分的决定因素。与其定义类似：
$$
\hat{\eta}{i^{\prime} \cdot \ell}=\frac{\hat{w}{\ell} N\left(x{i} ; \hat{\mu}{\ell}, \hat{\mathbf{\Sigma}}{\ell}\right)}{\sum_{\ell^{\prime}=1}^{m} \hat{w}{\ell} N\left(\boldsymbol{x}{i} ; \hat{\boldsymbol{\mu}}{\ell}, \hat{\mathbf{\Sigma}}{\ell}\right)}
$$
　　事实上，关于以上参数目前还没有较好的方法来求解。但我们可以通过数值计算的方法来迭代求解。

高斯混合模型的梯度上升算法

这种算法其实就是梯度下降。通过让目标函数值升高来迭代更新参数。当约束条件较好时能够保证算法收敛至局部最优（local maxima）。算法通过计算梯度来选择参数更新方向。但该算法存在超参数。我们需要初始化其算法更新步长（学习率）。而步长的设定不太容易，我们往往通过设置不同的初始值放入算法进行迭代来选择最优值或通过模拟退火来逐步降低学习率来解决。

算法的Python实现