统计信号处理 (估计三)： Cramer-Rao下限

Cramer-Rao下限

未知参数的信息都是通过观测的数据以及那些数据的PDF得来的，所以如果PDF对参数依赖性较弱，或在极端情况下PDF与参数不相关，那么所估计的参数精度很差，如果PDF受未知参数影响大，所得到的估计越好。

例：依赖于未知参数的PDF。

如果观测到单个样本，即x[0]=A+w[0]x[0]=A+w[0]x[0]=A+w[0]，其中w[0]∼N(0,σ2)w[0]\sim N(0,\sigma^2)w[0]∼N(0,σ2)。

实际上，一个好的无偏估计是A^=x[0]\hat{A}=x[0]A^=x[0]，其方差刚好是σ2\sigma^2σ2，所以估计量的精度随着方差的减少而改善。

如上图所示，展示了两个具有不同方差的PDF与未知参数A的关系。可以看出(a)图的PDF更加集中，可以更为准确的估计未知参数A。

当把PDF看作观测值x已知，且是未知参数的函数时，称其为似然函数。
pi(x[0];A)=12πσi2exp[−12σi2(x[0]−A)2]p_i(x[0];A)=\frac{1}{\sqrt{2\pi\sigma_i^2}}exp[-\frac{1}{2\sigma_i^2}(x[0]-A)^2] pi(x[0];A)=2πσi21exp[−2σi21(x[0]−A)2]
根据上图可知，似然函数的"尖锐"性决定了所估计的位置参数的精度，所以为了定量的衡量似然函数的"尖锐"性。因此可以用对数似然函数在其峰值处的负的二阶导数，求出对数似然函数的曲率来衡量"尖锐"性。计算步骤如下：
p(x[0];A)=12πσ2exp[−12σ2(x[0]−A)2]∂lnp(x[0];A)∂A=1σ2(x[0]−A)−∂2lnp(x[0];A)∂A2=1σ2∴var(A^)=1−∂2lnp(x[0];A)∂A2p(x[0];A)=\frac{1}{\sqrt{2\pi\sigma^2}}exp[-\frac{1}{2\sigma^2}(x[0]-A)^2]\\ \frac{\partial lnp(x[0];A)}{\partial A}=\frac{1}{\sigma^2}(x[0]-A)\\ -\frac{\partial^2 lnp(x[0];A)}{\partial A^2}=\frac{1}{\sigma^2}\\ \therefore var(\hat A)=\frac{1}{-\frac{\partial^2 lnp(x[0];A)}{\partial A^2}} p(x[0];A)=2πσ21exp[−2σ21(x[0]−A)2]∂A∂lnp(x[0];A)=σ21(x[0]−A)−∂A2∂2lnp(x[0];A)=σ21∴var(A^)=−∂A2∂2lnp(x[0];A)1
所以估计量的方差随着曲率的增加而减少。而一般可以用对数似然函数的平均曲率来度量。
−E(∂2lnp(x[0];A)∂A2)-E(\frac{\partial^2 lnp(x[0];A)}{\partial A^2}) −E(∂A2∂2lnp(x[0];A))

求解Carmer-Rao下限

如果PDF p(x;θ)p(\pmb x;\theta)p(xxx;θ)满足其一阶导数的均值为0，
E[∂lnp(x;θ)∂θ]=0E[\frac{\partial lnp(x;\theta)}{\partial \theta}]=0 E[∂θ∂lnp(x;θ)]=0
那么任何无偏估计量θ^\hat \thetaθ^的方差必定满足
var(θ^)≥1−E(∂2lnp(x;θ)∂θ2)var(\hat \theta)\geq \frac{1}{-E(\frac{\partial^2 lnp(x;\theta)}{\partial\theta^2})} var(θ^)≥−E(∂θ2∂2lnp(x;θ))1
其中导数是在θ\thetaθ的真值处计算的，数学期望是对p(x;θ)p(\pmb x;\theta)p(xxx;θ)求取的，而且对于某个函数ggg和III，当且仅当
∂lnp(x;θ)∂θ=I(θ)(g(x)−θ)\frac{\partial lnp(x;\theta)}{\partial \theta}=I(\theta)(g(x)-\theta) ∂θ∂lnp(x;θ)=I(θ)(g(x)−θ)
时，对所有θ\thetaθ达到下限的无偏估计量就可以求得。这个估计量为θ^=g(x)\hat\theta=g(x)θ^=g(x)，它是MVU估计量，最小方差是1I(θ)\frac{1}{I(\theta)}I(θ)1。又因为
E[(∂lnp(x;θ)∂θ)2]=−E[(∂2lnp(x;θ)∂θ2)]E[(\frac{\partial lnp(x;\theta)}{\partial\theta})^2]=-E[(\frac{\partial^2 lnp(x;\theta)}{\partial\theta^2})] E[(∂θ∂lnp(x;θ))2]=−E[(∂θ2∂2lnp(x;θ))]
所以方差也可表示为
var(θ^)≥1E[(∂lnp(x;θ)∂θ)2]var(\hat\theta)\geq \frac{1}{E[(\frac{\partial lnp(x;\theta)}{\partial\theta})^2]} var(θ^)≥E[(∂θ∂lnp(x;θ))2]1
而公式5中的分母也被称为数据x的Fisher信息I(θ)I(\theta)I(θ),即
I(θ)=−E[(∂2lnp(x;θ)∂θ2)]var(θ^)≥1I(θ)I(\theta)=-E[(\frac{\partial^2 lnp(x;\theta)}{\partial\theta^2})]\\ var(\hat\theta)\geq\frac{1}{I(\theta)} I(θ)=−E[(∂θ2∂2lnp(x;θ))]var(θ^)≥I(θ)1

所以信息越多，CRLB越低，所以Fisher信息具有信息测度的基本性质：

