克拉美-罗下界

估计量的衡量标准

对于参数估计问题,目前存在着很多估计算法。那么如何去衡量一个估计器(estimator, 也称估计量或估计算法)的性能,我们主要考量以下三个方面

  1. 无偏性(unbiased)。对于参数估计问题,设未知参数$\theta$,估计器模型$\hat{\theta}$。则有$\mathbb{E}[\hat{\theta}]=\theta$。对于估计对象为随机变量,则有$\mathbb{E}[\hat{\theta}]=\mathbb{E}[\theta]$。我们称满足这个条件的估计量为无偏估计量
  2. 有效性(availability)。有效性刻画估计量到真实值的偏离程度,$D(\hat{\theta})=\mathbb{E}[(\hat{\theta}-\mathbb{E}[\hat{\theta}])^2]$,即若存在多种无偏估计器,我们称估计量方差最小的估计器是有效的。
  3. 一致性(consistency)。设$\hat{\theta}$为未知参数$\theta$的估计量,若当样本数$N\rightarrow \infty$时,对于任意$\epsilon>0$,有$\lim\limits_{N\rightarrow \infty} P\left\{ {|\hat{\theta}-\theta|<\epsilon}\right\}=1$。我们称$\hat{\theta}$与$\theta$是一致的。一致性所体现的是,当样本总数逐渐增加时,估计量逐渐收敛于真实值 基于这三点考量,那么很自然我们会问,如何衡量一个无偏估计器是否是有效的。统计信号处理理论中的克拉美罗下界(Cramer-Rao Lower Bound,CRLB)就是衡量一个无偏估计器的有力工具。

克拉美-罗下界(Scale Parameter 标量参数)

