我们继续接上一篇讲Kalman filter的基础预备知识。没有看过上篇的朋友，传送门在此：
写着写着发现东西还挺多的，为了避免文章过长。所以就把基础部分扩展成三篇吧。这篇算第二篇，我们主要讨论经典估计理论，下一篇我们主要讨论贝叶斯估计理论和Wiener filter。值得一提的是，我们应该从后者的角度更多地去理解Kalman filter——从他1960年写的关于KF的文章中，你能明显找到Minimum Mean Square Error Estimator(MMSE)的使用，而这正是Bayesian estimation里面一个常用的标准。而且状态估计和参数估计不是一样的，状态作为被估计量，并不是一个固定值，也是一个随机变量。
除了扩充了很多扩展参考资料供你们参考，Kalman本人的原文[6]被加入到了reference之中，感兴趣或者着急的朋友可以直接去看他文章中是怎么推导的。（一般都读的懂，前提是你懂概率，统计和信号处理等。）不过你要是没什么基础，就跟着文章一点点了解基础知识吧。有很多方法和思路可以引导我们最终得到经典的Kalman filter方程，但是我会把它放到后面几篇文章中去讨论。在此之前，我希望带领大家回顾和熟悉一些信号处理和经典的估计算法的内容作为知识的储备。主要介绍一些思路，帮助理解，他日深入的研究就会容易。如有笔误和疏漏等，也请诸位不吝指出。
才发现距离上次发文居然一个月了....喜欢的话多多点赞，可能能让我多抽点时间加快写写(你猜点赞后加速的概率有多大？)
本篇目录
1. Supplement to linear system theory and signal processing 线性系统和信号处理理论补充
2. Classic estimation theory 经典估计理论
- 2.1 Unbiasedness and consistency 无偏性与相合性
- 2.2 method of moments 矩法
- 2.3 Maximum Likelihood Estimation 最大似然估计(MLE)
- 2.4 Mean Square Error and Minimum-Variance Unbiased Estimator 均方差、最小方差无偏估计（MVUE）
- 2.5 Efficiency 有效性
- 2.6 Best Linear Unbiased Estimator 最佳线性无偏估计(BLUE)
- 2.7 Ordinary Least Squares 最小二乘算法 (OLS)
- 2.8 Weighted least squares 加权最小二乘算法 (WLS)
- 2.9 Recursive least squares 迭代最小二乘法 (RLS)
3. Summary 小结
1. Supplement to linear system thoery and signal processing 线性系统和信号处理理论补充
线性系统理论我在这个专栏里目前已经讲了很多了，一直关注我的朋友应该也已经了解这部分知识了。这里我们就不把state-space representation, controllability, observability这种基础知识再翻出来了，还不是很明白的朋友可以回去看之前的文章。
1.1 信号处理知识补充
这里要讲一些以前还没提及或者没有详细展开的知识。我们考虑简单一点的线性时不变系统：
这个系统的全解我们曾经说过，是以下面的式子来计算的：
前面一部分是初值影响的状态量，后一部分是输入引起的状态量。我们之前定义过LTI系统的状态转移矩阵(state transition matrix)为 ，所以(2)也就被重写成了：
我们以前也讲过LTI的零状态响应在时域中可以由以下convolution表示：
所以对比(3)和(4)，我们可以定义脉冲响应函数 . 注意卷积中
的积分变量
是可以和
的
互换而保持卷积值不变的。
现在假设系统的输入是一个random process，统一把输入改为宽平稳(Wide-sense stationary，WSS)随机过程 ，因为实际中我们多数信号都满足WSS要求。我们已经知道WSS随机过程的均值是常数，自相关函数(autocorrelation function) 只与延时
有关，而与起始时间无关，且二阶矩存在(为有限值)。两个信号
被称为jointly WSS, 如果它们都是WSS的，并且互相关函数(cross-correlation function)是时不变的，即只与信号之间的时间间隔有关。
注：我们再强调一下这个定义与stationary process的差别。 前者要求更弱，并不要求二阶矩以上的矩是否存在，只需要均值为常数，自相关函数只与时间间隔有关；后者要求过程的概率密度函数，包括一维和联合分布函数，都必须是严格时不变的，而对各阶矩的存在性没有任何要求。故而是平稳随机过程的期望，方差等等是不一定存在的（如柯西分布），但是如果存在的话，应当也是时不变的。
知道了 的统计特性和功率谱函数，下面来计算
经过系统
后的输出
的统计特性以及功率谱函数(power spectral density, PSD)：
我们考虑LTI系统的输出
输出的期望 满足：
注意在频域中就是频率特性或者传递函数
，利用终值定理可知其稳态时的Gain为
输出 和输入
的互相关函数(cross-correlation function)
:
注：不同定义中关于的定义有略微不同，同理也可以定义
，但本质上我们需要指定一个固定的公式，比如
，那么
。这里假设了互相关函数仅与时间间隔有关。(7)中使用了
的定义。
由于 是WSS过程，我们知道其自相关函数
只与时间间隔
有关(这是个dummy variable，只要当做是一个变量决定了整个函数值就行了)。因此上式也就是：
的自相关函数(autocorrelation function) 根据定义为：
我们发现：
借用(8)，我们知道：
根据维纳-辛钦定理(Wiener–Khinchin theorem)， 都是WWS random process，我们把(8)两端进行Fourier transformation，得到
与
的互功率谱密度函数(cross power spectral density, CPSD)
：
根据(11) 我们就得到了 的自功率谱密度函数(power spectral density)的两种表达式：
注：不就是所谓的我们在经典控制理论中的频率特性吗？换言之，也就是传递函数，或者在system identification领域叫frequency response function, FRF。(12)表明，如果我们知道
的PSD,
，以及频率特性或者传递函数
，那么我们就可以得到CPSD
。同时这个结论还有一个非常重要的应用，那就是传递函数估计。如果能通过序列估计出CPSD和PSD，我们可以通过：
来估算传递函数，这种FRF的估计方法称为estimator。有
，即采用(13)来计算FRF，这两种估计有不同的假设，即前者假设噪音与输入
是uncorrelated。(14)我们放着后面再解释。互谱和自谱的估计，我们在MATLAB中就有cpsd命令帮我们计算，而专门用于传递函数估计的tfestimate函数，也是调用了cpsd，以及采用了上面的公式作为其中一种估计方法。这是一种估计线性系统FRF的非常简单又实用的办法，只需要两个时间序列
注：下面是Stanford EE278的一些lecture notes，讲了更多细节方面的东西，之后有些不一定会用到，我就先略过了，供扩展阅读参考：
输出Y(t)的平均功率，满足各态历经假说的前提下，由二阶中心矩给出：
当我们把序列看成是一个n维实数向量 时，我们有自相关矩阵(auto-correlation matrix)是任意两个时间
之间定义的矩阵：
而自协方差矩阵(autocovariance matrix)定义为：
这两者的差别其实就在于autocovariance matrix减掉了 的期望。如果它的均值为0，即
，那么自相关矩阵和自协方差矩阵是等价的。至于为什么要写成矩阵的形式，那当然是因为矩阵运算有其独特的优越性。同理我们可以两个随机过程之间的互相关矩阵和协方差矩阵，这个很容易类推就不赘述了。更多阅读见[1,p116-119]。
1.2 随机线性系统、白化处理
一个连续的线性随机系统(linear stochastic differential equation)一般可以由以下方程表示：
这是一个时变线性系统，所以各矩阵都是时变的。其中 是测量输出，与我们以前定义的
，名义上的输出，还不一样（虽然这里还是用
矩阵符号了）。在第一行中的
就是一个噪音向量，我们考虑 zero-mean uncorrelated noise process，即均值为0的不自相关的噪音随机过程。显然高斯白噪声(AGWN)是满足条件的。
是噪声影响动态的耦合矩阵。而
是加入在测量中的测量噪声，我们同样假设是zero-mean uncorreleated的，同样也可以是AGWN。
而实际中我们经常会碰到离散的(discrete)线性随机系统，根据上面的结论，我们可以写出：
其中 为transition matrix，
，
.
为什么我们在模型中总假设噪音都是Gussian white noise呢？实际中的噪声常常都是有色噪声，即power spectral density并不总是为常数的随机噪声。上一篇中我提到了，当白噪声通过一个线性系统时，其输出是有色噪声。因此，已知有色噪声的PSD或者autocorrelation function，假定它是高斯白噪声通过某个线性系统后产生的输出，是否可以反向推出这个线性系统的模型呢？答案是肯定的。通过一个特定的线性系统，称之为"shaping filter"[1,p131],使得白噪声变成了需要的有色噪声，我们就可以把这部分动态的模型包含到原有系统的建模中。这使得我们可以把原来系统动态中的有色噪声，替换成白噪声和shaping filter的模型。
我们举一个时域连续系统的例子来说明。
当高斯白噪声通过一个一阶稳定线性系统时，其输出为一个自相关函数为指数函数的有色噪声。现在我们假设有色噪声的autocorrelation为：
其对应的PSD，满足宽平稳过程WSS的条件下，就是它的Fourier Transformation：
我们把这个有色噪声当做 ，把高斯白噪声当做
，服从标准正态分布，即方差
为1，均值
为0。那么
的自功率谱密度函数
，其自相关函数为
。根据公式(11)：
其中 代表了convolution运算。这个结论正是公式(14)的来源。我们把(22)两端进行Fourier Transform，得到了线性系统输出的自谱和输入的自谱关系：
那么现在代入 ，我们得到了：
取稳定的频率特性，便有了这个我们需要的shaping filter，即一个想象中假定存在的系统：
这就是把标准正态分布的高斯白噪声转化为autocorrelation为(20)，或者自谱为(21)的特定filter，我们可以把它写成时域微分方程：
我们说shaping filter (26)的输出 正是有色噪声
，而输入则是高斯白噪声
。注意这里我们替换了符号，采用了熟悉的
作为状态，
作为噪声。有了(26)之后，我们便可以对有色噪声进行替换，比如状态空间方程：
其中 是有色噪声，其autocorrelation满足(20)。那么我们 incorporate (26) into (27)，得到：
把(27)重新写成成新的状态空间表达式就不是什么难事了，所谓的augmented state space。我们这种将有色噪声输入转化为等效白噪声输入的过程也称之为白化处理。离散系统的做法并不一样，可以参考[1,p129, Example 4.4]，不过背后的想法是一样的。
2. Classic estimation theory 经典估计理论
估计是用过去到现在的数据来计算现在的某个量的大小。估计理论主要还是属于概率与统计学领域。我们知道统计学有两大学派：频率学派(frequentist statistics)和贝叶斯学派(bayesian statistics)。估计理论也可以根据这两个学派不同的观点分为经典估计理论和贝叶斯估计理论。 经典的估计理论我们在以前上统计学课程的时候应该已经接触过了，主要在参数估计那一块，包括点估计和区间估计两大类，其中有比如最大似然估计算法，一致最小方差无偏估计等。我们抽里面的重点进行回顾和讲解，目的是为了和后面的贝叶斯估计有一个对比。Kalman filter的推导用到的估计知识，主要是后者。对比两种估计理论的不同点，有利于加深记忆。统计推断(statistical inference)的任务是要通过研究样本数据来分析总体的特征。总体是由分布描述的，现实中很多时候容易知道分布的类型，但却不知道其参数。那么如何从样本数据中得到分布的参数是关键问题，即参数估计(parameter estiamtion)。值得注意的，参数估计和状态估计问题还不能完全等同起来，这点随着知识不断扩展慢慢体会。
2.1 Unbiasedness and consistency 无偏性与相合性
我们下面都以 来作为未知参数。频率学派认为，未知参数
的真值是一个确定的值(deterministic)。根据总体中的抽样样本
，构造一个统计量(statistic)：
称之为 的估计量(estimator)，或叫
的点估计(point estimator)，简称估计——点估计的名字是因为最后估计量的取值是一个数，与区间估计相对。统计量
是一个样本的函数(统计量的定义)，也是一个随机变量(因样本具有随机性)。
注：统计量本身是一个函数，但很多地方符号本身也表示了其取值，类似于
这么个意思。注意英语中的用词，estimator表示的是估计量，即统计量，其为样本的函数。estimate是estimator在某个样本处的取值，是根据取得的样本给出的估计值。
的估计可以构造无数种，估计之间孰优孰劣，需要一些标准来评判。我们说一个估计是unbiased(无偏的)，是指这个估计的期望值与对所有可能的被估计量
的真实值一样：
是
可能的取值空间，称参数空间。无偏性(unbiasedness)告诉了我们这个估计
是在真值
周围波动的，但如果我们求平均(期望)之后，最后等于真值。除非是
的线性函数，一般
的函数
是不会继承其无偏性的，比如修正样本方差
是无偏的，但是开平方根的
却是有偏的(biased)。 在其他评判标准满足时，存在无偏估计时，应当尽量使用无偏估计，这意味着估计从原理上就不存在偏差。有些估计，如样本方差
本身是有偏的，但随着样本数增加，它可以满足渐进无偏性(asymptotic unbiasedness)，在大样本时是可以近似认为无偏的。一个拥有较大bias的估计认为是有缺陷的，一般是不会采用的。不是所有参数都存在无偏估计，若
存在unbiased estimator，则称其是可估的(estimable)
随着样本量 的增加，我们希望估计量不断逼近参数真值，这种性质称为consistency(相合性)。把
看成是随着n变化形成的随机变量序列，则
依概率收敛到
。判断相合性一个常用的充分条件是：若
则
是
是consistent estimator. 相合性是统计量的一种大样本性质。无偏和相合是两个不同的概念。相合性被认为是很基本的要求，如果某些场合只存在有偏估计，那么至少也要采用一个相合的估计。
2.2 method of moments 矩法
常用和常见的一种估计方法是矩法(method of moments)，其思想使用样本矩(sample moment)去替换总体矩(population moment)，然后通过样本矩反推总体分布的参数值。其背后的思想基于格里汶科定理——经验分布函数随着样本增加接近总体分布。根据这个方法，即便有些时候分布形式是未知的，我们依旧可以估计一些参数：比如用样本均值 和样本方差
估计总体均值
(一阶原点矩
)与总体方差
(二阶中心矩
) 。在分布和概率函数形式已知时，我们可以写出该随机变量
的
阶矩表达式(以原点矩为例)：
，它们应当是概率函数中未知参数
的函数。用样本原点矩
代替
就可以列出
个方程来求解
。此法简单实用，不过不一定得到无偏估计量，但一般满足相合性。
2.3 Maximum Likelihood Estimation(MLE) 最大似然估计
MLE(最大似然估计，或极大似然估计)是非常传统的频率学派估计方法，当然学过了Bayesian estimation之后，我们又把它当作了一种maximum a posteriori (MAP) estimation的一种特例。MLE背后的思想非常简单易懂的：认为已经发生的事情是最有可能发生的。
相信各位已经看厌了投硬币的例子，那我编个让各位感同身受的例子：你通常喜欢在打游戏时反锁家门。上个月由于天天把家门反锁打游戏，被你爸发现骂的狗血淋头。这个月开始你洗心革面，决定不打游戏了，天天在家里看知乎学习。结果今天不小心你又把家里门给反锁了，傍晚爸爸回家看到了门被反锁，问他是相信你当天在打游戏还是没打游戏？
这个问题的题设有很多隐藏的信息需要解读。“当天家里门反锁状态”是一个随机变量 ，满足伯努利分布
，
是家里门反锁的概率。根据你打游戏喜欢锁家门的习惯，我们认为
的取值受到“当天你是否在打游戏”这个参数
的影响，
参数空间只有两个取值“打”或者“没打”。很明显，当天打没打游戏成为了当天家门反锁状态
分布列的一个至关重要的参数 。假设真实的概率情况：你当天打游戏99%概率门是反锁的，1%门不是反锁的。不打游戏则1%概率是反锁的，99%不反锁。那么真实的分布列只有两种可能：
当天打游戏 ——
="当天反锁”
；
="当天没反锁”
当天没打游戏 ——
="当天反锁”
；
="当天没反锁”
老爸当天晚上来到门前，这个参数 “打”还是“没打”已经是确定的事情了，只不过在他看来是未知的（其实你认真在看知乎）。今天门是否反锁的概率为
。他从上个月整整30天的观察中猜测了
与
之间的关系：打游戏，家里门一定反锁。这意味着他认为
"打"，
"当天反锁"
，据此此时可知，打游戏时门反锁概率
。而如果
“不打”， 他想了想，以前也没有遇到过其他门反锁的原因，于是认为
"当天反锁"
，此时
。他伸手去开门后，发现今天门被反锁了！综合以上，这两种情况，哪种可能性大一目了然！他微微一笑，晚餐的氛围当然不会很融洽——这小子又TM在打游戏。结果巧了，超小概率事情今儿发生了——咱们真的在认真看知乎（那你锁个毛线门啊，喂！）
爸爸的推断用到了最大似然估计的思想。我们假定一个随机变量 的分布函数包含
个未知参数需要估计，它们组成了向量
，它们都是确定值。所以这
取值必然也直接影响到
的
个样本
的联合概率密度函数，即
。将这一个joint probability density function看作是
的函数，我们把定义为likelihood function(似然函数)
：
最大似然估计就是取estimator满足似然函数的取值最大化：
老爸运用了MLE给了判断，不过运气不好小概率事情发生了！我们再次分一下整个过程：当天家门反锁情况 的真实分布是受到
这个参数的影响。当天你没有打游戏，所以
的真值其实是“no game today”。爸爸根据伯努利分布知道今天门被反锁的概率为
，马上检查了一下门，被反锁了！据此写出了样本的联合概率分布函数（这里只有当天的一个样本可检查，所以实际上只是单变量概率分布函数），亦为此时的似然函数
:
我们寻找一个 使得(32)最大化，则它就是我们给出的MLE。这里有一个很重要的问题，
与
之间的关系是什么？往往在丢硬币的例子中，我们把这个
当做了
，这是我们对确定分布形式得到的一个必然结果，因此我们只要选择大的
就可以令(32)最大化。这里“打游戏”和“当天家门是反锁的”的关系不是数值关系。对于既定的问题，我们必须要明确给出或者假定分布函数和参数的关系，才能让我们选择参数来最大化似然函数。爸爸在对
关系的猜测是基于过去30天经验给出的—— 所以他得到的结论是：
=“打”会让
，
="不打"，则反锁概率
。显然“今天儿子打游戏”这个判断能够使得“今天门上锁”这件事的概率达到最大。所以他自信推开门，想要开始教训你。很有意思的是，真实的情况，你打游戏时有99%的可能会锁门（剩下的1%不知道你在干什么锁门，难道是在看...？），所以真实的关系应该是，这样得出的概率
=“打”
，
="不打"，则反锁概率
。无论是那种情况，只要爸爸相信打游戏时门上锁的概率更大，那么他一发现当天门的反锁了，就会得出一样的结论：你今天打了游戏。
从条件概率的角度看，MLE试图找到最大化条件概率密度 的条件
。
这个例子与常见的丢硬币用MLE估计正面概率有微小区别。因为成功概率和参数
的关系是可能改变的，取决于爸爸如何看待，而不像投硬币时
是确定的。不管是什么关系，在做出MLE前，它们之间的关系应该是明确的，这样才有优化似然函数的方向，就像爸爸相信打游戏那天门反锁的概率远比不打时的大。这其实意味着，这个问题中“当天家门反锁”这个随机变量的分布形式不仅仅受到
关系的影响。这与一般MLE中我们假定分布形式确定，而参数未知还是有区别的。一旦
分布的具体形式。
贝叶斯估计中，认为参数，这是参数的先验分布函数。老爸可以根据过去经验猜测，预先猜测了今天我大概率在打游戏。我们后面会再提到。
最后总结一下，MLE逻辑是这样的：
- 有一个或多个固定的参数
影响随机变量
出现的概率。
- 既然如此，我们写出样本的联合概率分布函数或者PDF，那它就是样本和参数
- 对于我们抽样得到的某个样本，我们认为既然它发生了，它联合概率就应该是最大的。把这个样本代入到联合概率分布函数或者PDF后，它就只是参数
。
- 既然联合概率要最大，找到
使得
对数似然函数 用的比较多，因为它本身也是单调函数，且解析解求取比较方便。实际中似然函数多是用数值求极值的多，会用到牛顿法。根据似然函数对参数的一阶导数为0，MLE的解并不一定存在(解跑出参数空间之外)，也可能存在多解的问题。理论上MLE并不总是能给出无偏估计，不过仍可能是渐进无偏的。再加上MLE实际上需要你对统计模型有一个很好的建立，实际中的无偏也没法实现。MLE通常被认为是consistent的，虽然因为模型误差实际也不会得到，而且某些概率模型还真会让MLE不是consistent的[c]。但是consistency依旧是能保证就要保证的，即便是理论上的。关于MLE我们就说到这里，我们在后面的贝叶斯估计中可能还会call back。
关于样本的联合密度概率函数，样本抽样不一定要满足i.i.d，即独立同分布。当然很多时候满足的话，可以让我们直接将所有样本的概率密度函数直接相乘，是个不错的情况。我看很多人写文章都假定了i.i.d，尤其是讲机器学习的。故而提示一下，其实这不是一般的情况。
2.4 Mean Square Error and Minimum Variance Unbiased Estimator 均方差、最小方差无偏估计
估计量的unbiasedness和consistency并不意味着满足了我们对估计的全部要求。我们看来自[b]中的下图：
上面左图中的点可以代表一个参数，来自不同样本的无偏估计量分布。我们可以看到，取不同样本时，无偏估计的波动可以非常大，虽然我们可以大概看出其在红点(真值)附近比较均匀的波动的，即期望仍然是真值。右图是可以代表有偏估计的估计量分布。我们明显能看到这些估计量都十分靠近彼此，意味着估计量随着样本采样的变化波动很小，虽然我们一眼就看出来了其与真值之间总是存在某个方向的偏差。我想有学过机器学习基础课的朋友看到这里的图，应该一下就想起来了，这很像是讨论所谓的 Bias–variance tradeoff 问题，可能会用到accuracy来描述bias大小，用precision来描述variance大小。没错，统计学和机器学习本来就有很多内容是一样的，这里就是个很典型的例子。
评价点估计，除了使用无偏性之外，我们还引入了均方误差(Mean square error, MSE)来评价参数估计的优劣。我们发现有偏还是无偏，考虑的是估计误差 的一阶矩(期望)，而MSE考虑的就是
的二阶原点矩：
注意！ 此处我们假定了成为一个const了，于是也不会去讨论将MSE分成Var与Bias的平方两部分。事实上，
都是随机变量时，MSE的计算就要考虑它们之间的联合概率密度函数。
MSE的意义是明确的：估计参数 与真值
之间的距离的平方。当认为
是固定值时，它是由方差与偏差的平方两项组成的。同理，我们可以定义均方根误差，直接开平方就可以了。当估计量本身是无偏的时候，我们就有
，即估计量的均方差在估计无偏时就等于估计量的方差。MSE是一个非常一般的估计评判标准，我们自然希望MSE能够越小越好。 比较有偏的估计，方差和偏差的大小都应该考虑进去。这也启示我们：无偏估计虽然只有方差一项，但仍然可能使得MSE大于有偏估计。[see a page 324 ,example 6.4.1, b example 2.13] 此时有偏估计将比无偏估计在MSE评判下更优，或称之为更有效。
看到MSE，很多在学或者学过机器学习基础知识的同学可能会好奇，这里的MSE定义和我最后看到的形式有点不一样啊。比如在regression问题中检验data fitting好坏时常定义的MSE是：
其中是训练集或者测试集中
的target，
是根据predictor预测的值。通过计算预测值与target的MSE，我们评判predictor的优劣。所以实际上那里的MSE是求已知样本值与predictor预测值之间误差的二阶原点矩，或者说误差的平方的样本均值。理论分析一下，已知训练集或者测试集中的target值
外加uncorrelated噪声
，那么理论上MSE应该写成：
对比(54)与(52)，只是因为我们使用的不是真值，而是观测的值，所以产生了一项。(53) 采用误差平方的样本均值来近似计算MSE，是一种矩估计，是相合的。所以(53)这样的MSE定义是针对predictor的。给定训练集中的target的真值还是确定的，因此仍然有(52)那样分解得到Var和Bias的平方，这里只不过多考虑了噪声。
讲到这里你会发现这些都是统计学里的参数估计内容，而我们最终关心的还是偏状态估计。有三个问题我们要注意：。我们在统计学中定义估计量的特性时，通常会针对参数空间中任意的
确定——产生的样本去估计
![[公式]](https://www.zhihu.com/equation?tex=%5Ctheta%5Cin%7B%5B0%2C1%5D%7D)
时产生一组样本
，我们对此样本数据，使用
计算的估计量是无偏的。而
确定的一组样本
产生的估计量
也应该是无偏的。
既然有了MSE这个标准，我们自然想知道，从众多可能的estimator中是否能找到一个这样的 ？ 对
都，都有
，其中
是任意一个
的估计。换句话说
的对所有
可能取值的估计值，总有比其他任何estimator的MSE，即它在参数空间
中所有取值上都最小化了MSE。在频率学派的观点下，当
是确定值时，对所有可能形式的estimator，这样能够全局最小化MSE的estimator通常是不存在的！
这是一个很容易与稍后Bayesian MMSE混淆的地方，在那里我们明明定义了一个Minimum MSE的estimator！在经典估计理论中，因为我们并不可能知道真实的参数值，MSE中的Bias平方项的值将取决于真实的参数值。假设存在一个能够全局最小化MSE的estimator
，现在我们取一个特殊的estimator
使得其恒等于参数空间中的某一个参数
，此时对
，
。由于
。仔细一想，你就发现这是不可能的。这种最小均方差estimator在参数确定的情况下，只存在于规定形式的一类估计器中。
当我们只考虑unbiased estimator时，MSE中的Bias项就消失了，MSE就退化为估计量的方差。很多时候我们确实可以找到一个所谓Uniformly Minimum Variance Unbiased Estimator(UMVUE, or MVUE, 一致最小方差无偏估计，或最小方差无偏估计)——它相比于所有其他的unbiased estimators，对所有 的取值都有最小化的
。UMVUE并不总是存在的，有时候一个参数的unbiased estimator压根不存在；或者即便是存在的，也有可能没有一个estimator能够对所有参数实现最小Variance，如下图：
Cramer-Rao Lower Bound(CRLB,C-R下界) 指明了任何一个unbiased estimator的方差的下界，可以帮助我们有时确定UMVUE。为了得到CRLB，首先要介绍Fisher Information（费希尔信息量）。在满足一些Regularity conditions时(这里先不展开)，Fisher Information 被定义为：
其中 是log-likelihood function(对数似然函数)。此处我们没有进行任何样本独立同分布的假设。为了对(55)产生一点直观的理解，让我们考虑一个十分简单的参数估计问题：
，
是一个未知参数待估计，
是测量噪声，其满足正态分布
。如果我们只抽样了一个样本
，它应该满足正态分布
。那么其log-likelihood function是：
显然对(56)求导，我们得到：
令(57)为0，我们得到了 的MLE为
。如果我们继续求导，就可以得到：
根据(58)知道最大似然函数应该是一个开口朝下的曲线，要寻找最大值。考察函数的曲率，由于此时一阶导数(57)非常接近于0（事实上求期望可知E(x)=E(A)，即 ），那么二阶导数的绝对值就可以基本等同于曲率。我们发现在这个例子中，log-likelihood的最大值附近的曲率等于
。这意味着如果方差很大，即测量时噪声功率很大，那么log-likelihood曲线在最大值附近的曲率就小，导致整个函数看起来“又矮又平”；反之，噪声功率小，样本方差小，曲线在最大值附近看起来会非常陡峭。一个平缓的似然函数并有利于我们缩小参数真值的范围（自变量是参数
），这也是合乎常理的，因为噪声功率大了，自然就很难分辨真值了。所以我们看到
的二阶偏导的越大，表示目前这些样本能够更好地反映真值，即已经观察到的样本包含了更多的
的信息。可以证明最大似然估计MLE在特定的regularity condition下满足渐进正态性(asymptotic normality), 即随着样本量增加分布函数向正态分布
渐进，此时还满足相合性以及渐进无偏性。此时Fisher Information的倒数构成了它的渐进方差。一个回答，帮你理解：
让我们回到log-likelihood function ，对待估计参数
的一阶偏导被称为score function，可以记作
。可以证明满足概率函数满足某些运算条件时，
。score的方差被定义为Fisher Information：
可以看到Fisher Information也可以表示为似然函数二阶偏导期望的相反数，可以视作似然函数curvature(曲率)。
注意最后计算时这里的似然函数中的
取的观察到的样本,。
Cramer-Rao Lower Bound指出当 满足时，任何
的无偏估计
的方差必须满足：
并且以上不等式的等号成立，当且仅当 . 满足等号的unbiased estimator就是UMVUE，最小方差，即CRLB为
。对于i.i.d的样本，可以把总体密度函数写成各样本密度函数乘积，则可以定义单个observation的Fisher Information
, CRLB此时为
。
注意Fisher Information有些书上定义时采用population密度，有些假设了i.i.d采用了单个observation的密度。在茆诗松的书[a]上p329，他明明定义的population的FI，但是证明定理6.4.3时又采用了observation的FI，应该是个笔误。
我们可以看到有一些无偏估计确实达到了CRLB，比如用样本均值估计正态分布
的
。根据(57)有
. 所以
。而我们知道sample mean的方差
. 达到了CRLB，因此此时
是
此外，还有vector form的CRLB定理，适用于多于一个参数的情况，这里就不介绍了。既然已经知道了概念和用途，我们需要时自己去查阅资料就好了。
2.5 Efficiency 有效性
定义CRLB之后，我们可以更好解释efficiency这个概念。一个unbiased estimator的efficiency被定义为：
所以实际上就是理论最小方差比上实际方差，CRLB和efficiency给了一个判断unbiased estimator之间优劣的标准。当达到了CRLB，我们就称这个estimator是efficent的。更多关于efficiency参考wiki里的一些内容做一个简单了解[e]. 这些涉及细节定义上的东西，尽量不去多展开了，有需要自己做扩展阅读会比较好。
2.6 Best Linear Unbiased Estimator 最佳线性无偏估计(BLUE)
CRLB可能可以验证一些UMVUE(当它们达到CRLB时)，但是更多的时候这个下界是达不到的，一个UMVUE并不容易找到。如果我们能够把unbiased estimator的形式限制在样本的线性函数，那么我们便可以在linear unbiased estimator中找到UMVUE了。但是这种Variance的比较仅仅针对linear unbiased estimators之间而言，并不意味着它们达到了CRLB。所以BLUE的结果其实说是best，其实算是suboptimal的。
回顾MLE和一些矩估计中我们通常需要对概率密度 或者分布形式
有一个假设或者了解，概率模型的确定对结果会有很大影响。在上面计算Fisher information时也用到了这些信息。而对于BLUE，我们只需要知道
的first moment(一阶矩)和second moment(二阶矩)信息，连具体是什么分布都不需要知道，规避了完全掌握密度函数形式的要求，非常实用。
前面一直是从统计学的角度看待估计问题。我们尝试从信号处理的角度去求解BLUE。假设我们观测到了的一串受到噪音污染的信号序列，记为向量 。真实的无污染信号序列为，记为向量
。
的值受到
个参数组成的向量
的影响。
受到了噪声序列
的影响，其均值为
，协方差矩阵
。
的维度就是观测的数量
，待估计参数的数量记为
(
往往小于
)。对于第
个观测值，它们之间一般有如下量测关系：
其中 是给定的一般的非线性函数。我们的任务是通过
构造
的估计量。
为了得到 的BLUE，我们列出所需满足的条件限制：
- 对
的线性组合。
是所有Linear unbiased estimator中最小的，即BLUE中的Best的含义。
针对第一条：我们得到 的形式，其中
矩阵被定为
，
是每个参数
的线性项系数组成的列向量。故对每一个参数都有：
针对第二条：对 求期望
为了满足(64)，我们必须构造一个矩阵满足 。这样我们就有
这意味着 与
之间必须是仿射关系：
我们可以看到，x的随机性在这里完全是由噪声引起的。而
的值由于参数真值是固定的不存在随机性。我们看到噪声的协方差矩阵
就是测量信号的协方差矩阵
针对第三条：我们计算 的协方差矩阵：
协方差矩阵的对角线元素之和就是我们需要的方差之和。则我们要最小化 同时使得
。令
，则
对
求偏导得到：
所以BLUE已经找到：
的协方差矩阵此时变成了估计误差的协方差矩阵：
注意
整个过程求解BLUE的过程，我们用到了几个信息：
- 样本观察
- 我们需要掌握噪声的均值(一阶矩)与协方差矩阵（二阶矩），而不需要对噪声的分布和密度函数有所限定。
(69)的解是在我们认为噪声均值为0，协方差矩阵 已知时求出来的。如果我们假定噪声是符合正态分布的呢？假定此时我们知道噪声是高斯白噪声，其联合分布符合
。由于白噪声互不相关，于是满足i.i.d，其互协方差均为0。于是噪声
的协方差矩阵只有对角线上拥有非零元素，它们都等于噪声的方差
,即：
将(71) 代入(70)，则得到此时的BLUE拥有协方差矩阵：
我们知道CRLB告诉了我们UMVUE的最小方差，而在多变量时我们有Fisher Information Matrix(FIM)来告诉我们一个估计的协方差的下界：
。此时的。因而当噪声是高斯白噪声时，此时的BLUE(线性无偏最佳估计)：
并且它达到了CRLB，同时也是UMVUE！即没有其他任何无偏估计能比它产生更小的方差 Gaussian-Markov定理也表明了此时(74)的成立，但定理不保证UMVUE。所以UMVUE只是正态分布下的特殊结果。
细心的你可能已经发现了，只要噪声满足方差相同，不相关，均值为0，都会产生(74)的结果。没错！ 这就是Gaussian-Markov定理的内容。但是如果不满足正态分布的话，(74)就只是BLUE，并达不到CRLB。而不是UMVUE了！你可以仔细学习一下G-M定理的内容。
有色噪声的协方差矩阵还有对角线之外的元素。那么问题又来了，有色噪声不是白噪声怎么办？ 请复习前文：白化。对于高斯有色噪声，如果我们知道其协方差矩阵
2.7 Ordinary Least Squares 最小二乘算法
Ordinary least squares(OLS, or least squares, 最小二乘法)是最为经典的最优估计算法。实际上，经过前面一连串的分析，我们得到了一个有意思的结论：
Gauss–Markov theorem（高斯-马尔科夫定理）： 在线性回归模型中，如果误差满足零均值、同方差且互不相关，则回归系数的最佳线性无偏估计(BLUE, Best Linear unbiased estimator)就是普通最小二乘法估计。
什么是线性回归模型呢？请戳这里：https://en.wikipedia.org/wiki/Linear_regression。关键是因变量和参数为线性关系。
根据定理，我们可以立刻写出OLS的表达式：
正是(74)。但是与之前任何估计方法不同的是，它没有用到概率密度函数的任何信息，因此适用的范围就大大扩大了。我们常见的线性回归问题中，OLS是最基础且常用的。注意到(74)是BLUE的条件，那就是：
- 线性模型
- 零均值
- 同方差、且互不相关
我们在上面已经说过了，协方差矩阵此时就只有对角线上的数值，即各分量的方差。各方差都是相同的。定理中甚至都不同满足i.i.d条件。现在假设有高斯分布的噪声，我们知道这种情况下各条件均是满足的，于是就有了最小二乘法此时是BLUE的结论（虽然在高斯分布下还是UMVUE）。换句话说，如果此时噪声各分量方差不同，且还互相相关（比如有色噪声），那么OLS必然不是BLUE。
那BLUE是哪个呢？ 你想想？ 答案前面翻一翻。下面也提到了。
好了，下面就举一个OLS参数估计的例子。
OLS算法是不需要概率模型的，也不必一定要满足Gauss–Markov theorem的要求成为BLUE。它需要一个线性测量模型和足够且合适的样本观察。典型的应用是可以为一个方程数比未知数多的线性模型，称为overdetermined system，估计出近似最优解——这组解可以使得近似解的残差(residual)的平方和是最小的。举一个例子来说明。考虑一个参数估计问题：
这是一个典型的单自由度振动运动方程， 是质量(mass)，
是阻尼系数(damping coefficient)，
是刚度系数(stiffness)，
为作用在物体上的外力。假设我们手上有加速度、速度和位置的数据，即
，以及外力的信息
，试求未知常系数
。
从解方程的角度来看，这个问题非常容易，因为有3个unknowns，只要3组不同的数据建立3个方程就可以解出来了：
当然手上这些收集的数据其实都是有噪音的，所以任意选3组数据解出来的 值可都不一样。现在手上的数据严重冗余，比方说有
吧，可以任意选择其中3组来构成方程组(76)，每一种选择都有对应的一种估计
，那么问题来了： 哪一种选择是比较好呢？ 有没有办法把所有的数据都尽量用起来，不造成信息浪费呢？有的，比如把所有可能的解都求出来，求个平均。是不是很麻烦？ Least Squares就很好地解决了这个问题。考虑把所有的数据代入(75)得到矩阵方程：
记(31)为 . 我们定义一个参数估计，也就是(77)的近似解
，并定义其代入方程后与对应测量值
的残差，或者说误差：
其中 是
的第
行。求令所有
的平方和最小时对应的解，即令：
其中 是
的第
行，第
列的元素。(79)对
每一个元素
求偏导，并且使偏导为0，得到：
因为 , 那么就有：
等式的左边是 第
列的转置与
的乘积，右边是
第
列的转置与
的乘积。当
时，我们就有[1,p8]：
方程(82)称之为normal equation. 式(82)的解，当 是nonsingular时，为：
这时 是least squares意义下的最优参数估计。可见我们只要把加速度信号，速度信号和位置信号排列成矩阵
，外力测量值排列成向量
，我们就可以得到一个唯一的解，当然前提是
是可逆的。
的可逆性取决于
是否是列满秩的，实际上
也就是当
列满秩时的Moore-Penrose pseudoinverse(广义逆矩阵，或称之为伪逆矩阵)。实际遇到的矩阵
基本都是列满秩的(要是刚刚不列满秩，运气得有多好)
注：是可逆的，那么就有
，我们就得到了
的左逆为：
。这种超定方程
本来可以通过求逆
的方法来直接得到。对于非方的
，也就是(83)。所以在Matlab中，你可以用 y=inv(H'*H)*H^T*f 也可以直接用 y=pinv(H)*f, 结果是一样的。Matlab是通过SVD来求解的，这里的结论采用了满秩分解。关于线代不深究了，用到再说。
最小二乘法作为一种线性估计方法，其适用于估计值是观测值的线性函数。我们把之前的观测关系放到这里就能看到它们是一样的：
在这个参数估计问题中，在模型中加不加上噪声并不重要，因为我也没有噪声的任何信息可以利用。我们给出的根据误差平方和最小化的估计结果是：
注意到这里的优化标准和BLUE并不一样，误差平方和并不是MSE，因为我们没有求取期望。其实 中各误差分量的权重是一样，下面我们细讲。
2.8 Weighted least squares 加权最小二乘算法
如果我们观察OLS的指标：
(怎么定义 的符号并不影响，这里与我上面的例子中正负不相同)在最小化
的过程中，每一项误差的平方
在最小化过程中是有相同地位的，即它们的权重是一样的。这样的设定暗示着我们把每一个观测值的(决定估计量时的)重要程度都是一样的。然而不是所有时候观测值都是靠谱的；有多传感器融合时，不同观测值的统计特性有很大差别，比如方差大小不同。这一块我们在下一篇中写误差分析总结时会再次提到。总之，给不同的误差赋予不同的权重可以把这种差异主动考虑进估计结果之中，即产生了所谓的Weighted least squares(WLS, 加权最小二乘算法)。考虑WLS的指标为：
其中 是加权矩阵，
是
对角线上的元素。
矩阵的形式应该很容易想象吧？就是只有对角线元素的方阵。我们直接给出WLS的估计公式：
推导过程和LS是一样的，没有什么难度。
其实(88)我们已经很熟悉了，没错！ 当 时，它和BLUE的结果是一样的！这意味着权重矩阵并不是随意的，当权重矩阵和噪音的协方差矩阵互逆时，能使其成为线性无偏最优估计BLUE。
2.9 Recursive least squares 迭代最小二乘算法
传统的LS和WLS说起来算是一种“离线”估计，或者说批处理——每次都会把之前所有的观测值一起使用并直接给出估计值。计算机执行这种算法时，需要不断开辟存储空间来满足不断增加的观测量，以及矩阵计算会越加费力。在线估计的想法是，根据当前的观测值计算出现有的估计量，再根据新的观测值去更新原有的观测值，而不再需要重新使用历史中的观测量。这样大大降低了计算机的计算负担和存储负担，有利于算法的应用。LS和WLS都可以根据这种思想改为recursive least squares，即RLS算法。
考虑使用 个观测值进行WLS估计，我们有观测方程
：
其中 是矩阵
的第
行，
是行向量。考虑一个
维权重方阵
，根据(42)有WLS估计量:
其中 。考虑在此基础上，新添加一个观测方程，使用
个观测值进行WLS估计，同样有：
其中 。因为
,
,
，我们代入(45)中得到：
又因为其中:
即：
所以我们其实就是把两个公式进行了某些项的拆分并建立前后的联系。把(44)和(47）代入(46)，我们有：
最终有Recurive WLS 迭代公式为：
其中 为迭代的增益向量，会随着迭代的进行而发生改变最终收敛。相信做过优化的朋友应该对(50)的形式不会感到陌生，非常经典的形式。其中
正是采用前一次估计量
时产生的residual(残差)。
我们考虑一个十分简单的问题。如果有 个离散数据，我们尝试用一个常数
来估计这堆数据，也就是说我们假设这对数据是围绕一个常值
去波动的。那么每个量测方程就满足
很自然的，我们就知道这里的 。LS算法给出的估计结果就是：
就是所有观测量的平均数！如果我们考虑用RLS来计算的话，
根据(50）有最终结果：
可见这里的增益 会随着迭代进行不断衰减，最后收敛的结果和批处理的LS结果是一致的。
2.10 About other method 其他方法
还有不少算法其实在机器学习中经常出现，我可以对目标函数加以限时，得到LASSO或者ridge regression等等都是非常实用的参数估计算法。感兴趣的同学可以多搜集搜集，并且学着编编程序来实现他们。这里和KF关系有点远了，就不继续往外扩展了。
3. Summary 小结
本篇开始回顾了线性系统和信号处理的一些知识，重点是自功率谱函数和互功率谱函数，它们和FRF之间的关系，以及有色噪声的白化。
第二部分我们回顾了统计学中经典估计理论和众多实用的结果。我们从无偏和相合的基本概念出发，重点介绍了MLE，MVUE，BLUE以及最小二乘法和它的迭代算法。这一部分主要针对的还是参数估计，记住其重要的假设就是参数是确定的。这也是所说的为什么没有通常没有minimum MSE的estimator存在。在后续的贝叶斯估计中我们看到，参数也可以作为随机变量，我们可以给定或者猜想先验概率 。贝叶斯估计又将带领我们进入另外一个全新的世界，以及我们要去那里迎接Wiener filter——Kalman当初的文章也是为了解决Wiener Problem。
了解来龙去脉，方可游刃有余。
Main Reference
[1] Grewal, M. S., & Andrews, A. P. (2015).Kalman filtering: Theory and Practice with MATLAB Fourth Edition. John Wiley & Sons.
[2] Anderson, B. D., & Moore, J. B. (1979).Optimal filtering. Prentice-Hall, Englewood Cliffs, New Jersey.
[3] Gelb, A. (1974).Applied optimal estimation. MIT press.
[4] Musoff, H., & Zarchan, P. (2009).Fundamentals of Kalman filtering: a practical approach. American Institute of Aeronautics and Astronautics.
[5] 王可东(2019).Kalman滤波基础及MATLAB仿真，北京航空航天大学出版社.
[6] Kalman, R. E. (March 1, 1960). "A New Approach to Linear Filtering and Prediction Problems." ASME.J. Basic Eng. March 1960; 82(1): 35–45.https://doi.org/10.1115/1.3662552
Extended Reference for this article
[a] 茆诗松等，概率论与数理统计教程 第二版, 高等教育出版社
[b] 何志坚，数理统计讲义 数理统计讲义
[c] http://personal.psu.edu/drh20/asymp/fall2002/lectures/ln12.pdf
[d] http://web.eecs.umich.edu/~cscott/past_courses/eecs564w11/07_mvue.pdf
[e] https://en.wikipedia.org/wiki/Efficiency_(statistics)
[f] Bayesian Inferen 非常不错的学习基础概率与统计知识的网站
[g] http://engr207b.stanford.edu/lectures/mmse_estimation_2011_02_02_01.pdf
[h] ELEG 5633 Detection and Estimation
[j] https://www.ece.iastate.edu/~namrata/EE527_Spring08/l3.pdf
评论（0）
您还未登录，请登录后发表或查看评论