引入-最大似然估计

当MVUE难以找到时，MLE 提供了一种寻找实用估计量的替代方法。寻找MLE 有 “按部就班” 的特点，即步骤相对固定。但缺点是，某些情况下推导得到的方程会非常复杂，不易实现。

对于很多实际应用场景，当数据量很大时（10的指数级别），MLE 的性能是最优的（渐近最优）。

标量参数下的最大似然估计

最大似然估计的原理

对于标量参数(\theta)，观测数据为(\mathbf{x})时，有关于参数(\theta)的概率密度函数(p(\mathbf{x};\theta))。MLE 的定义是：找到使似然函数(p(\mathbf{x};\theta))最大的(\theta)值，即：

$\hat{\theta}_{MLE}=\arg\{\max_{\theta} (p(\mathbf{x};\theta))\}$

举个例子来理解，假设对于WGN下振幅为A的直流电平，他的PDF将服从均值$\mu=A$的高斯分布。PDF的极值点就在$A$附近。若估计量为$\hat \theta$，对其直流电平进行最大似然估计，那么能让似然函数$p(\mathbf{x};\theta)$达到极大值的$\hat \theta$就是其均值A

最大似然估计的原理，其实就是PDF的极值点代表观测值在这附近发生的概率最大，去求这个发生概率最大的点。

在真实的操作中，只要观测量$x[n]$足够多，就可以使用离散的点将PDF画出来，此时去寻找极值点即可。

由于对数函数(\ln(\cdot))是单调递增的，最大化似然函数(p(\mathbf{x};\theta))等价于最大化对数似然函数(\ln p(\mathbf{x};\theta))，因此也可写为：

$\hat{\theta}_{MLE}=\arg\max_{\theta} \ln p(\mathbf{x};\theta)$

这样改写之后会对后续处理方便很多。

示例：求解最大似然估计

现有直流信号的观测量(x[n]=A + w[n])（(n = 0,1,\dots,N - 1)）。其直流电平A未知（(A>0)），w[n])是WGN，且噪声方差也为A，未知。估计信号的DC电平

此时概率密度函数（PDF）为：

$p(\mathbf{x},A)=\frac{1}{(2\pi A)^{N/2}}\exp\left(-\frac{1}{2A}\sum_{n = 0}^{N - 1}(x[n]-A)^2\right)$ $\ln p(\mathbf{x};A)=-\frac{N}{2}\ln(2\pi)-\frac{N}{2}\ln A-\frac{1}{2A}\sum_{n=0}^{N-1}(x[n]-A)^2$

步骤1：分析CRLB

先试着用CRLB的方式求出MVUE：

$\frac{\partial \ln p(\mathbf{x},A)}{\partial A}=\frac{1}{2A^2}\sum_{n = 0}^{N - 1}(x[n]-A)^2+\frac{1}{A}\sum_{n = 0}^{N - 1}(x[n]-A)-\frac{N}{2A}$

发现无法将其整理成 “(I(A)(\hat{A}-A))” 的形式，因此不存在高效无偏估计。但是仍可以通过二阶偏导拿到Fisher信息，找到CRLB

$\frac{\partial^2 \ln p}{\partial A^2}= -\frac{1}{A^3}\sum_{n=0}^{N-1}(x[n]-A)^2-\frac{1}{A^2}\sum_{n=0}^{N-1}(x[n]-A)$ $-\frac{1}{A^2}\sum_{n=0}^{N-1}(x[n]-A)-\frac{1}{A^2}\sum_{n=0}^{N-1}1+\frac{N}{2A^2}$ $=\frac{N}{2A^2}-\frac{1}{A^3}\sum(x[n]-A)^2-\frac{2}{A^2}\sum(x[n]-A)-\frac{N}{A}$

Fisher信息为：

$-E\left\{\frac{\partial^2 \ln p}{\partial A^2}\right\}=-[\frac{N}{2A^2}-\frac{1}{A^3}\cdot NA+\frac{2}{A^2}\cdot 0-\frac{N}{A}]$ $=-[\frac{N}{2A^2}-\frac{N}{A^2}-\frac{N}{A}]=\frac{N}{2A^2}+\frac{N}{A}$

CRLB为：

