估计理论
估计理论是统计学和信号处理中的一个分支,主要是通过测量或经验数据来估计概率分布参数的数值。这些参数描述了实质情况或实际对象,它们能够回答估计函数提出的问题。
例如,估计投票人总体中,给特定候选人投票的人的比例。这个比例是一个不可观测的参数,因为投票人总体很大;估计值建立在投票者的一个小的随机采样上。
又如,雷达的目的是物体(飞机、船等)的定位。这种定位是通过分析收到的回声(回波)来实现的,定位提出的问题是“飞机在哪里?”为了回答这个问题,必须估计飞机到雷达之间的距离。如果雷达的绝对位置是已知的,那么飞机的绝对位置也是可以确定的。
在估计理论中,通常假定信息隐藏在包含噪声的信号中。噪声增加了不确定性,如果没有不确定性,那么也就没有必要估计了。
使用估计理论的领域
有非常多的领域使用参数估计理论。这些领域包括(当然不局限于以下列出的领域):
信号处理 X射线断层成像 脑电图 心电图 核磁共振 医学超声波扫描术 雷达、声纳、地震学——物件的定位 噪声方差 参数化(例如周期图和相关图谱)分析 非参数化(例如MUSIC、Root-MUSIC和ESPRIT)谱分析 维纳滤波 粒子滤波器 临床试验 民意调查 质量控制 通讯 信道参数 DC增益(请看下边的例子) 控制理论 卡尔曼滤波 随时间改变的执行器(英文:Actuator) 网络入侵侦查系统测量参数包含噪声或者其他不确定性。通过统计概率,可以求得最优化的解,用来从数据中提取尽可能多的信息。
估计过程
估计理论的全部目的都是获取一个估计函数,最好是一个可以实现的估计函数。估计函数输入测量数据,输出相应参数的估计。
我们通常希望估计函数能最优,一个最优的估计意味着所有的信息都被提取出来了;如果还有信息没有提取出来,那就意味着它不是最优的。
一般来说,求估计函数需要三步:
为了实现一个预测单个或者多个参数的所期望的估计器,首先需要确定系统的模型。这个模型需要将需要建模的过程以及不确定性和和噪声融合到一起,这个模型将描述参数应用领域的物理场景。 在确定模型之后,需要确定估计器的限制条件。这些限制条件可以通过如Cramér-Rao不等式这样的方法找到。 下一步,需要开发一个估计器或者应用一个已知的对于模型有效的估计器。这个估计器需要根据限制条件进行测试以确定它是否是最优估计器,如果是的话,它就是最好的估计器。 最后,在估计器上运行试验或者仿真以测试性能。当实现一个估计器之后,实际的数据有可能证明推导出估计器的模型是不正确的,这样的话就需要重复上面的过程重新寻找估计器。不能实现的估计器需要抛弃,然后开始一个新的过程。总的来说,估计器根据实际测量的数据预测物理模型的参数。
基础
对于给定模型,估计器需要若干统计 “成分”才能实现。第一,统计样本从长度为 N 的随机向量(英语:Multivariate_random_variable)(Random Variable,RV)中采样获得,观测值构成向量:
x = [ x [ 0 ] x [ 1 ] ⋮ x [ N − 1 ] ] . {\displaystyle \mathbf {x} ={\begin{bmatrix}x[0]\\x[1]\\\vdots \\x[N-1]\end{bmatrix}}.}第二,有 M 个参数:
θ = [ θ 1 θ 2 ⋮ θ M ] , {\displaystyle {\boldsymbol {\theta }}={\begin{bmatrix}\theta _{1}\\\theta _{2}\\\vdots \\\theta _{M}\end{bmatrix}},}它们的值需要被估计。第三,用于生成连续数据的概率密度函数(Probability density function,PDF)或离散数据的概率质量函数(Probability mass function,PMF)以参数值为条件(这些概率函数潜在存在),即条件概率为:
p ( x | θ ) . {\displaystyle p(\mathbf {x} |{\boldsymbol {\theta }}).\,}参数自身可能也存在概率分布(如贝叶斯统计),此时就需要定义贝叶斯概率:
π ( θ ) . {\displaystyle \pi ({\boldsymbol {\theta }}).\,}模型形成后,目标是估计参数,估计的参数通常表示为 θ ^ {\displaystyle {\hat {\boldsymbol {\theta }}}} ,其中 ⋅ ^ {\displaystyle {\hat {\cdot }}} 表示估计值。
常用的估计器包括最小均方误差(Minimum mean squared error,MMSE)估计器,它利用了估计参数和参数实际值之间的误差:
e = θ ^ − θ {\displaystyle \mathbf {e} ={\hat {\boldsymbol {\theta }}}-{\boldsymbol {\theta }}}作为优化的基础。该误差项平方的期望对MMSE估计器来说是最小的。
估计函数(估计子)
以下是一些相关的估计函数以及相关的主题
最大似然估计(Maximum likelihood estimation,简称MLE) 贝叶斯估计器(英语:Bayes_estimator)(Bayes estimator) 矩估计(Method of moments estimators,简称MME) Cramér-Rao界(英语:Cramér–Rao_bound) 最小二乘法(Least squares) 最小均方差(Minimum mean squared error,简称MMSE) 最大后验概率(Maximum a posteriori probability,简称MAP) 最小方差无偏估计(Minimum variance unbiased estimator,简称MVUE) 非线性系统识别(英语:Nonlinear_system_identification)(Nonlinear system identification) 最佳线性非偏估计(BLUE) 非偏估计,见偏差 (统计学)。 粒子滤波器(Particle filter) 马尔可夫链蒙特卡洛(Markov chain Monte Carlo,简称MCMC) 卡尔曼滤波 维纳滤波例子:高斯白噪声中的直流增益
考虑由 N {\displaystyle N} 个独立采样点构成的离散信号 x [ n ] {\displaystyle x[n]} ,它由常数 A {\displaystyle A} 和零均值、方差为 σ 2 {\displaystyle \sigma ^{2}} 的加性高斯白噪声 w [ n ] {\displaystyle w[n]} (即 N ( 0 , σ 2 ) {\displaystyle {\mathcal {N}}(0,\sigma ^{2})} )构成。方差已知,未知参数为 A {\displaystyle A} 。
x [ n ] = A + w [ n ] n = 0 , 1 , … , N − 1 {\displaystyle x[n]=A+w[n]\quad n=0,1,\dots ,N-1}参数 A {\displaystyle A} 的两个可能的估计器是:
A ^ 1 = x [ 0 ] {\displaystyle {\hat {A}}_{1}=x[0]} A ^ 2 = 1 N ∑ n = 0 N − 1 x [ n ] {\displaystyle {\hat {A}}_{2}={\frac {1}{N}}\sum _{n=0}^{N-1}x[n]} ,即采样平均(Sample mean)通过计算两个估计器的期望可以发现,它们的均值均为 A {\displaystyle A} :
E [ A ^ 1 ] = E [ x [ 0 ] ] = A {\displaystyle \mathrm {E} \left[{\hat {A}}_{1}\right]=\mathrm {E} \left[x[0]\right]=A}和
E [ A ^ 2 ] = E [ 1 N ∑ n = 0 N − 1 x [ n ] ] = 1 N [ ∑ n = 0 N − 1 E [ x [ n ] ] ] = 1 N [ N A ] = A {\displaystyle \mathrm {E} \left[{\hat {A}}_{2}\right]=\mathrm {E} \left[{\frac {1}{N}}\sum _{n=0}^{N-1}x[n]\right]={\frac {1}{N}}\left[\sum _{n=0}^{N-1}\mathrm {E} \left[x[n]\right]\right]={\frac {1}{N}}\left[NA\right]=A}两个估计器的均值没有差异,然而它们的方差不同:
v a r ( A ^ 1 ) = v a r ( x [ 0 ] ) = σ 2 {\displaystyle \mathrm {var} \left({\hat {A}}_{1}\right)=\mathrm {var} \left(x[0]\right)=\sigma ^{2}}和
v a r ( A ^ 2 ) = v a r ( 1 N ∑ n = 0 N − 1 x [ n ] ) = 1 N 2 [ ∑ n = 0 N − 1 v a r ( x [ n ] ) ] = 1 N 2 [ N σ 2 ] = σ 2 N {\displaystyle \mathrm {var} \left({\hat {A}}_{2}\right)=\mathrm {var} \left({\frac {1}{N}}\sum _{n=0}^{N-1}x[n]\right)={\frac {1}{N^{2}}}\left[\sum _{n=0}^{N-1}\mathrm {var} (x[n])\right]={\frac {1}{N^{2}}}\left[N\sigma ^{2}\right]={\frac {\sigma ^{2}}{N}}}当1}”> N > 1 {\displaystyle N>1} 1″ src=”/media/math_img/82695/001be86c3d8ec97f9b8194bb9ee5c282a5602d40.svg” usemap=”undefined” style=”width: 6.325ex; height: 2.176ex; vertical-align: -0.338ex;”>时,是一个更好的估计器。
最大似然估计
使用最大似然估计继续上面的例子,噪声在采样点 w [ n ] {\displaystyle w[n]} 上的概率密度函数(pdf)为:
p ( w [ n ] ) = 1 σ 2 π exp ( − 1 2 σ 2 w [ n ] 2 ) {\displaystyle p(w[n])={\frac {1}{\sigma {\sqrt {2\pi }}}}\exp \left(-{\frac {1}{2\sigma ^{2}}}w[n]^{2}\right)}此时 x [ n ] {\displaystyle x[n]} 的概率为( x [ n ] {\displaystyle x[n]} 服从分布 N ( A , σ 2 ) {\displaystyle {\mathcal {N}}(A,\sigma ^{2})} ):
p ( x [ n ] ; A ) = 1 σ 2 π exp ( − 1 2 σ 2 ( x [ n ] − A ) 2 ) {\displaystyle p(x[n];A)={\frac {1}{\sigma {\sqrt {2\pi }}}}\exp \left(-{\frac {1}{2\sigma ^{2}}}(x[n]-A)^{2}\right)}由于相互独立, x {\displaystyle \mathbf {x} } 的概率为:
p ( x ; A ) = ∏ n = 0 N − 1 p ( x [ n ] ; A ) = 1 ( σ 2 π ) N exp ( − 1 2 σ 2 ∑ n = 0 N − 1 ( x [ n ] − A ) 2 ) {\displaystyle p(\mathbf {x} ;A)=\prod _{n=0}^{N-1}p(x[n];A)={\frac {1}{\left(\sigma {\sqrt {2\pi }}\right)^{N}}}\exp \left(-{\frac {1}{2\sigma ^{2}}}\sum _{n=0}^{N-1}(x[n]-A)^{2}\right)}对上式取自然对数:
ln p ( x ; A ) = − N ln ( σ 2 π ) − 1 2 σ 2 ∑ n = 0 N − 1 ( x [ n ] − A ) 2 {\displaystyle \ln p(\mathbf {x} ;A)=-N\ln \left(\sigma {\sqrt {2\pi }}\right)-{\frac {1}{2\sigma ^{2}}}\sum _{n=0}^{N-1}(x[n]-A)^{2}}于是最大似然估计器为:
A ^ = arg max ln p ( x ; A ) {\displaystyle {\hat {A}}=\arg \max \ln p(\mathbf {x} ;A)}计算对数-最大似然函数的一阶导数:
∂ ∂ A ln p ( x ; A ) = 1 σ 2 [ ∑ n = 0 N − 1 ( x [ n ] − A ) ] = 1 σ 2 [ ∑ n = 0 N − 1 x [ n ] − N A ] {\displaystyle {\frac {\partial }{\partial A}}\ln p(\mathbf {x} ;A)={\frac {1}{\sigma ^{2}}}\left[\sum _{n=0}^{N-1}(x[n]-A)\right]={\frac {1}{\sigma ^{2}}}\left[\sum _{n=0}^{N-1}x[n]-NA\right]}令其为0:
0 = 1 σ 2 [ ∑ n = 0 N − 1 x [ n ] − N A ] = ∑ n = 0 N − 1 x [ n ] − N A {\displaystyle 0={\frac {1}{\sigma ^{2}}}\left[\sum _{n=0}^{N-1}x[n]-NA\right]=\sum _{n=0}^{N-1}x[n]-NA}得到最大似然估计器:
A ^ = 1 N ∑ n = 0 N − 1 x [ n ] {\displaystyle {\hat {A}}={\frac {1}{N}}\sum _{n=0}^{N-1}x[n]}它是一个简单的采样平均。从这个例子中可以发现,被独立同分布的加性高斯白噪声污染的、由未知常数构成的 N {\displaystyle N} 点信号的最大似然估计其就是采样平均。
Cramér-Rao下限
为了找到采样平均估计器的Cramér-Rao下限(CRLB),需要找到Fisher information数
I ( A ) = E ( [ ∂ ∂ θ ln p ( x ; A ) ] 2 ) = − E [ ∂ 2 ∂ θ 2 ln p ( x ; A ) ] {\displaystyle {\mathcal {I}}(A)=\mathrm {E} \left(\left[{\frac {\partial }{\partial \theta }}\ln p(\mathbf {x} ;A)\right]^{2}\right)=-\mathrm {E} \left[{\frac {\partial ^{2}}{\partial \theta ^{2}}}\ln p(\mathbf {x} ;A)\right]}从上面得到
∂ ∂ A ln p ( x ; A ) = 1 σ 2 [ ∑ n = 0 N − 1 x [ n ] − N A ] {\displaystyle {\frac {\partial }{\partial A}}\ln p(\mathbf {x} ;A)={\frac {1}{\sigma ^{2}}}\left[\sum _{n=0}^{N-1}x[n]-NA\right]}取二阶导数
∂ 2 ∂ A 2 ln p ( x ; A ) = 1 σ 2 ( − N ) = − N σ 2 {\displaystyle {\frac {\partial ^{2}}{\partial A^{2}}}\ln p(\mathbf {x} ;A)={\frac {1}{\sigma ^{2}}}(-N)={\frac {-N}{\sigma ^{2}}}}发现负的期望值是无关紧要的(trivial),因为它现在是一个确定的常数
− E [ ∂ 2 ∂ A 2 ln p ( x ; A ) ] = N σ 2 {\displaystyle -\mathrm {E} \left[{\frac {\partial ^{2}}{\partial A^{2}}}\ln p(\mathbf {x} ;A)\right]={\frac {N}{\sigma ^{2}}}}
最后,将Fisher information代入
v a r ( A ^ ) ≥ 1 I {\displaystyle \mathrm {var} \left({\hat {A}}\right)\geq {\frac {1}{\mathcal {I}}}}得到
v a r ( A ^ ) ≥ σ 2 N {\displaystyle \mathrm {var} \left({\hat {A}}\right)\geq {\frac {\sigma ^{2}}{N}}}将这个值与前面确定的采样平均的变化比较显示对于所有的 N {\displaystyle N} 和 A {\displaystyle A} 来说采样平均都是等于Cramér-Rao下限。
采样平均除了是最大似然估计器之外还是最小变化无偏估计器(MVUE)。
这个直流增益 + WGN的例子是Kay的统计信号处理基础中一个例子的再现。