对于估计参数$\theta$为标量时,假定PDF满足“正则”条件
\begin{align}
\mathbb{E}\left[{\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} }\right]=0\quad (\ \text{for any }\theta \ )
\end{align}
其中数学期望对$p(\boldsymbol{x};\theta)$取。那么无偏估计量$\hat{\theta}$的方差必然满足
\begin{align}
D(\hat{\theta}) \geq \frac{1}{-\mathbb{E}\left[{ \frac{\partial ^2\ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]}=\frac{1}{\mathbb{E}\left[{ \left(\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\right)^2}\right]}
\end{align}
其中导数是在$\theta$的真实值处求,数学期望是对$p(\boldsymbol{x};\theta)$取。因此,我们可以说一个无偏估计量$g(\boldsymbol{x})$达到CRLB,当且仅当,该估计量满足
\begin{align}
\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} =\mathbf{I}(\theta)(g(\boldsymbol{x})-\theta)
\end{align}
其中,$\mathbf{I}(\theta)=-\mathbb{E}\left[{ \frac{\partial ^2\ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]$,称为Fisher information。证明见博客第四部分。

Remarks: CRLB是衡量一个无偏估计器是否有效的重要工具,也就是说,给定一个无偏估计器,我们可以利用克拉美-罗下界去判断这个估计器是否是最优的。

Example:线性高斯模型(Linear Gaussian model)

\begin{align}
\boldsymbol{x}=\boldsymbol{h}\theta+\boldsymbol{w}, \quad \boldsymbol{w}\sim \mathcal{N}(\boldsymbol{0},\boldsymbol{C}_{\boldsymbol{w} })
\end{align}
其中$\theta$是未知参数,$\boldsymbol{x}\in \mathbb{R}^p$是观测值(observed signal),$\boldsymbol{w}$是均值为$\boldsymbol{0}$,协方差矩阵为$\boldsymbol{C}_{\boldsymbol{w} }$的高斯噪声。

我们考虑如下估计器
\begin{align}
\hat{\theta}=(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{x}
\end{align}
对于该模型,其似然函数$p(\boldsymbol{x};\theta)$为
\begin{align}
p(\boldsymbol{x};\theta)=\frac{1}{(2\pi)^{p/2}|\boldsymbol{C}_{\boldsymbol{w} }|^{1/2} } \exp \left[{-\frac{1}{2}(\boldsymbol{x}-\boldsymbol{h}\theta)^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}(\boldsymbol{x}-\boldsymbol{h}\theta)}\right]
\end{align}
因此

  1. 无偏性
    \begin{align}
    \mathbb{E}[\hat{\theta}]=\int_{\boldsymbol{x} } \hat{\theta} p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}
    \end{align}
    我们可以将$p(\boldsymbol{x};\theta)$看作为自变量为$\boldsymbol{x}$均值为$\boldsymbol{h}\theta$,协方差矩阵为$\boldsymbol{C}_{\boldsymbol{w} }$的高斯PDF,即$\int_{\boldsymbol{x} }\boldsymbol{x}p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}=\boldsymbol{h}\theta$。因此$\mathbb{E}[\hat{\theta}]=(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h}\theta=\theta$,即$\hat{\theta}$为无偏估计量。

  2. 有效性
    \begin{align}
    \frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}&=(\boldsymbol{x}-\boldsymbol{h}\theta)^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h}\\
    \frac{\partial^2 \ln p(\boldsymbol{x};\theta)}{\partial \theta^2}&=-\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h}
    \end{align}
    关于矩阵求导不太熟悉的朋友可以看下这个网站:https://en.wikipedia.org/wiki/Matrix_calculus
    基于上述表述,该系统模型的CRLB为
    \begin{align}
    -\frac{1}{-\mathbb{E}\left[{ \frac{\partial ^2\ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]}=\frac{1}{\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h} }
    \end{align}
    而估计器$\hat{\theta}$的方差为
    \begin{align}
    D(\hat{\theta})=\left({(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1} }\right) \boldsymbol{C}_{\boldsymbol{w} } \left({(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1} }\right)^T=(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}
    \end{align}
    由于$\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h}$是一维的,有$(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}=\frac{1}{\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h} }$,因此,该估计量是有效的,即该无偏估计量$\theta$的方差可以达到CRLB。

  3. 一致性
    将系统模型$\boldsymbol{x}=\boldsymbol{h}\theta+\boldsymbol{w}$代入估计器中,有
    \begin{align}
    \hat{\theta}
    &=(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}(\boldsymbol{h}\theta+\boldsymbol{w})\\
    &=\theta+(\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{w}
    \end{align}
    若假设噪声能量一定,即$\boldsymbol{C}_{\boldsymbol{w} }$元素值固定,随着观测样本$p\rightarrow\infty$,则噪声的方差
    \begin{align}
    D((\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{h})^{-1}\boldsymbol{h}^T\boldsymbol{C}_{\boldsymbol{w} }^{-1}\boldsymbol{w})=\frac{1}{\boldsymbol{h}^T\boldsymbol{c}_{\boldsymbol{w} }^{-1}\boldsymbol{h} }
    \end{align}
    从公式可以看出,假设噪声$\boldsymbol{w}$的每个元素具有相同的方差,则必然$\lim\limits_{p\rightarrow \infty}\boldsymbol{h}^T\boldsymbol{c}_{\boldsymbol{w} }^{-1}\boldsymbol{h}\rightarrow \infty$。因此,当$p\rightarrow \infty$时,我们可以将估计量$\hat{\theta}$看作
    \begin{align}
    \hat{\theta}=\theta+n,\quad n\sim\mathcal{N}(0,(\boldsymbol{h}^T\boldsymbol{C}_\boldsymbol{w}^{-1}\boldsymbol{h})^{-1}) \ \ \text{and} \ \lim\limits_{p\rightarrow \infty}\boldsymbol{h}^T\boldsymbol{c}_{\boldsymbol{w} }^{-1}\boldsymbol{h}\rightarrow \infty
    \end{align}
    因此,对于任意$\epsilon>0$,有
    \begin{align}
    \lim\limits_{N\rightarrow \infty} P\left\{ {|\hat{\theta}-\theta|<\epsilon}\right\}=1
    \end{align}
    即,该估计量满足一致性。

CRLB证明

由于$\hat{\theta}$是无偏估计,即
\begin{align}
&\int_{\boldsymbol{x} } (\hat{\theta}-\theta)p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}=0\\
\end{align}
等式两边对$\theta$求偏导有
\begin{align}
&\int (\hat{\theta}-\theta)\frac{\partial p(\boldsymbol{x};\theta)}{\partial \theta}\text{d}\boldsymbol{x}=1\\
\Rightarrow& \int (\hat{\theta}-\theta)\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}=1\\
\Rightarrow & \int (\hat{\theta}-\theta)\sqrt{p(\boldsymbol{x};\theta)}\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\sqrt{p(\boldsymbol{x};\theta)}\text{d}\boldsymbol{x}=1
\end{align}
由于柯西-施瓦茨不等式
\begin{align}
\int f^2(x)\text{d}x \int g^2(x)\text{d}x \geq\left({\int f(x)g(x)\text{d}x}\right)^2
\end{align}
当且仅当$f(x)=g(x)$时,取等号。