$\text{var}(\hat{A})\geq\frac{1}{N\left(\frac{1}{2A^2}+\frac{1}{A}\right)}=\frac{2A^2}{N(2A + 1)}$

步骤2：计算MLE

MLE是$\ln p$的极值点，因此需要使其一阶偏导数=0，解出这个估计量就可以。

通过最大化对数似然函数（令偏导为 0），得到方程：

$\frac{\partial \ln p(\mathbf{x},A)}{\partial A}=\frac{1}{2A^2}\sum_{n = 0}^{N - 1}(x[n]-A)^2+\frac{1}{A}\sum_{n = 0}^{N - 1}(x[n]-A)-\frac{N}{2A}=0$

消除分母，化简方程：

$\sum_{n = 0}^{N - 1}(x[n]-A)^2+2A\sum_{n = 0}^{N - 1}(x[n]-A)-AN=0$

利用平方展开公式((a-b)^2 = a^2 - 2ab + b^2)，展开求和项：

$\sum(x[n]-A)^2 = \sum x^2[n] - 2A\sum x[n] + NA^2$ $2A\sum(x[n]-A) = 2A\sum x[n] - 2A\sum A = 2A\sum x[n] - 2NA^2$

代入方程并化简：

$\left[\sum x^2[n] - 2A\sum x[n] + NA^2\right] + \left[2A\sum x[n] - 2NA^2\right] - NA = 0$ $\sum x^2[n] - NA^2 - NA = 0$ $A^2 + A - \frac{1}{N}\sum_{n=0}^{N-1}x^2[n] = 0$

一元二次方程(ax^2 + bx + c = 0)求根公式为：

$x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}$

使用求根公式解这个方程：

$\hat{A}_{MLE} = \frac{-1 + \sqrt{1 + 4 \cdot \frac{1}{N}\sum_{n=0}^{N-1}x^2[n]}}{2} = -\frac{1}{2} + \sqrt{\frac{1}{N}\sum_{n=0}^{N-1}x^2[n] + \frac{1}{4}}$

MLE的大样本最优与无偏性

对上例计算期望(E\{\hat{A}_{MLE}\})，发现(E\{\hat{A}_{MLE}\}\neq A)，因此这个 MLE 是有偏估计。

但是其在大样本下时（当观测数N很大时）是无偏的，证明如下：

令(u=\frac{1}{N}\sum_{n = 0}^{N - 1}x^2[n])，则 MLE 可表示为(\hat{A}_{MLE}=g(u)=-\frac{1}{2}+\sqrt{u+\frac{1}{4}})

计算u的期望(u_0 = E\{u\}=A + A^2)，并对(g(u))在(u_0)处做泰勒展开（一阶近似）：

$g(u)\approx g(u_0)+\left.\frac{dg}{du}\right|_{u = u_0}(u - u_0)$

代入(g(u))的表达式并化简后，可得：

$\hat{A}_{MLE}\approx A+\frac{1}{2A + 1}\left\{\frac{1}{N}\sum_{n = 0}^{N - 1}x^2[n]-(A + A^2)\right\}$

分析渐近期望与方差：

渐近期望：(E\{\hat{A}_{MLE}\}\approx A)，说明大样本下 MLE 是无偏的（渐近无偏）。
渐近方差：(\text{var}(\hat{A}_{MLE})\approx\frac{2A^2}{N(2A + 1)}=CRLB)，因此大样本下 MLE 是高效的（渐近高效）

MLE的渐进性

若数据的(p(\mathbf{x};\theta))满足一些 “正则性条件”，则未知参数(\theta)的 MLE 在大样本下，渐近服从正态分布：

$\hat{\theta}_{MLE} \stackrel{a}{\sim} \mathcal{N}\left(\theta, I^{-1}(\theta)\right)$

其中：

(\stackrel{a}{\sim})表示 “渐近服从”；
(I(\theta))是在真实参数值处计算的Fisher 信息（衡量数据中包含的关于参数(\theta)的信息量）；
(I^{-1}(\theta))是 Fisher 信息的逆，作为渐近分布的方差（体现估计量的精度）。

正则条件包括：

对数似然函数的一阶和二阶导数存在（可通过求导计算 Fisher 信息）；
Fisher 信息(I(\theta) \neq 0)（确保方差存在且合理）。

