Cramer-Rao下界

CRLB的作用和意义

在引入章节中介绍效率时已经介绍：克拉默-拉奥下界（Cramer-Rao Lower Bound, CRLB） 是无偏估计量方差的理论下界，它告诉我们在满足一定正则条件下，任何无偏估计量的方差都不可能低于这个界限。如果一个无偏估计的方差等于CRLB，则称其为最优的。因此，CRLB就是来评判估计还有没有优化的空间的。

虽然还有很多其他的方差下界，但 CRLB最容易确定，计算和推导相对简便。

如何求CRLB下界

（a）满足正则条件的PDF具有的性质

记一个随机变量的PDF为$p(\mathbf{x},\theta)$（又称为似然函数），这个$\theta$表示的是被估计值，是一个为标量。例如对方差为已知的$\sigma^2$的DC电平进行幅度估计时，$\theta$就是电平的振幅$A$。

如果$p(\mathbf{x},\theta)$满足正则条件（regularly conditions）的话，那么就会有：

正则条件（Regularity Conditions）目的是排除 “病态” 的概率模型。它的英文名称更直观，它其实就是限制在常规的范围内，这些”常规”通常包括：

  1. 可微性；
  2. 积分与求导的交换：若被积函数$p(x;\theta)$关于$\theta$的偏导数$\frac{\partial p(x;\theta)}{\partial \theta}在$积分区间内一致连续（或满足其他温和条件，常见概率分布都满足），则积分和求导的顺序可以交换，即$ \frac{\partial}{\partial \theta} \int_{a(\theta)}^{b(\theta)} f(x, \theta) dx = \int_{a(\theta)}^{b(\theta)} \frac{\partial f(x, \theta)}{\partial \theta} dx + f(b(\theta), \theta) \frac{db(\theta)}{d\theta} - f(a(\theta), \theta) \frac{da(\theta)}{d\theta} $
  3. 概率密度的规范性保持（$\int p(\mathbf{x},\theta)dx=1$）
  4. 支撑集与参数无关
  5. 参数空间的紧致性（即$\theta \in \Theta$，$\Theta$做为集合闭且有界）

关于$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$的证明如下：

已知概率密度函数有性质$\int p(\mathbf{x},\theta)dx=1$，现在对式子两边关于$\theta$求偏导。由于1求导等于0，因此得到：

根据求导与积分的交换性，上式可以写成：$\int_{-\infty}^{\infty} \frac{\partial}{\partial \theta}p(\mathbf{x},\theta)dx=0$ （这就来自于正则条件的基础）。

将期望内的函数进一步改写，因为 $\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta} = \frac{1}{p(\boldsymbol{x},\theta)} \cdot \frac{\partial p(\boldsymbol{x},\theta)}{\partial\theta}$（根据链式求导法则，$(\ln u)^\prime=\frac{u^\prime}{u}$，这里$u=p(\boldsymbol{x},\theta)$），通过移项得到：

而期望$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}$的求法是：

已知$\int_{-\infty}^{\infty} \frac{\partial}{\partial \theta}p(\mathbf{x},\theta)dx=0$，那么$\int_{-\infty}^{\infty}\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}\cdot p(\boldsymbol{x},\theta)dx=0$，故此证明了$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$

（b）由期望为零转换为方差下限

对于$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$，且$\hat \theta$为无偏估计量，可推导成：

推导过程如下：

对于无偏估计量$\hat \theta$，$E(\hat \theta)=\theta$。对它两边求导，就有：

将$E(\hat \theta)$展开，使用积分和求导的交换律：

在步骤（a）中推导过$\frac{\partial p(\boldsymbol{x},\theta)}{\partial\theta}=\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}\cdot p(\boldsymbol{x},\theta)$。将其代入其中：

首先将上式的$\hat \theta$写成$(\hat \theta + \theta - \theta)$的形式，来方便使用$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$消项：

对于上式后半截的$\int_{-\infty}^{\infty}(\theta) \frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}\cdot p(\boldsymbol{x},\theta)$，有：

因此：