根据柯西-施瓦茨不等式(Cauchy-Schwarz inequality),有
\begin{align}
&\left({ \int (\hat{\theta}-\theta)^2{p(\boldsymbol{x};\theta)}\text{d}\boldsymbol{x} }\right)
\left({\int \left(\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\right)^2p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x} }\right)\geq 1\\
\Rightarrow &\int (\hat{\theta}-\theta)^2{p(\boldsymbol{x};\theta)}\text{d}\boldsymbol{x}\geq \frac{1}{\left({\int \left(\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\right)^2p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x} }\right)}
\end{align}

\begin{align}
D(\hat{\theta})\geq \frac{1}{\mathbb{E}\left[{\left(\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\right)^2}\right]}
\end{align}
现在只需证明
\begin{align}
\mathbb{E}\left[{\left(\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\right)^2}\right]=-\mathbb{E}\left[{ \frac{\partial ^2\ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]
\end{align}
证:由正则条件$\mathbb{E}\left[{\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} }\right]=0$,等式两边对$\theta$求偏导,有
\begin{align}
&\frac{\partial }{\partial \theta} \int \frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}=0\\
\Rightarrow& \int \left[{\frac{\partial^2 \ln p(\boldsymbol{x};\theta)}{\partial \theta^2}p(\boldsymbol{x};\theta)+\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta}\frac{\partial p(\boldsymbol{x};\theta)}{\partial \theta} }\right]\text{d}\boldsymbol{x}=0\\
\Rightarrow & \int \frac{\partial^2 \ln p(\boldsymbol{x};\theta)}{\partial \theta^2}p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}=-\int \left({\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} }\right)^2p(\boldsymbol{x};\theta)\text{d}\boldsymbol{x}
\end{align}

现在证明,若估计量$\hat{\theta}=\text{g}(\boldsymbol{x})$可以达到CRLB,则有
\begin{align}
\frac{\partial \ln p(\boldsymbol{x};\theta)}{\partial \theta} =\mathbf{I}(\theta)(g(\boldsymbol{x})-\theta)
\end{align}
其中,$\mathbf{I}(\theta)=-\mathbb{E}\left[{ \frac{\partial ^2\ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]$。

证:等式两边同时对$\theta$求偏导,有
\begin{align}
\frac{\partial^2 \ln p(\boldsymbol{x};\theta)}{\partial \theta^2} =\frac{\partial \mathbf{I}(\theta)}{\partial \theta}(g(\boldsymbol{x})-\theta)-\mathbf{I}(\theta)
\end{align}
等式两边同时对乘上$p(\boldsymbol{x};\theta)$,并对$\boldsymbol{x}$积分,得
\begin{align}
\mathbb{E}\left[{\frac{\partial^2 \ln p(\boldsymbol{x};\theta)}{\partial \theta^2} }\right]=-\mathbf{I}(\theta)
\end{align}
证毕。