“如果你不能把它解释地足够简洁,说明你还不够理解它。” ————阿尔伯特·爱因斯坦

基本概念

以雷达探测一架飞机为例,我们可以使用$[x,y,z,v_x,v_y,v_z,a_x,a_y,a_z]$这样一组参数来描述飞机的状态,这一组参数成为系统状态
而我们又知道,由牛顿运动方程可知: \(x=x_0+v_{0x}t+\frac{1}{2}a_{0x}t^2\) \(y=y_0+v_{0y}t+\frac{1}{2}a_{0y}t^2\) \(z=z_0+v_{0z}t+\frac{1}{2}a_{0z}t^2\) 这一组方程描述了预测算法的输入和输出之间的关系,称之为动态模型状态空间模型
测量噪声:雷达显然无法绝对精确的测量出飞机的位置,其测量的噪声或不确定性成为测量噪声。
过程噪声:飞机受到风阻干扰,其状态转换并不总是完全按照状态方程,称之为过程噪声。

  考虑这样一个过程:一个人希望得到他的真实体重,他测量了五次,那么我们将他的真实体重成为隐藏状态(Hidden State),也即真实存在但是无法直接测得的状态。不过,虽然我们无法直接测得,但是我们仍然可以通过多次测量取平均值的方法来对这个真值进行估计。这个平均值称为期望。所谓期望,可以直观理解为对一个隐藏状态进行足够多次测量以后它所应表现出来的值。

方差

  显然,这个东西值得单独拿出来说。
  方差,是为了度量数据的分散程度,那么我们求出均值后,理所当然的想到去作差。但是我们发现,作差得到的有正有负,因为有的值大于均值,有的小于均值。
  那么我们理所当然的想到两个东西:平方,或者绝对值。选择平方有两个好处,其一,平方后能够扩大更加远离均值的点的影响;其二,平方可导,或者说绝对值并不是基本初等函数。
  说实话,这个理由是我到现在为止听的最充分的了。然后,把差值的平方加起来取均值,这就是方差!

\[\sigma^2=\frac{1}{N}\sum^{N}_{n=1}(x-\mu)^2\]

标准差的量纲和样本是一致的:

\[\sigma=\sqrt{\frac{1}{N}\sum^{N}_{n=1}(x-\mu)^2}\]

我们知道,方差还有一个公式,是这样的:

\[\sigma^2=\frac{1}{N-1}\sum^{N}_{n=1}(x-\mu)^2\]

关于这里为什么是N-1而不是N,我们不妨做一个详细的推导。考研的时候这一块我就没仔细看,我想,这里的推导应该是全网最详细的了,我每一步全都有注解,不可能有人写的比我还详细了。因为我当时就是看不懂这里的推导,死记硬背的,我不希望后面还有人再步我的后尘。
  首先,假定我们有10w个样本,这些样本作为总体的一部分,服从某种分布。约定$\mu$是总体的真实均值,$\overline{X}$是样本的均值,$S^2$是样本的方差,$\sigma^2$是总体方差。我们先抽样100次,计算得到均值$\overline{X_1}$,这似乎很好,但是相对于10w显然有点少了,于是我们又抽样了100个样本,计算得到均值$\overline{X_2}$,以此类推,多次重复后,我们得到了$\overline{X_1}, \overline{X_2}, \overline{X_3}, …, \overline{X_n}$,我们再对其计算均值,也就是期望,那么随着抽样次数的增多,最终$E(\overline{X})\to\mu$,这样的过程成为无偏估计。无偏估计的意义是:在多次重复下,它们的平均数接近所估计的参数真值。
  回到方差的讨论,假定我们的方差是这样的:

\[S^2=\frac{1}{N}\sum^{N}_{i=1}(x_i-\overline{X})^2\]

根据无偏估计的定义,也就是套上期望运算符,考察这个式子的期望是不是正在的方差。则有:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\overline{X})^2)\]

这一步就是把上面的S方套了个期望,然后我们在$(x_i-\overline{X})^2$里$+\mu-\mu$,有:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\overline{X}+\mu-\mu)^2)\]

整理为:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}((x_i-\mu)-(\overline{X}-\mu))^2)\]

至于为什么这样,我也不知道,应该是为了引入$\mu$。然后使用完全平方公式展开这个平方,有:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}((x_i-\mu)^2-2(x_i-\mu)(\overline{X}-\mu)+(\overline{X}-\mu)^2))\]

