53 – 贝叶斯滤波到粒子滤波

1、 状态估计

估计，就是从带随机噪声干扰的观测信号中提取有用信息。
估计分状态估计和参数估计两类：
状态估计是随时间变化的随机过程，
参数估计是不随时间变化或变化极慢的随机变量。

状态估计，就是根据受干扰的观测数据来估计关于随机变量或系统的特性。
像贝叶斯状态估计，及其派生的卡尔曼系列的滤波、粒子滤波都可以解决状态估计问题：
如果系统具有线性和高斯性，则考虑卡尔曼滤波；
如果处理非线性滤波问题，有扩展卡尔曼滤波；
如果非线性、非高斯，则考虑粒子滤波。

2、 最优估计

应用中，总希望估计出来的参数或状态愈接近真实值愈好。这个好得程度如何度量，就是最优估计问题。
最优估计是针对某一估计准则而言的，选择合理的估计准则很重要。
估计准则可分为非贝叶斯估计准则和贝叶斯估计准则：
前者非贝叶斯估计准则常见有最小二乘准则和极大似然准则；
后者贝叶斯估计准则常见有极大后验准则和最小方差准则。

2.1、非贝叶斯估计准则
非贝叶斯估计准则常见的有最小二乘准则和极大似然准则

2.1.1、最小二乘估计
由高斯提出：
对于非时变的n维随机向量X，对其进行k次线性观测得m维观测向量Y，有Yi=Hi*Xi+Vi，则可能由Y来估计X。对于一个X的估计X~，如果（Y-HX`）^t*（Y-
HX`）达到极小，则称X`为X的最小二乘估计或加权最小二乘估计。
注意点：必须保存和处理所有观测数据Y。
2.1.2、极大似然准则
由费希尔首次使用：
对于n维待估计量X，Y为X的k次独立同分布的观测数据，则定义似然函数为L（X）= ∏[ p(Yi|X) ]，如果样本集上的一个函数X`=f（Y1…Yk）使得L（X`）最大，则X`为X的极大似然估计。

2.2、贝叶斯估计准则
设X为待估计量，X的先验分布为 p(X)，k个观测值Y，Y的条件概率密度分布为p(yi|X)，则由贝叶斯公式知后验概率p(X|y1…yk)=p(y1…k|X)p(X)/p(1…k)，对此后验概率估计X`用于估计X的散失，定义为损失函数为L（x，x`），则此损失函数的期望E（L（x，x`））即贝叶斯后验风险，对于使得此风险最小的估计为贝叶斯估计。
根据损失函数的形式，有不同的贝叶斯估计结果，常见的有极大后验准则和最小方差准则。

2.2.1、最大后验估计
如果设损失函数为0-1风险，即损失函数L(x`,x)使得x`=x为1，x！=x则为0，则显然使贝叶斯后验风险达到最小即要求后验概率密度p(X|Y1…yk)最大，这种估计称为最大后验估计。
如果对X没有先验知识，最大后验估计就退化为极大似然估计，即极大似然估计是一种特殊的最大后验估计。
一般情况下由于考虑了X的先验分布），一般的最大后验估计优于极大似然估计。
2.2.2、最小方差准则
如果损失函数为二阶欧式距离，即L(x`,x)=(x`-x)^t*(x`-x)，则此准则为最小方差估计。

3、蒙特卡洛采样

采样的目的，是评估一个函数在某个分布之上的期望值。

例如抛硬币：
期望的结果，是符合一个出正面的概率为p反面概率为1-p的伯努利分布。
现实中，可以通过不断地进行抛硬币进行采样动作，从而评估这个概率p。
这个概率p，就是用出正面次数除以总次数的频率近似。

这个方法叫做蒙特卡洛。

4、蒙特卡洛积分

蒙特卡洛法就是：
对于一个分布来说，如果想要知道这个分布的期望，可以通过撒粒子的样本方式。

例如：
要求一个稀奇古怪的平面形状的面积。如果没有解析的表达方法，那就用蒙特卡洛法，就是在一个包裹了这个形状的更大图形例如外接圆或者外接正方形内随机的撒点，然后统计落进图形里的个数，当撒的点够多时，待求小面积/已知大面积=(落在图形内的点的个数/撒出去的点的总个数)，即可求出待求小面积=(落在图形内的点的个数/撒出去的点的总个数)*已知大面积。

再例如：
要求一个稀奇古怪的函数f(x)在ab区间的积分。如果没有解析办法或者解析方法很难，蒙特卡洛法就是在积分区间随机撒点，然后计算每个点xi的f(xi)值————这个是不费力的，当点足够多的时候， 点个数/(b-a)=∑[f(xi)]/∫[f(x)dx]，即ab区间的积分∫[f(x)dx]值为 (b-a)/点的个数*∑[f(xi)]。

推广一般：
要求随机变量x的期望 E(x)=∫[xp(x)]dx；
或者求某随机变量x的函数f(x)的期望 E(f(x))=∫[f(x)p(x)]dx；
而通常这个随机变量x的已知的分布函数p(x)形式非常复杂的，通常很难进行几分。
此时可以使用经典的蒙特卡洛方法：
比照这个分布，采样某些样本，取其平均，近似所求积分。
即用∑[ f(xi) ] —> ∫[f(x)p(x)]dx
使用这种蒙特卡洛方法计算的积分，就称为蒙特卡洛积分。