将$p(\boldsymbol{x},\theta)$拆成两半，方便变成施瓦茨不等式的样子：

施瓦茨不等式（Schwarz’s Inequality）：（对于实数，可以去掉绝对值。）

取最大值的条件是：$f_1(x)=kf_2^*(x)$，即$f_1$是$f_2$的共轭的常数k倍。

将式子两侧平方，来凑齐施瓦茨不等式的样子：

根据不等式，有：

方差$Var(\hat \theta)=E\{[\hat \theta(x) - \theta]^2\}=\int_{-\infty}^{\infty}(\hat \theta(x) - \theta)^2p(\boldsymbol{x},\theta)dx$，

期望$E\bigg\{\bigg(\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg)^2\bigg\}=\int_{-\infty}^{\infty}\bigg[\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}\bigg]^2p(\boldsymbol{x},\theta) dx$

因此：

接着来处理$E\bigg\{\bigg(\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg)^2\bigg\}$。先通过$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$来建立联系，对它求导：

使用积分与求导的交换律：

对于复合函数求导，有：$\frac{\partial}{\partial\theta}(A\cdot B)=\frac{\partial A}{\partial\theta}B+\frac{\partial B}{\partial\theta}A$，那么上式就可以写成：

由于步骤（a）中推导的$\frac{\partial p(\boldsymbol{x},\theta)}{\partial\theta}=\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}{p(\boldsymbol{x},\theta)}$，上式右侧可以改变为：

因为$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$，因此$\frac{\partial}{\partial \theta}E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=\frac{\partial}{\partial \theta}0=0$，因此：

将上式写成期望的形式：

至此，就得到了：

（C）Fisher信息与CRLB

将$E\bigg\{\bigg(\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg)^2\bigg\}=-E\bigg\{\frac{\partial^2 \ln p(\mathbf{x},\theta)}{\partial \theta^2}\bigg\}$记为Fisher信息$I(\theta)$

由步骤（b）的不等式可以看出，最小的方差出现在不等式取等时，即CRLB是：

Fisher信息一定是非负的，因为$E\bigg\{\bigg(\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg)^2\bigg\}\geq0$

达到CRLB的无偏估计量

如果估计量$\hat\theta(x)$是通过估计函数$g(x)$得到的，那么，对于$g(x)$这个估计器，当且仅当其满足下式时，$g(x)$可以达到CRLB。

推导如下：

对于施瓦茨不等式，取等号时方差达到CRLB，而不等式成立的条件为$f_1(x)=kf_2^*(x)$，即$f_1$是$f_2$的共轭的常数k倍。在推导过程中，我们凑的不等式形式是：

由于这里$f_1,f_2$都是实函数，因此可以忽略共轭。那么，满足如下等式时，Var最小：

将两侧的$\sqrt{p(\boldsymbol{x},\theta)}$消掉，k移项，可以写成：

将$\hat\theta(x)$记作是通过估计函数$g(x)$得到的，那么上式就可以写成：$\frac{1}{k}(g(x) - \theta)$

单在施瓦茨不等式中，“常数k” 是证明施瓦茨不等式时引入的任意实数辅助量。然而这个地方还需要使得估计量的方差匹配 CRLB，即要满足$Var(\hat \theta)_{\min}=\frac{1}{I(\theta)}$，因此$k$并非可以为任意实数。这是由于在上面步骤（b）进行推导时，使用了步骤（a）的$E\bigg\{\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta}\bigg\}=0$来进行化简。这相当于为不等式引入了新的限定条件，而这个限定条件传导至结果就是$Var(\hat \theta)_{\min}=\frac{1}{I(\theta)}$，故，现在来确定这个常数$k$时，也需要把这个条件算进去。

这里的$k$是一个仅与$\theta$有关，与样本$x$无关的常数，先姑且记为$c(\theta)=\frac{1}{k}$，那么：

由于方差$Var(\hat \theta)=E[(\hat \theta-\theta)^2]$。因此我们对两侧平方后求期望来求方差：