把$\frac{1}{N}\sum^{N}_{i=1}$给乘进去,有:

\[E(S^2)=E(\frac{1}{N} \sum^{N}_{i=1} (x_i-\mu)^2-2\frac{1}{N} \sum^{N}_{i=1}(x_i-\mu)(\overline{X}-\mu)+\frac{1}{N}\sum^{N}_{i=1}(\overline{X}-\mu)^2)\]

我们把$2\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)(\overline{X}-\mu)$单独拿出来看。我们知道,$\overline{X}-\mu$与求和符号无关,因为$\overline{X}$表示样本均值,当N固定的时候(只要N固定了),$\overline{X}$是一个固定的值,不论$i$遍历1到n都不会影响$\overline{X}$的取值。那么根据求和符号的性质,常数可以提到求和符号的外面来看:

\[2\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)(\overline{X}-\mu) = 2(\overline{X}-\mu)\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)\]

而后面的$\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)$又可以打开,有:

\[\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)=\overline{X}-\mu\]

所以:

\[2(\overline{X}-\mu)\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)=2(\overline{X}-\mu)^2\]

代回刚刚有三个求和符号的式子里面,有:

\[E(S^2)=E[\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2-2(\overline{X}-\mu)^2+\frac{1}{N}\sum^{N}_{i=1}(\overline{X}-\mu)^2]\]

最后一项同样与求和符号无关,求和符号和$\frac{1}{N}$抵消,有:

\[E(S^2)=E[\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2-2(\overline{X}-\mu)^2+(\overline{X}-\mu)^2]\]

即:

\[E(S^2)=E[\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2-(\overline{X}-\mu)^2]\]

根据期望的性质拆开,有:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2)-E((\overline{X}-\mu)^2)\]

考察这两项,显然第一项$E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2)$是样本的方差,但是又减去了一个什么,那么我们在令

\[S^2=\frac{1}{N}\sum^{N}_{i=1}(x_i-\overline{X})^2\]

时,$E(S^2)$竟然不是总体方差

\[\sigma^2=\frac{1}{N}\sum^{N}_{n=1}(x-\mu)^2\]

前面的$E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2)$是真实的总体方差,那么这个减掉的东西到底是什么呢?我们推导一下:

\[E((\overline{X}-\mu)^2)=Var(\overline{X})\]

这个是很显然的,你把方差公式里的N令为1就是这样。展开$\overline{X}$,有:

\[Var(\overline{X})=Var(\frac{1}{N}\sum^{N}_{i=1}{x_i})\]

根据方差的性质,提出来$\frac{1}{N}$,有:

\[Var(\frac{1}{N}\sum^{N}_{i=1}{x_i})=\frac{1}{N^2}Var(\sum^{N}_{i=1}{x_i})\]

然后这里需要注意,随机抽样的样本反映总体的特性,我的意思是,当总体服从方差为$\sigma^2$的时候,单个样本也服从方差为$\sigma^2$的分布。所以呢,

\[\frac{1}{N^2}Var(\sum^{N}_{i=1}{x_i})=\frac{1}{N^2}N\sigma^2=\frac{1}{N}\sigma^2\]

综合上面的推导,有:

\[Var(\overline{X})=\frac{\sigma^2}{N}\]

也就是,均值的方差也是服从总体的方差。那么考察刚刚推导中样本方差和总体方差有误差的部分:

\[E(S^2)=E(\frac{1}{N}\sum^{N}_{i=1}(x_i-\mu)^2)-E((\overline{X}-\mu)^2)\] \[=Var(X)-\frac{1}{N}Var(\overline{X})\] \[=\sigma^2-\frac{1}{N}\sigma^2\] \[=\frac{N-1}{N}\sigma^2\]

所以,要想让$E(S^2)$刚好等于$\sigma$,就需要上式两边同乘$\frac{N}{N-1}$,则有:

\[E(\frac{N}{N-1}S^2)=\sigma^2\]

其中

\[\frac{N}{N-1}S^2=\frac{N}{N-1}\frac{1}{N}\sum^{N}_{i=1}(x_i-\overline{X})^2\] \[=\frac{1}{N-1}\sum^{N}_{i=1}(x_i-\overline{X})^2\]

此即修正后的方差,其期望是等于总体方差的。 所以,当使用样本估计总体的时候,