4、重要性采样

以上蒙特卡洛积分很好很优秀，可是存在一个关于采样的问题，怎么才能按照p(x)这个分布生成样本？

如果概率分布p(x)具有标准固定的形式，比如高斯分布，那是可以直接从p(x)中抽取粒子的。
但是通常p(x) 的分布未知，或者已知但是十分复杂。此时可以生成N个随机数ti…n，然后根据随机变量X的分布函数f(x)计算得到累计概率分布F(x)，通过每个F(xi)=ti就可以求得n个采样的样本点xi，i=1…n。 但这个方法，从随机变量x的pdf计算cdf这本身就是个几分问题。

所以此时需要其他的采样策略进行采样，例如重要性采样BIS：

如果由于变量X的分布函数f(x)不具有解析性、或者可计算性较差，那么可以引入一个已知的、形态较好的分布函数q(x)，这样就很容易求随机变量x或其函数f(x)在q(x)上的期望了。

变量x的期望：
E(x)
= ∫ [xp(x)dx] <---同乘q(x) = ∫ [x * p(x) * q(x)/q(x) dx] <---移位 = ∫ [p(x)/q(x) * x*q(x) dx] 函数f(x)的期望： E(f(x)) = ∫ [f(x)p(x)dx] <---同乘q(x) = ∫ [f(x) * p(x)q(x)/q(x) dx] <---移位 = ∫ [p(x)/q(x) * f(x)*q(x) dx] 可见，求p(x)分布下f(x)的期望，变成了，求q(x)分布下p(x)/q(x)*f(x)的期望。

重要性采样的关键，是把不好求的p(x)分布下f(x)的期望，经过重要性采样，变成了在另一个分布q(x)下相对好求的另外形式的期望。
就是：
由q(x)得累计概率分布Q(x)，然后用上面提到的采样方法通过Q(x)获得样本点xi…n，而p(xi)/q(xi)*f(xi)容易计算，q(xi)也容易计算，
那么∑[ p(xi)/q(xi)*f(xi) * q(xi)]， 就可以近似E(f(x))。

5、样本权重

重要性采样里面分布函数q(x)的选择，是个重要的问题。
如何能够尽量接近随机变量X的分布函数p(x)，需要根据具体的情况斟酌选择。
像前面图，从p(x)采样的粒子点对应的f(x) 数值都较小，那么在采样数量N有限的情况下很可能会导致无法获得f(x) 的值较大的那些样本，这样评估出来的期望值就导致偏大。

因此，可以：
通过在更大范围内按照均匀分布随机采样，然后从均匀分布随机生成的样本按照一种比例于p(xi) 的权重修正。

或者：
用某种距离/相似度来构造度量函数，然后再把这些距离通过某种方式变换成概率分布。
例如视频跟踪，在某给定区域内寻找特定的已知的人像模板，则粒子的状态是位置(x, y)和大小(w,h)，那么计算每个粒子权重的方法是：
对每个粒子，求其直方图的向量，然后和模板的向量求余弦距离t；
通过通过常用似然的指数族函数W(t)=K*exp^(-alpha*t)构造权重函数，距离远值小距离近值大，其中K和alpha保证累计和为1；
这样对于每个样本点就都有了一个权重值W(t)；
对此权重值归一化处理，使所有样本粒子的权重值之和为1，则权重值处理成概率值WB(t)；
这个WB(t)就是需要的q(x)。

或者：
像粒子滤波pf使用的是特征。
在机器人定位中，绿点机器人被随机放到已知地图的某位置，蓝线是其激光测量线，求自身在地图位置。
撒N个红点代表机器人对自己位置的猜想，即蒙特卡洛的样本粒子；这粒子的分布是均匀分布，也就是∫ [p(x)/q(x) * f(x)*q(x) dx]
里面的q(x)；那么如果简单把所有粒子进行平均，则机器人定位要永远处于地图中心点。不过有权重值p(x)/q(x)存在，可以修正不同的红点粒子。但是，这个权重比值里面，只有q(x)明了，而p(x)是关于机器人在地图中的位置的概率、这个作为定位问题的关键企却是不知道的。
粒子滤波器pf对于这里的不知道这个问题，不像上面那样，先选择某种度量，例如欧氏距离，然后构造似然，例如指数族函数=K*exp^(-alpha*t)；

它采用的是特征作为权重值来修正均匀分布的q(x)的：
地图已知、激光探测的环境已知，所以可以把p(x)/q(x)理解为：当前激光所探测到的环境的特征、猜测的地图红点附近的环境的特征、这两者的相似程度，这样p(x)/q(x)的计算转为特征相似度的计算。

为了更好的性能，自然希望p(x)/q(x)尽量接近1，也就是使蒙特卡罗预设的q(x)要尽可能的接近原始随即变量的概率密度p(x)，
这是重采样的工作意义。

6、 重采样