由渐近正态分布可推导出 MLE 的两个关键渐近性质：

渐近无偏：大样本下，(E\{\hat{\theta}_{MLE}\} \approx \theta)（估计量的期望趋近于真实参数）；
渐近高效：大样本下，MLE 的方差达到克拉美 - 罗下界（CRLB），是 “最优” 的估计量（在所有渐近无偏估计中，方差最小）。

举个例子来观察：假设对上面DC电平进行估计，真实的A=1，下图分别展示了N=5,20,50,100时，观测函数的直方图

可以看到，观测量越多时，直方图的peak越接近真实的1，因此当它是渐进无偏的。

多参数（向量）下的MLE估计

当待估计参数是向量(\boldsymbol{\theta})时，MLE 的核心逻辑与标量参数一致，依旧是寻找PDF的极值$\hat{\boldsymbol{\theta}}_{MLE} = \arg\max_{\boldsymbol{\theta}} p(\mathbf{x};\boldsymbol{\theta})$，只是需要用向量求导（梯度）来处理。只不过，在求$\frac{\partial \ln p(\mathbf{x};\boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \boldsymbol{0}$时，每个分量对应对(\boldsymbol{\theta})中单个元素的偏导

向量估计例子：求解WGN下的DC电平

假设WGN的方差是未知数，DC电平的振幅A是未知数，那么估计向量就是(\boldsymbol{\theta} = [A\ \sigma^2]^T)。单个样本的PDF是：

$p(x[n];A,\sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x[n]-A)^2}{2\sigma^2}\right)$

观测量为$x[n]$，N个样本的联合PDF为：

$p(\mathbf{x};A,\sigma^2) = \prod_{n=0}^{N-1} p(x[n];A,\sigma^2) = \frac{1}{(2\pi \sigma^2)^{N/2}} \exp\left(-\frac{1}{2\sigma^2} \sum_{n=0}^{N-1} (x[n]-A)^2\right)$

求$\ln p$：

$\ln p(\mathbf{x};A,\sigma^2) = -\frac{N}{2} \ln(2\pi) - \frac{N}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum_{n=0}^{N-1} (x[n]-A)^2$

分别对A和$\sigma^2$求偏导：

$\frac{\partial \ln p}{\partial A} = \frac{1}{\sigma^2} \sum_{n=0}^{N-1} (x[n]-A)$ $\frac{\partial \ln p}{\partial \sigma^2} = -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4} \sum_{n=0}^{N-1} (x[n]-A)^2$

令两个偏导数等于0，解出极值点：

$\frac{1}{\sigma^2} \sum_{n=0}^{N-1} (x[n]-A) = 0 \implies \sum (x[n]-A) = 0 \implies A = \frac{1}{N} \sum_{n=0}^{N-1} x[n] = \bar{x}$ $-\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4} \sum_{n=0}^{N-1} (x[n]-\bar{x})^2 = 0 \implies$ $-N\sigma^2 + \sum (x[n]-\bar{x})^2 = 0 \implies \sigma^2 = \frac{1}{N} \sum_{n=0}^{N-1} (x[n]-\bar{x})^2$

最终，向量参数(\boldsymbol{\theta} = [A\ \sigma^2]^T)的 MLE 为：

$\hat{\boldsymbol{\theta}}_{MLE} = \begin{bmatrix} \bar{x} \\ \frac{1}{N} \sum_{n=0}^{N-1} (x[n]-\bar{x})^2 \end{bmatrix}$

线性模型下MLE的最优性

对于前面介绍的线性模型(\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w})，且(\mathbf{w} \sim \mathcal{N}(\boldsymbol{0}, \mathbf{C}))，因此(\mathbf{x})的条件分布为：(\mathbf{x} \mid \boldsymbol{\theta} \sim \mathcal{N}(\mathbf{H}\boldsymbol{\theta}, \mathbf{C}))。尝试对这个线性模型使用MLE估计，将会看到它是最优的：

PDF为：