\[S^2=\frac{1}{N-1}\sum^{N}_{i=1}(x_i-\overline{X})^2\]

其它基本概念

估计(Estimation):对系统的隐藏状态的估计。例如某飞行器的真实位置对观测者而言是不可见的,我们可以用传感器,例如雷达,来估计飞行器的位置。通过使用多个传感器并且使用先进的估计和跟踪算法(例如卡尔曼滤波)能够大幅提升对飞行器位置估计的效果。每次这样的测量和计算都是一次估计。
准度(Accuracy):描述测量值与真值的接近情况。高准度就是估计的均值接近真实值。
精度(Precision):描述一系列测量值相对同一个真值的偏差分布情况。高精度就是低方差。 有偏(Biased)系统:低准度系统。 无偏系统:测量(估计)值和隐藏状态完全相同。

α−β−γ滤波器

这类滤波器常用来对时间序列数据进行平滑。α−β滤波器、α−β−γ滤波器在原理上和卡尔曼滤波高度相关。

例1. 给金条称重

  本例对一个静态系统的状态进行估计。所谓静态系统,是指在合理时间范围内系统状态不会自发改变的系统。例如一座塔便是一个静态系统,高度便是其状态之一,它不随时间改变而变化。
  本例中,我们对金条进行称重。假定我们用来称金条的秤是无偏的,即称重结果没有系统性偏差,但是有随机噪声。
  首先,我们约定符号:
$n$为某个时刻
$z_n$是n时刻金条的测量值
$x$是金条重量的真实值
$\hat{x}{n,n}$是在n时刻,使用了n时刻的测量值,对x的估计值。
$\hat{x}
{n+1,n}$ 是在n时刻,对n+1时刻状态的预测。称之为外插。不过鉴于我们这个是静态系统,金条的重量并不会变化,所以显然 $\hat{x}{n+1,n}=\hat{x}{n,n}$ ,也就是说,金条在n+1时刻的重量预测值,等于当前时刻的估计值。

$\hat{x}{n-1,n-1}$ 是n-1时刻用了n-1时刻的测量值,产生的估计值。
$\hat{x}
{n,n-1}$ 是一个先验估计,也即在n-1时刻对n时刻的系统状态进行的预测。对于第n个时刻而言, $\hat{x}{n,n-1}$ 是先验估计,$\hat{x}{n+1,n}$ 是预测。

那么,首先想到的是,使用

\[\hat{x}_{n,n}=\frac{1}{n}(z_1+z_2+...+z_n)\]

进行估计当前时刻的系统状态。不过这需要存储历史上所有的测量值。那么我们不妨只存储上一时刻的估计值$\hat{x}_{n-1,n-1}$,并在新的测量完成之后更新一下。

\[\hat{x}_{n,n}=\frac{1}{n}\sum^n_{i=1}(z_i)\]

把这一时刻的提出来,写成之前的均值加上这一时刻的测量值

\[\hat{x}_{n,n}=\frac{1}{n}(\sum^{n-1}_{i=1}(z_i)+z_n)\]

为了推导状态转移方程,我们需要凑出上一时刻的估计值$\hat{x}_{n-1,n-1}$:

\[\begin{aligned} \hat{x}_{n,n} &= \frac{n-1}{n}\frac{1}{n-1}\sum^{n-1}_{i=1}z_i+z_n \\ &= \frac{n-1}{n}\hat{x}_{n-1,n-1}+\frac{1}{n}z_n \end{aligned}\]

整理一下,将$\frac{1}{n}$提出来,有

\[\begin{aligned} \hat{x}_{n,n} &= (1-\frac{1}{n})\hat{x}_{n-1,n-1}+\frac{1}{n}z_n \\ &= \hat{x}_{n-1,n-1} + \frac{1}{n}(z_n - \hat{x}_{n-1,n-1}) \end{aligned}\]

而系统模型是静态的,对第n时刻 x 的预测就等于n-1时刻对 x 的估计,有:

\[\hat{x}_{n-1,n-1} = \hat{x}_{n,n-1}\]

代入上面的式子,得到状态更新方程