将方差代入，同时$E\bigg\{\bigg(\frac{\partial \ln p(\boldsymbol{x},\theta)}{\partial\theta}\bigg)^2\bigg\}$为Fisher信息$I(\theta)$，即：

代入$Var(\hat \theta)_{\min}=\frac{1}{I(\theta)}$，来进行方差匹配

因此，只有当$\frac{\partial \ln p(\mathbf{x},\theta)}{\partial \theta} = I(\theta)[g(\mathbf{x}) - \theta]$时，$g(x)$的估计达到了CRLB。

在WGN中求CRLB并进行估计

0均值WGN中，对DC电平进行VMUE估计

有直流信号在信道传输，信道有均值为0，方差为$\sigma^2$的WGN $w$，采集信号为$x[n]=A+w[n]$。采集样本有N个，估计直流信号的电平值。

1.写出联合概率密度（似然函数）

被估计量为$A$，采集的每一个$x[n]$都是一个高斯随机变量，且他们相互独立。因此他们的联合概率密度为：（对于相互独立的随机变量，联合概率密度=各自的边缘概率密度相乘）

2. 求Fisher信息

对$p(\boldsymbol{x}; A) $求两次关于$A$的偏导，以此得到Fisher信息。

$\sum_{n=0}^{N-1} \left( x[n] - A \right)$可以写成N个x的平均值$N\times (\overline x-A)$的形式，因此有：

再次求导：

Fisher信息$I(\theta)=-E\{\frac{\partial^2 \ln p(\boldsymbol{x}; A)}{\partial A^2}\}$

此时的CRLB就是$\frac{1}{I(\theta)}=\frac{\sigma^2}{N}$

3. 通过以此偏导拿到估计函数

很明显，这里的MVUE估计函数$g(x)=\overline x$，即取采样值的平均值就是对其的MVUE估计。

0均值WGN中，对信号相位进行估计

有信号$Acos(2\pi f_0+\phi)$在信道传输，$A$和$f_0$为已知参数。信道有均值为0，方差为$\sigma^2$的WGN 。采集信号为$x[n]=Acos(2\pi f_0+\phi)+w[n]$，采集样本有N个，估计信号相位$\phi$。

1.写出联合概率密度

首先对观测的$x[n]$进行分析，$x[n]$应当是一组余弦序列叠加有高斯噪声。那么对单个采样点而言，它是一个均值为$Acos(2\pi f_0+\phi)$的高斯随机变量。

因此联合概率密度为：

2. 求Fisher信息与CRLB

先对$\ln p(\boldsymbol{x};\phi)$处理一下，方便后面求导：

然后对$\ln p(\boldsymbol{x};\phi)$求一阶偏导数。上式前半部分不含$\phi$，偏导直接为0。后半部分由链式法则：$F=f(g(x)),\partial F/\partial x=f’(g(x))g’(x)$，记$\left( x[n] - A\cos(2\pi f_0 n + \phi) \right) =g,\ f=g^2$。再使用求导与求和的交换律把求导变到求和内：

其中

代入原式：

再次求导，拿到Fisher信息，这里同样需要使用链式法则：

由于$x[n]=Acos(2\pi f_0+\phi)+w[n]$，$E\{x[n]\}=Acos(2\pi f_0+\phi)$，代入上式：

由于当N足够大时，$\cos(4\pi f_0 n + 2\phi)\approx 0 $，因此

CRLB为$\frac{1}{I(\theta)}=\frac{2\sigma^2}{NA^2}$

很明显，此时$\frac{\partial \ln p(\boldsymbol{x};\phi)}{\partial \phi}$不能被表示成$I(\phi)(g(x)-\phi)$的形式。因此通过此方法找到达到CRLB的MVUE，但是这并不代表MVUE不存在，而是需要通过其他方式找到

0均值WGN中CRLB的通式

假设信号的函数为$s[n;\theta]$，那么经过0均值方差为$\sigma^2$的WGN后，观测值为$x[n]=s[n;\theta]+w[n]$。

其联合概率密度函数为：

$\ln p$的一二阶倒数为：

Fisher信息为:

由于$x[n]-s[n;\theta]=w[n]$，均值为0。同时$ - \left( \frac{\partial s[n;\theta]}{\partial \theta} \right)^2$不含随机变量因此期望就是他自己。因此：

前面正则条件的约束才此时转化为：要求 $\frac{\partial s[n;\theta]}{\partial \theta}$是有限的

对正弦波的频率进行预测的CRLB

对于$s[n;\theta]=Acos(2\pi f_0n+\phi)$，$0<f_0<0.5$，假设幅度与相位已知，求预测$f_0$的CRLB

假设$N=10,\phi=0$，那么CRLB与频率$f_0$的关系如下图：

image-20250926164204008

这张图可以体现出信号对频率变化的响应强度。当$f_0$趋近于0或者0.5时，信号对频率变化的敏感性很低（曲线接近平坦）。此时，频率估计的方差下界会较高（因为分母的求和项小），意味着在端点处频率更难精确估计。

预测量多为参数时的CRLB（只是给个idea的程度）

假设被估计值$\theta$有多个参数，那么就要把它写成矩阵的形式：$\boldsymbol{\theta}=[A,B,C…M]^T$。其中ABC等等每一个都对应一个被估计的参数。假设$\theta$是一个M维的矩阵。

此时的似然函数$p$就变成了传入矩阵参数：$p(\boldsymbol{x}; \boldsymbol{\theta})$，$ \boldsymbol{\theta} \in \Theta$。

正则条件表示为（注意这个0是矩阵）：（求$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}}$时，是对每个元素分别求一次偏导，最后记录成M维的矩阵）

Fisher信息也变为矩阵的形式，维度为$M\times M$：

估计量的协方差矩阵$C_{\hat \theta}$满足：

其中第$i$个估计量的方差：

即，第$i$个元素的方差等于协方差矩阵对角线上第$ii$个元素（自己与自己的协方差）；它大于等于Fisher信息逆矩阵的第$ii$个元素。

可以找到一个估计函数矩阵$\boldsymbol{g(x)}$，如果满足下式，则$\boldsymbol{g(x)}$提供VMUE估计。

一个例子：直线拟合

假设有一条直线，它的采样收到了WGN的影响，有N个观测量：$x[n]=A+Bn+w[n]$。对AB为已知常数，根据观测量对这条曲线的参数A和B进行拟合。

此时的联合PDF：

求$\ln p$：

分别对A和B求$\ln p$的一阶偏导数：

因此$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}}$就是：

对它再次求偏导：

二阶混合偏导数连续时对称，即，先对 x 求偏导，再对 y 求偏导和先对 y 求偏导，再对 x 求偏导相等。而且这个只要所有的低阶偏导数连续，这个定理对更高阶也适用。

平方和求和公式：$\sum_{n=0}^{N-1} (n^2)=\frac{(N-1)N(2N-1)}{6}$

因此，$\ln p$的二阶偏导矩阵为：

对它挨个求期望，拿到Fisher矩阵

凑凑预测矩阵：

把这个矩阵切割成几个部分来写，试着凑凑$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \mathbf{I}(\boldsymbol{\theta}) [\mathbf{g}(\boldsymbol{x}) - \boldsymbol{\theta}]$的形式：

首先把AB分割出来，这样方便构成单独的$\theta$：

这构成了$\frac{\partial \ln p(\boldsymbol{x}; \boldsymbol{\theta})}{\partial \boldsymbol{\theta}}=\boldsymbol{h(x)}-\boldsymbol{I(\theta) \theta}=\mathbf{I}(\boldsymbol{\theta}) [\mathbf{g}(\boldsymbol{x}) - \boldsymbol{\theta}]$，那么：

使用预测矩阵估计出的CRLB

由于$\mathbf{I}^{-1}(\boldsymbol{\theta})= \sigma^2 \begin{bmatrix}
\frac{2(2N - 1)}{N(N + 1)} & -\frac{6}{N(N + 1)} \\
-\frac{6}{N(N + 1)} & \frac{12}{N(N^2 - 1)}
\end{bmatrix}$

对于两个估计量，估计出的CRLB满足：