$p(\mathbf{x} \mid \boldsymbol{\theta}) = \frac{1}{(2\pi)^{N/2} \det(\mathbf{C})^{1/2}} \exp\left( -\frac{1}{2} (\mathbf{x} - \mathbf{H}\boldsymbol{\theta})^T \mathbf{C}^{-1} (\mathbf{x} - \mathbf{H}\boldsymbol{\theta}) \right)$

为简化计算，对似然函数取自然对数（对数似然函数）：(\ln p(\mathbf{x} \mid \boldsymbol{\theta}) = -\frac{N}{2} \ln(2\pi) - \frac{1}{2} \ln \det(\mathbf{C}) - \frac{1}{2} (\mathbf{x} - \mathbf{H}\boldsymbol{\theta})^T \mathbf{C}^{-1} (\mathbf{x} - \mathbf{H}\boldsymbol{\theta}))

MLE 的目标是找到使(\ln p(\mathbf{x} \mid \boldsymbol{\theta}))最大的(\boldsymbol{\theta})。由于前两项与(\boldsymbol{\theta})无关，只需最小化二次型项：(Q(\boldsymbol{\theta}) = (\mathbf{x} - \mathbf{H}\boldsymbol{\theta})^T \mathbf{C}^{-1} (\mathbf{x} - \mathbf{H}\boldsymbol{\theta}))

对(Q(\boldsymbol{\theta}))关于(\boldsymbol{\theta})求导并令导数为 0（极值条件）：

先展开(Q(\boldsymbol{\theta}))：

$Q(\boldsymbol{\theta}) = \mathbf{x}^T \mathbf{C}^{-1} \mathbf{x} - \mathbf{x}^T \mathbf{C}^{-1} \mathbf{H}\boldsymbol{\theta} - \boldsymbol{\theta}^T \mathbf{H}^T \mathbf{C}^{-1} \mathbf{x} + \boldsymbol{\theta}^T \mathbf{H}^T \mathbf{C}^{-1} \mathbf{H}\boldsymbol{\theta}$

对(\boldsymbol{\theta})求导（利用矩阵求导规则：(\frac{\partial (\boldsymbol{a}^T \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \boldsymbol{a})，(\frac{\partial (\boldsymbol{\theta}^T \mathbf{A} \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = 2\mathbf{A}\boldsymbol{\theta})，若(\mathbf{A})对称）：

$\frac{\partial Q}{\partial \boldsymbol{\theta}} = -\mathbf{H}^T \mathbf{C}^{-1} \mathbf{x} - \mathbf{H}^T \mathbf{C}^{-1} \mathbf{x} + 2\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H}\boldsymbol{\theta}$

令导数为 0，解方程：

$-2\mathbf{H}^T \mathbf{C}^{-1} \mathbf{x} + 2\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H}\boldsymbol{\theta} = \boldsymbol{0}$ $\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H}\boldsymbol{\theta} = \mathbf{H}^T \mathbf{C}^{-1} \mathbf{x}$

若(\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})可逆（实际中通常满足），则：

$\hat{\boldsymbol{\theta}}_{MLE} = (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{x}$

计算这个估计的协方差：

$\text{cov}(\hat{\boldsymbol{\theta}}_{MLE}) = (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \cdot \text{cov}(\mathbf{x}) \cdot \mathbf{C}^{-1} \mathbf{H} (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1}$

由于(\text{cov}(\mathbf{x}) = \mathbf{C})（因(\mathbf{x} = \mathbf{H}\boldsymbol{\theta} + \mathbf{w})，(\text{cov}(\mathbf{w}) = \mathbf{C})），代入得:

$\text{cov}(\hat{\boldsymbol{\theta}}_{MLE}) = (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1} \mathbf{H}^T \mathbf{C}^{-1} \mathbf{C} \mathbf{C}^{-1} \mathbf{H} (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1} = (\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H})^{-1}$

线性模型在一般噪声下的Fisher信息是$\mathbf{I}(\boldsymbol{\theta}) = -E\left\{ -\mathbf{H}^T \mathbf{C}^{-1} \mathbf{H} \right\} = \mathbf{H}^T \mathbf{C}^{-1} \mathbf{H}$，可以看到$\text{cov}(\hat{\boldsymbol{\theta}}_{MLE})$达到了$\mathbf{I}(\boldsymbol{\theta})^{-1}$，因此它是最优的。