\[\hat{x}_{n,n} = \hat{x}_{n,n-1} + \frac{1}{n}(z_n - \hat{x}_{n,n-1})\]

  关于这个例题,最后需要说的是,用这个公式对于静态系统,最终估算的结果是和平均数是一样的,但是我们改写成了递推式,这样可以少存储许多无用数据,减少许多计算。

  这个例题结束了,需要尤其注意的是,预测值和估计值之间的区别。在我的理解中,对于本例而言,也即对于一个静态系统,若当前时刻是n,那么因为$\hat{x}{n-1,n-1} = \hat{x}{n,n-1}$,第n-1时刻的估计值,就是对第n时刻的预测值,或者是先验估计值。而在n时刻得到了测量值$z_n$之后,经过某种计算公式推导得到的第n时刻的值,称为第n时刻的估计值。
  GPT给出的先验估计值的定义:这是基于前一个时刻(n-1)得到的对当前时刻(n)的估计值。它是在没有当前时刻观测数据的情况下,基于系统模型和先前状态进行的预测。
  总之,第n时刻的先验估计值,是不包含第n时刻的观测数据的,后验估计值是结合了当前时刻的观测数据的。
  第n时刻的先验估计值是基于第n-1时刻的后验估计值(已经结合了第n-1时刻的观测数据)通过系统模型预测得到的(这里静态系统就是相等),且不包含第n时刻的观测数据。因此称为“先验估计值”。而后在第n时刻的观测数据到达后,再计算出第n时刻的后验估计值。
  回到状态更新方程的讨论,这样的形式是卡尔曼滤波的五大方程之一,其中

\[K_n=\frac{1}{n}\]

  称为卡尔曼增益,符号为$K_n$
  状态更新方程的通式为:

\[\hat{x}_{n,n}=\hat{x}_{n,n-1} + K_n(z_n - \hat{x}_{n,n-1})\]

  其中$(z_n - \hat{x}_{n,n-1})$被称为“测量残差”,也叫更新量。更新量包含新的信息。
  我们还发现,第一次测量的权重是最大的,越往后随着n的增大,测量的权重会越来越小。

例2. 跟踪匀速直线运动的飞行器

  考虑如下情景:一个一维数轴上,有一个雷达站在坐标轴的零点,现有一架飞机正在朝一个方向匀速远离/接近雷达。
  那么,设n时刻飞机到雷达的距离为$x_n$,那么显然,微分之可得飞机当前的速度:

\[x'=v=\frac{dx}{dt}\]

  但是我们知道,雷达作为一个测量系统,其终究是有测量间隔的,就像无论多么高频率的示波器,也是对真实世界波形的采样。所以,设雷达相邻两次测量之间的时间间隔为$\Delta t$,那么由匀速直线运动的模型,我们可以得到以下公式

\[x_{n+1}=x_n+\Delta tv_n\] \[v_{n+1} = v_n\]

  上述方程成为状态外插方程,或者转移方程、预测方程。考虑这样一个情况,现在通过预测的位置发现并不准确,怎么办?可能是由两个原因:飞机并不是匀速直线运动,或者是雷达本身的测量误差。那么我们考虑使用状态更新方程修正速度:

\[\hat{\dot{x}}_{n,n}=\hat{\dot{x}}_{n,n-1}+\beta (\frac{z_n-\hat{x}_{n,n-1}}{\Delta t})\]

  也就是,下一时刻速度的预测值,应该是当前时刻的速度值,加上一个误差修正。这个误差修正是误差速度乘以一个系数得到的,所谓误差速度,就是第n时刻的实际位置,减去第n-1时刻预测的第n时刻的位置。为什么是这样一个误差速度?因为这个误差位置体现的是,使用n-1时刻的速度,得到的误差,但是我们现在发现了真正的位置并不是按照n-1时刻预测的那样,所以应当予以修正。
  关于修正系数$\beta$的取值,取决于你有多大程度相信你的位置测量值。例如,若你的雷达精度是20m,位置预测误差是90m,那么大概率是你的速度估计值出了问题,但是如果翻过来,那么可能就是系统误差,也就是测量值除了问题。这时候修正系数可以取小一点例如0.1。
  统共起来,有下面两个公式:

\[\hat{x}_{n,n} = \hat{x}_{n,n-1} + \alpha(z_n - \hat{x}_{n,n-1})\] \[\hat{\dot{x}}_{n,n}=\hat{\dot{x}}_{n,n-1}+\beta (\frac{z_n-\hat{x}_{n,n-1}}{\Delta t})\]

  其中的alpha和beta都是修正系数,根据雷达的测量精度调整。这个方程称之为α−β跟踪滤波方程。