由公式8可知，它是非负的。
对独立观测的可加性。

其中性质2可以得到一个结论对N个IID观察的CRLB是单次观察的1N\frac{1}{N}N1倍。而对于非独立的N次观察，得到的信息小于单次观察的N倍。

例：在高斯白噪声干扰中，信号为高斯白噪声的CRLB计算，可得观测值PDF
x[n]=s[n;θ]+w[n]p(x;θ)=1(2πσ2)N2exp{−12σ2∑n=0N−1(x[n]−s[n;θ])2}x[n]=s[n;\theta]+w[n]\\ p(x;\theta)=\frac{1}{(2\pi\sigma^2)^\frac{N}{2}}exp\{{-\frac{1}{2\sigma^2}\sum_{n=0}^{N-1}(x[n]-s[n;\theta])^2}\} x[n]=s[n;θ]+w[n]p(x;θ)=(2πσ2)2N1exp{−2σ21n=0∑N−1(x[n]−s[n;θ])2}
此时对其求一次导和二次导：
∂lnp(x;θ)∂θ=1σ2∑n=0N−1(x[n]−s[n;θ])∂s[n;θ]∂θ∂2lnp(x;θ)∂θ2=1σ2∑n=0N−1{(x[n]−s[n;θ])∂2s[n;θ]∂2θ−(∂s[n;θ]∂θ)2}\frac{\partial lnp(x;\theta)}{\partial \theta}=\frac{1}{\sigma^2}\sum_{n=0}^{N-1}(x[n]-s[n;\theta])\frac{\partial s[n;\theta]}{\partial \theta}\\ \frac{\partial^2 lnp(x;\theta)}{\partial \theta^2}=\frac{1}{\sigma^2}\sum_{n=0}^{N-1}\{(x[n]-s[n;\theta])\frac{\partial^2 s[n;\theta]}{\partial^2 \theta}-(\frac{\partial s[n;\theta]}{\partial \theta})^2\} ∂θ∂lnp(x;θ)=σ21n=0∑N−1(x[n]−s[n;θ])∂θ∂s[n;θ]∂θ2∂2lnp(x;θ)=σ21n=0∑N−1{(x[n]−s[n;θ])∂2θ∂2s[n;θ]−(∂θ∂s[n;θ])2}

取数学期望后，计算CRLB：
E(∂2lnp(x;θ)∂θ2)=−1σ2∑n=0N−1(∂s[x;θ]∂θ)2var(θ^)≥σ2∑n=0N−1(∂s[x;θ]∂θ)2E(\frac{\partial^2lnp(x;\theta)}{\partial\theta^2})=-\frac{1}{\sigma^2}\sum_{n=0}^{N-1}(\frac{\partial s[x;\theta]}{\partial \theta})^2\\ var(\hat\theta)\geq\frac{\sigma^2}{\sum_{n=0}^{N-1}(\frac{\partial s[x;\theta]}{\partial \theta})^2} E(∂θ2∂2lnp(x;θ))=−σ21n=0∑N−1(∂θ∂s[x;θ])2var(θ^)≥∑n=0N−1(∂θ∂s[x;θ])2σ2
CRLB中展示了信号s[x;θ]s[x;\theta]s[x;θ]依赖θ\thetaθ的重要性。如果信号随着未知参数的改变而迅速变化的话，那么其CRLB更小，能够产生更精确的估计量。

参数的变换

如果所要估计的量是某个基本参数的函数情况。例如对基本参数θ\thetaθ，想要得到估计量θ2\theta^2θ2的CRLB，如果估计量α=g(θ)\alpha=g(\theta)α=g(θ)，那么CRLB为
var(α^)≥(∂g∂θ)2−E[∂2lnp(x;θ)∂θ2]var(\hat\alpha)\geq\frac{(\frac{\partial g}{\partial \theta})^2}{-E[\frac{\partial^2lnp(x;\theta)}{\partial\theta^2}]} var(α^)≥−E[∂θ2∂2lnp(x;θ)](∂θ∂g)2
如果取α=g(θ)=θ2\alpha=g(\theta)=\theta^2α=g(θ)=θ2
var(A2^)≥(2A)2N/σ2=4A2σ2Nvar(\hat{A^2})\geq\frac{(2A)^2}{N/\sigma^2}=\frac{4A^2\sigma^2}{N} var(A2^)≥N/σ2(2A)2=N4A2σ2

扩展到矢量参数

如果将上述结果扩展到估计矢量参数θ=[θ1,θ2,⋯,θp]T\pmb \theta=[\theta_1,\theta_2,\cdots,\theta_p]^Tθθθ=[θ1,θ2,⋯,θp]T，假设估计矢量θ^\hat{\pmb \theta}θθθ^是无偏估计，那么矢量参数的CRLB可以对每个元素的方差放置一个下限，其中CRLB可以通过一个矩阵的逆的[i,i][i,i][i,i]元素求出，即
var(θi^)≥[I−1(θ)]iivar(\hat{\pmb\theta_i})\geq[\pmb I^{-1}(\pmb \theta)]_{ii} var(θθθi^)≥[III−1(θθθ)]ii
其中I(θ)\pmb I(\pmb \theta)III(θθθ)是p×pp\times pp×p的Fisher矩阵，其定义如下
[I(θ)]ij=−E[∂2lnp(x;θ)∂θi∂θj]i=1,2,...,p;j=1,2,...,p[\pmb I(\pmb \theta)]_{ij}=-E[\frac{\partial^2lnp(x;\theta)}{\partial\theta_i\partial\theta_j}]\quad i=1,2,...,p;j=1,2,...,p [III(θθθ)]ij=−E[∂θi∂θj∂2lnp(x;θ)]i=1,2,...,p;j=1,2,...,p