Finance_时间序列分析_贝叶斯滤波&卡尔曼滤波

卡尔曼滤波的原理

https://space.bilibili.com/287989852?from=search&seid=10742836651514724908

https://www.bilibili.com/video/BV1QE411y7jb?from=search&seid=8397173496726970618

卡尔曼滤波配对交易

https://zhuanlan.zhihu.com/p/47471981

卡尔曼滤波pykalman的example：

一维：https://zhuanlan.zhihu.com/p/23670413

二维：

https://blog.csdn.net/coming_is_winter/article/details/79048700

卡尔曼滤波及其在配对交易中的应用: https://cloud.tencent.com/developer/article/1358311

https://xbuba.com/questions/55228609

总结：

1）卡尔曼滤波浓缩版

贝叶斯滤波：

贝叶斯滤波--> 卡尔曼滤波：

--》

卡尔曼滤波是一种利用线性系统状态方程，通过系统输入输出观测数据，对系统状态进行最优估计的算法。如图所示，x_{k_hat}是下一个状态的最优估计，每个最优估计都是有期望和方差，其中某些估计是很贴合真实的状态。但每次预测不能是特别准确的，需要用观测来优化我们的预测，这就是卡尔曼滤波总更新的操作。测量的方差代表我们测量的噪声，同样测量的部分数据也是在真实值里面的。所以我们可以利用预测和测量来找的更准确的估计。由于卡尔曼滤波的前提是高斯分布，所以我们可以融合这两个分布，得到更新后的下一状态的最优估计。

卡尔曼的核心思想就是：上面的那5个公式。

F是状态转移矩阵，μ是随机变量的期望，σ是随机变量x的方差/协方差矩阵，预测/转换方程为：

u_k^-=F*u_k-1⁺

(σ_k^-)² = F*(σ_k-1^-)² + Q_k　　Q~N(0,1) 若是协方差矩阵：Σ_k= F Σ_k-1F+ Q_k

我们预测时，通过F(状态矩阵）转换为新的状态，期望/均值 u 是不变的，但由于外界的不确定性，状态转移时会有偏差，这些不确定性就叫做噪声Q_k。

预测完之后，我们还需要用观测值来优化我们的预测。观测就需要用到传感器来观测数据。矩阵 Hk 来表示传感器的数据。同样，由于外界因素，传感器的信息也是完全准确，所以传感器读到的数据也有一定的方差。通过观测值我们可以大致猜到真实状态。

现在我们有了两个高斯分布，一个是预测值的高斯分布，一个传感器读数的高斯分布。这样我们就可以通过这两组数据/分布，来找到最优解。由于两个高斯分布相乘还是一个高斯分布。通过化简后的到一个相同部分K，这个K就是卡尔曼增益。

μ‘是更新后的均值，σ’是更新后的方差， μ₀和σ₀是预测的均值和方差，而 μ₁和σ₁是观测的均值和方差

K = σ₀²/σ₀²+σ₁² = hσ_k²/(hσ_k²+R) 　　R = σ₁² 是观测值的方差，hσ_k

μ‘= μ₀+K( μ₁-μ₀) = 》 μ⁺ = μ_k^-+K( y_k - hμ_k^-)

σ’ = σ₁²-Kσ₀² =》 σ_k⁺ = (1-Kh)σ_k^-²

2）卡尔曼滤波详细版

离散型随机变量的贝叶斯公式：P(X|Y) = P(Y|X)*P(X) /P(Y) 　　

连续型随机变量的贝叶斯公式：f(X|Y) = f(Y|X)*f(X) / f(Y) -- f(Y) =η --> η = 1/ ∫f(Y|X)*f(X) ，所以f(y)是常数η ，与三大概率无关。

后验概率 = η * 似然概率 * 先验概率 --> f(X|Y) = η*f(Y|X)*f(X)

预测方程：X_k = f(X_k-1) + Q

观测方程：Y_k = h(X_k) + R

贝叶斯滤波：

将预测方程带入到先验概率得到先验概率：P(X<x) = f_k^-(x) = ∫f_Q[x-f(v)] * f_k-1(v)dv
将观测方程带入到似然概率得到后验概率：P(X|Y) = f_k⁺(x) = η * ∫f_R[y_k-h(x)] * f_k^-(x)
从而得到常数η ：η = ( ∫f_R[y_k-h(x)] * f_k^-(x)dx )^-1

卡尔曼滤波：

1) 相对于贝叶斯滤波，多了几个假设：

f(x)与h(x)都是线性函数
f_Q(x)与f_R(x)都是服从正态分布

2）预测步（求先验概率）：根据贝叶斯滤波，第一步是写出预测步公式，然后将相应得正态分布的概率密度带入积分就行了。最后得到 u^-_k和σ^-_k。

3） 更新步（求后验概率）：已知 f_k^-(x) （先验概率）服从正态分布，那么f_k⁺(x)后验概率就是常数 η 乘以似然概率 f_R(y_k-h*x)（似然概率也服从正态分布）和先验概率f_k^-(x)。　　

最后得出卡尔曼滤波5公式：

卡尔曼滤波的应用：

模型1

# loding modules
import numpy as np
import pandas as pd
from pandas import *
import matplotlib.pyplot as plt

#X(k) = F*X(K-1)+Q
#Y(k) = H*X(K)+R
#第一个问题，生成一段随机信号，并滤波

#生成一段时间t，每隔0.01秒进行采样一次。
t = np.arange(0.1,1,0.01)
L = len(t)

#有了时间，就可以生成真实信号x，以及观测y。
#首先初始化
x=pd.Series(np.zeros(L))
y=pd.Series(np.zeros(L))    #不能将x赋值给y（即y=x），这样只会将x与y用做同一内存空间。
#生成信号，设x=t^2
for i in range(L):
    x[i]=t[i]**2
    y[i]=x[i]+np.random.normal(0,0.1)

plt.plot(t,x)
plt.plot(t,y)

#滤波算法
#预测方程怎么写：
#观测方程：Y(K)=X(K)+R R~N(0,1)，贝叶斯滤波也卡尔曼滤波的难点就是Q和R的方差的取值，它的难点在于没有规律性，只能不停的去试。
#预测方程不好写，画图（plt.plot(t,y)）后，可以看出曲线的一些规律性，可以猜一猜是线性增长，但是大多数问题，信号是杂乱无章的，怎么办。
#模型一，最粗糙的建模：
#X(k)=X(k-1)+Q
#Y(k)=X(k)+R
#猜Q~N(0,1)，弄成正态分布，相对好处理一些。
F1=1
H1=1
Q1=1
R1=1
#初始化X(k)+，要设置一个数组来装这个X(k)+，因为卡尔曼滤波是递推的，X(k)+需要X(k-1)+。
Xplus1=Series(np.zeros(L))
Pplus1=Series(np.zeros(L))
Xminus1=Series(np.zeros(L))
Pminus1=Series(np.zeros(L))
K1=Series(np.zeros(L))

#设置一个初值，假设Xplus1(1)~N(0.01,0.01^2)
Xplus1[0]=0.01
Pplus1[0]= 0.01**2

#卡尔曼滤波算法：它是期望与方差之间的递推关系，用前一时刻plus的期望和方差，推出这一时刻minus的期望和方差。以及这一时刻minus的期望和方差，推出这一时刻plus的期望和方差。再加上一个卡尔曼增益。X(0)是没有minus，且带plus的先验概率的，除了X(0)外所有X(k)的minus为先验概率，而plus为后验概率。
#X(k)minus = F*X(k-1)plus
#P(k)minus = F*P(k-1)plus*F+Q
#K = P(k)minus*H*inv(H*P(k)minus*H'+R)
#X(k)plus = X(k)minus+K*(y(k)-H*X(K)minus)
#P(k)plus = (1-K*H)*P(k)minus
for i in range(1,L):
    ##预测步：
    Xminus1[i]=F1*Xplus1[i-1]
    Pminus1[i]=F1*Pplus1[i-1]*F1'+Q1        #F1'是Matlab求转置。
    ##更新步：
    K1[i] = (Pminus1[i]*H1')/(H1*Pminus1[i]*H1'+R1)
    Xplus1[i] = Xminus1[i]+K1[i]*(y[i]-Xminus1[i]*H1)
    Pplus1[i]=(1-K1[i]*H1)*Pminus1[i]

plt.plot(figsize=(20,12))
plt.plot(t,x,'r')
plt.plot(t,y,'g')
plt.plot(t,Xplus1,'b')

#Q1,R1=0.01,控制噪声的方差。
#Xplus1[0]控制初始值，初始值的大小对滤波影响不大，最后滤波会调整。

模型1

模型2

###### 模型2 ##########
#模型2相对于模型1相比，观测方程还是一样，但预测方程一项，把X(k)写成用X(k-1)展开的一些项（泰勒展开式）。
#X(K)=X(K-1)*dt+X'(K-1)*dt+X''(K-1)*dt^2*(1/2!)
#Y(K)=X(K)+R R~N(0,1)
#此时状态变量不再是标量了，而是向量：X=[X(k) X'(K) X''(K)]T 列向量
#Y(K)=H*X+R H=[1 0 0] 行向量

#预测方程
#X(K)=X(K-1)+X'(K-1)*dt+X''(K-1)*dt^2*(1/2!)+Q2
#X'(K)=0*X(K-1)+X'(K-1)+X''(K-1)*dt+Q3
#X''(K)=0*X(K-1)+0X'(K-1)+X''(K-1)*dt+Q3 多项式信号多求多阶导数，总会比较平缓，而
#X'''(K)=X''(k-1)+Q3正是描述平衡的随机过程，这种建模相对精细一些，适用范围也较广。它的原理就是无论是什么函数，多求几次导，总会求到0，0就比较平缓了。对导数建随机平稳得模型，要比直接对多项式函数建平稳得模型相对来说精细一点。但有可能信号是指数，那么怎么建都不平缓。一般来说，信号不会增长太快，一般物理规律，信号不会增长太快。
#F= 1 dt 0.5*dt^2    #F是一个矩阵，在计算中要注意行向量和列向量。
#   0  1    dt
#   0  0    1
#H=[1 0 0]
#Q = Q2  0  0        #Q是一个噪声的协方差矩阵，Q协方差是用来描述他们随机变量的噪声的相互关系。由于状态变量X是由3个元素组成，Xk,Xk',Xk'',他们是多维随机变量，有各自的方差。而他们之间的关系，若两个随机变量的相关性为0，那么他们是不相关，若Xk增大了，Xk’是否会有变化，这需要用协方差来表述。
#    0   Q3 0
#    0   0  Q4
dt = t[2]-t[1]
F2 = np.mat([[1,dt,0.5*dt**2],[0,1,dt],[0,0,1]])    #此处要注意矩阵是否病态，若dt特别小，易导致矩阵病态或精度丢失
H2 = np.mat([[1,0,0]])
Q2 = np.mat([[1,0,0],[0,0.01,0],[0,0,0.001]])
R2 = 20            #只有一个观测，所以只有是标量，如果有多个观测，得写成矩阵形式。

#设置初始值
Xplus2 = np.matrix(np.zeros((3,L)))
Xplus2[0,0]=0.1**2
Xplus2[1,0]=0
Xplus2[2,0]=0
Pplus2=np.mat([[0.01,0,0],[0,0.01,0],[0,0,0.0001]])     #一般来说，初始得协方差矩阵都写成是对角矩阵形式。
#K2 = np.mat(np.zeros((3,1)))
for i in range(1,L):
    #预测步
    Xminus2[:,i]=F2*Xplus2[:,i-1]
    Pminus2=F2*Pplus2*F2.T+Q2
    #更新步
    K2=(Pminus2*H2.T)*np.linalg.inv(H2*Pminus2*H2.T+R2)
    Xplus2[:,i]=Xminus2[:,i] +K2*(y[i]-H2*Xminus2[:,i])
    Pplus2=(np.identity(3)-K2*H2)*Pminus2

Xplus2_0=np.zeros(L)
for i in range(L):
    Xplus2_0[i]=Xplus2[0,i]
plt.plot(figsize=(20,12))
plt.plot(t,x,'r')
plt.plot(t,y,'g')
plt.plot(t,Xplus1,'b')
plt.plot(t,Xplus2_0,'m')

# 可以进行在线滤波，实时滤波


### 问题2，两个传感器，进行滤波 ###
# Y1(K) = X(K)+R
# Y2(K) = X(K)+R
# H=[1 1]T （列向量） X=X(K)
# H=[[1 0 0],
     [1 0 0]]

模型2

模型2 中预测方程为：X(K)=X(K-1)+X'(K-1)*dt+X''(K-1)*dt^2*(1/2!)+Q2，可以将其分解为：

X'(K)=0*X(K-1)+X'(K-1)+X''(K-1)*dt+Q3
X''(K)=0*X(K-1)+0X'(K-1)+X''(K-1)*dt+Q3
X'''(K)=X''(k-1)+Q3

用矩阵表示为 X_k^- = F*X_k-1

#F= 1 dt 0.5*dt^2 　　#F是一个矩阵，在计算中要注意行向量和列向量。
# 　 0 1 　 dt
# 　 0 0 　 1

3）卡尔曼滤波的扩展

预测方程和观测方程的框架不变都是

预测方程：X_k = f(X_k-1) + Q

观测方程：Y_k = h(X_k) + R

在普通卡尔曼滤波时，假设f(Xk-1) 和h(Xk) 是线性的。f(X_k-1) = F*X_k 和 h(X_k) = H*X_k，其中F和H是系数。所以带入预测步，得到f(v) = F*v，以及带入更新步，得到h(x) = H*x。最终得到他们的期望，方差和增益，5个等式。

--》

而在扩展的卡尔曼滤波时，假设f(X_k-1) 和h(X_k) 是非线性的。f(x_k-1)=f'(x_k-1)X_k-1⁺f(x_k-1⁺)-f(x_k-1⁺)x_k-1⁺=AX_k-1+B，带入到了预测步，并利用软件计算

-软件->

EKF算法的总结：

设Xk-1服从N(xk-1+, Pk-1+)的正态分布。
预测步：先算A，再通过上一步的期望和方差，推出这一步的期望和方差
更新步：先算C，再算EKF增益，然后再算后验的期望和方差。

第二讲随机过程与主观概率

随机试验的定义，满足以下三个条件为随机实验：

　　1. 在相同条件下，试验可重复进行，即随机试验之间的独立性（如：抛硬币，第一次和第二次抛硬币之间没有任何影响，这样才能保证随机实验可重复进行，这样才可以通过多次实验来给概率赋值。若随机实验的独立性没有了，则无法再通过大数定律给概率赋值。）

　　2. 一次试验，结果不确定，所有可能的结果已知，

　　3. 试验之前，试验结果预先未知。所以无随机试验，就无法对概率赋值。

由于随机过程中的随机变量是不独立的，一般来说是无法做随机实验的，也就无法对概率赋值，但如果有了随机变量间的关系，且有了初始值，还是可以做随机试验的。即便如此，还得分情况，如果是与时间相关的随机过程是无法做随机试验，如股票，分子扩散，气温变化等（无法创造出随机实验的环境）。

概率分成两个学派：主观概率（贝叶斯学派）和大数定律（用频率代替概率，频率学派）

贝叶斯滤波的基本思想：引入外部的观测（证据、信息），来对主观的概率进行修正或更新。即主观概率（先验概率（先于实验之前的概率））通过外部观测（即贝叶斯公式），变为相对客观的概率，称这个概率叫做后验概率（实验之后的概率）

第三讲贝叶斯滤波的三大概率

先验概率靠猜（主观判断），后验概率靠测（测量结果推导出来），这两个概率分布之和都必须等于1，而似然概率则不需要，它只表示传感器的精度，与误差有关，而与真实温度。

1. 先验概率：没有做实验，通过学习等途径，得到的概率。包括主观概率（靠猜）

2. 后验概率：主观概率（先验概率（先于实验之前的概率））通过外部观测（即贝叶斯公式），变为相对客观的概率，称这个概率叫做后验概率（实验之后的概率）

3. 似然概率：似然概率代表着观测/传感器的精度，它衡量的是传感器的误差概率，叫似然概率。似然概率与先验/后验概率不一样，它在不同的真实值情况下的测量结果的概率无需加总为1。似然概率只是表示传感器的精度，只与传感器本身的误差有关，而与温度的真实值无关，不受其影响。例如：温度计的误差有可能是正负0.1，那么P(Tm=10.3|T=10)，真实温度为10的条件下，测量温度为10.3的概率就非常小，因为误差为正负0.1，那么真实值应该在10.2到10.4之间，所以Tm=10.3，而真实的值为10的概率很小。然而如果误差为正负0.5，那么P(Tm=10.3|T=10)就很大了。它是对两个不同真实值，进行不同的测量，他们两个不同的事件，所以他们之间的概率没有任何关系。

4. 最后，随机变量的取值与随机变量的分布概率无关，这对P(Tm=10.3)是一个常数起着至关重要的作用。因为P(Tm=10.3)，即测量值为10.3是一个结果，随机变量的概率分布对其无关（即，抛硬币，正面/反面朝上只是一次抛硬币的结果，与其概率分布（0.5）无关）。P(Tm=10.3)可以通过全概公式得到，如果先验概率定了，且似然概率也定了，那么P(Tm=10.3)的概率也就定了，是一个常数。

既然P(Tm=10.3)是一个常数（设这个常数为η），我们可以把上面这个贝叶斯公式进行改写：后验=η*似然*先验。所以常数η=1/Σ先验*似然。

第四讲连续随机变量的贝叶斯公式

状态看成原因，观测看成结果，要先有状态，才能有观测，所以说状态是因（温度，Y），而观测是果（测量值，X）。

后验概率就是由果推因（观测值->真实值的概率），看到了观测得结果，去推出最有可能的原因是什么？观测值是10度，那么它最有可能的状态（真实值），是11度的概率多大，是10度的概率是多大，还是9度的概率多大呢？

似然概率是由因去推果（真实值->观测值的概率），真实的温度是10度，看导致观测值是10度的概率，观测值是11度的概率呢？

连续随机变量的贝叶斯公式：

先验概率：f_x(x)

似然概率：f_y|x(y|x)，由因（X）推果（Y）

后验概率：先验概率*似然概率 = f_x|y(x|y) = f_y|x(y|x) * f_x(x) / f_y(y）= f(x,y)/f(y)

f(x,y) 的边缘密度为：f_x(x) = ∫ f(x,y)dy, f_y(y) = ∫ f(x,y)dx。

第五讲似然概率与狄拉克函数

引例：测温度

第一步，先猜一个先验概率密度，猜它是服从N(10,1)的正态分布。f_x(x)~N(10,1)，其中，期望是对温度的猜测，即猜测10度的概率最大，而方差是对主观猜测的误差。通过温度计测的的观测值为9度：y=9。如果有了先验概率，再得到似然概率，根据贝叶斯公式就可以算出后验概率。

所以说似然概率就是 （由中值定理可得：P(y<Y<y+ε|X=x)=f_Y|X( ξ₁|u) * ε）。

这个积分的意思：当真实值（X）取x的时候，你的观测值（Y）误差在正负0.2，即正负ε，之间，和似然概率密度f_y|x(y|x)围成的面积是1。红色的竖线代表f_y|x(y|x)，即似然概率密度。但要求似然概率密度还要通过假设似然概率的模型：等可能型，阶梯型和正态分布型（正态分布的好处：1.期望和方差独立，相对好控制（可以独立调整）。2. 中心极限定理：传感器的误差的来源相互独立的话，那么他们误差叠加的平均值就会满足正态分布）。

假设先验概率满足期望是10，方差为1的正态分布模型：

观测值为9，即u=9，精度为0.2，即σ=0.2，由于似然概率概率为

第三步，后验概率是通过先验概率与似然概率相乘，然后对后验概率做无穷积分，然后算出η ，最后的到后验概率。后验概率=(先验*似然概率)/f(y)= fx|y(x|y) = f(x,y)/f(y) = f_y|x(y|x)*f(x)/f(y) = η*f_y|x(y|x)*f(x) =1，其中f(y) = η（常数）。

由于归一化,所以：

最后得出后验概率为：

后验概率是通过先验概率与似然概率相乘，然后对后验概率做无穷积分，然后算出η ，最后的到后验概率

方差有明显的降低，不确定性减小，从而达到滤波的效果。

重要定理：若先验概率服从正态分布：f_x(x)~N(u₁,σ₁²)，且似然概率也服从正态分布：fy|x(y|x)~N(u₂,σ₂²)，那么后验概率就会服从正态分布。后验概率的重要定理中u₁与u₂前面的系数相加为1，说明他们的系数为权重，u₁和u₂哪个方差大，占的权重就小（可以通过暴力证明）。所以若σ₁²远大于σ₂²，那么后验概率就近似服从N(u₂,σ₂²)，反之亦然。

得到结论：σ₁²》σ₂² 时，后验概率倾向于观测值，而σ₁²《σ₁²时，后验概率倾向于预测值（先验概率）。并且，后验概率的方差中分别提取了σ₁²和σ₂²，新的滤波方差比原来的任何一个方差都要小，这带来一个很大的好处：尽可能削弱先验概率的影响，只要观测足够准确的话，你的后验概率的方差就可以非常的小，因为后验概率比先验概率和似然概率的方差任何一个都要小。哪怕预测和观测都不是很准，且似然概率的方差无论大/小，后验概率都可以通过滤波，都能得到相对准确，且比较小的方差。

方差越小，正态分布的区间就越短，当方差趋于0，曲线在期望上趋于无穷，而其他地方都是0，所以狄拉克函数是是一个陡峭的曲线。狄拉克函数 δ(x) 的定义是当x=0的时候δ(x)=oo，而x != 0的时候，δ(x)=o。而狄拉克函数的积分为1（归一化）：∫δ(x)dx=1。狄拉克函数最重要的特征是选择性：∫ f(x)*δ(x)dx =f(0)，看上去就像是用积分把f(0)选出来一样，因为f(x)只有在x=0的时候f(x)!=0，其他积分都等于0。

第六讲随机过程的贝叶斯滤波

1. 简介：随机过程的贝叶斯滤波。在随机过程中，有多个随机变量x，且不独立，有一个初始值x₀（靠猜），有k个观测值y₁,y₂,...y_k对应于每个变量。

2. 求解随机过程的解决方法：

1. 对所有的 x₀... x_k的先验概率都靠猜。
2. 只有x₀的概率式猜的，X₁....X_k的先验概率是递推的。

3. 贝叶斯滤波的步骤

随机过程的贝叶斯滤波的状态方程和观测方程

状态方程（预测方程）反映了X_k与X_k-1之间的关系， X_k=f(X_k-1) + Q，而Q是预测噪声。
观测方程就可以写成：Yk=h(Xk)+Rk，表示状态如何引起的传感器的读数Yk。

如何递推？贝叶斯滤波的基本思想：

首先设置初始值X₀ ~ N(0,1)，对这个状态方程做预测（预测步），得到X₁^-~N(0,2²)，若没有观测的话，直接利用X₁^-~ N(0,1)，去做预测，那么方差会越来越大，所以观测从这里开始引入（更新步），利用传感器所采集到的数据，对先验的数据进行更新（更新方法就是利用贝叶斯公式进行更新：后验=先验*似然（观测）），得到X₁⁺~ N(0,0.8)。再利用后验概率X₁⁺ 做为下一步得先验概率进行在预测，得到X₂^-~ N(0，3.2)，再根据观测数据Y₂，对X₂^-进行再更新，得到X₂⁺，以此类推。　　

X₀⁺~N(0,1) --预测步（通过状态方程X₁=2*X₀）--> X₁^--~N(0,2²) --更新步/后验步（观测方程：后验=先验*似然N(0,0.8)）-->X₁⁺~N(0,3.2) --再预测（状态方程）--> X₂^-~N(0,6.4)

递推分两步，第一步是预测步，是上一时刻的后验概率，通过状态方程，求出这一时刻的先验概率；第二步叫更新步，这一时刻的先验概率，通过观测方程以及观测数据，推出这一时段的后验概率，这一时刻的后验概率可以被当作下一时刻的先验概率。

4. 贝叶斯滤波算法：

贝叶斯滤波的前提（原料和条件）：

具体步骤：

1） 预测步（求先验概率）：

利用离散转连续的方式：P(X₁<x)=ΣP(X₁=u)，其中u是变量属于(-oo,x)。
而通过全概公式：P(X₁=u)=P(X₁=u|X₀=v)*P(X₀=v)。

由于X₀=v，所以可以条件概率的等式两边，同减去f(X₀)和f(v)。并将将状态方程带入，整理后，得到
已知Q1与X0是独立的，所以可以将X₀=v从条件概率中取出，得到P(Q1=u-f(v)|X₀=v)=P(Q₁=u-f(v)。
利用中值定理（P(Q₁=u-f(v)=f_Q1[u-f(v)]*ε，p(X₀=v)=f₀(v)*ε）
随机变量的概率公式：P=∫ f(x)*ε，其中ε 趋向于0。将概率化为概率密度乘以ε，然后再取极限。

所以P(X₁<x)的公式可得（上式漏掉一个ε，而下式中u属于(-oo,x)，而不是（-oo,+oo））：

2) 更新步（求似然概率）：

观测方程和更新步，主要是求似然概率

将观测方程带入，得P(y1-h(x)<Y1-h(x1)<y1-h(x)+ε)。

总结：递推公式推到过程。

f₀(x) --预测步（通过状态方程）--> f₁^-(x) --更新步/后验步（通过观测方程，且后验=先验*似然）--> f₁⁺(x) --预测步（通过状态方程）--> f₂^-(x)。

注：f1得后验概率与f2得后验概率是不相等的，所以他们的η也不一样。

预测步就是求先验概率，而更新步就是求似然概率，并求η。得到先验概率，似然概率和常数η，就可以得到后验概率。

贝叶斯滤波只是思想，虽然看上去算法很完美，但缺点就是无穷积分，在大多数情况下无法解析。

第七章贝叶斯滤波与卡尔曼滤波

1. 贝叶斯滤波：

1) 状态方程和观测方程，以及他们得前提假设：X₀，Q₁...Q_k，R₁...R_k相互独立

2) 贝叶斯滤波得三个公式：

2. 卡尔曼滤波

1）相对于贝叶斯滤波，多了几个假设：

f(x)与h(x)都是线性函数
f_Q(x)与f_R(x)都是服从正态分布。

2）预测步（求先验概率）：根据贝叶斯滤波，第一步是写出预测步公式，然后将相应得概率密度带入积分就行了。最后得到 u^-_k和σ^-_k。

3） 更新步（求似然概率）：已知f_k^-(x)服从正态分布，那么f_k⁺(x)后验概率就是常数 η 乘以似然概率 f_R(y_k-h*x)和先验概率f_k^-(x)。　　

利用数学公式去算，得到一个随机变量X_k⁺，但这个随机变量服从一个比较复杂的正态分布。通过变形得到u_k⁺和σ_k^-。

得到卡尔曼滤波的五个公式：预测步2个，更新步2个，外加卡尔曼增益。分别是代表从k-1的后验，通过两个独立的正态分布相加，推导出k的先验，再利用k的先验乘以似然概率（两个正态分布相乘再乘以常数η），推导到k的后验（其结果也是正态分布）。卡尔曼滤波之所以简单就是无需算积分，通过对期望和方差的计算，无需知道后验概率是什么形式。

第一讲绪论

1. 基本概率

HMM -- 隐马尔可夫模型

MCMC -- 可尔科夫-蒙特卡罗方法

马尔可夫链：未来与现在有关，而与过去无关。

MAP -- 最大后验估计

最小均方差估计

贝叶斯滤波

第二讲随机过程与主观概率

1. 基本概念

贝叶斯滤波的基本思想：通过概率统计的方法，利用贝叶斯公式，对随机信号进行处理，然后减小它的不确定性。

随机过程：x1,x2,x3....xn，其中x1,x2..xn..是一组随机变量，且彼此之间不独立的。（就无法做随机实验了。）

确定过程：v=gt（重力加速落体），通过重力加速度公式可以推出确定的过程：v1=g, v2=2g, v3=3g, .... vn

概率论和统计（大数定律、中心极限定理，数理统计（三大分布：卡方分布、T分布、F分布）等）都是以iid为前提的，但随机过程中的随机变量相互不独立，把独立性给去掉了，就无法做随机实验了。

2. 随机试验是什么？随机试验没了之后又什么影响？

随机试验是概率论的基础，它的最大的作用是个概率赋值的。比如，抛硬币，正面/反面朝上的概率是0.5。即怎么给一个概率事件赋一个数。赋值方法又两种：1.主观概率：概率是由主观定的，不是客观的东西，如：我自己主观认为抛硬币正面朝上的概率是0.5，2.大数定律：通过大量试验，得到正面朝上的次数与总次数之比，即正面朝上的频率，通过看次数增多的时候，看频率的趋势，抛硬币的正面朝上的频率在0.5附加波动，说明概率是0.5。第二种方法的基础就是随机试验。大数定律的例子如下（通过大量实验，当n（实验次数）趋向于+oo时，u/n（频率）依概率收敛于P1（0.5），所以说这个频率0.5是这么来的，是通过大量实验做出来的。）：

随机试验的定义，满足以下三个条件为随机实验：

　　1. 在相同条件下，试验可重复进行，即随机试验之间的独立性（如：抛硬币，第一次和第二次抛硬币之间没有任何影响，这样才能保证随机实验可重复进行，这样才可以通过多次实验来给概率赋值。若随机实验的独立性没有了，则无法再通过）

　　2. 一次试验，结果不确定，所有可能的结果已知，

　　3. 试验之前，试验结果预先未知。所以无随机试验，就无法对概率赋值。

随机过程：x1,x2,x3...xn 不独立，所以无法做随机试验，一般来说与时间相关无法做随机试验，如股票，分子扩散，气温变化等。这些随机过程是无法创造可以重复实验的环境出来。

随机过程：x1,x2,x3...xn.随机变量之间是不独立，如果我们找到了x_k与x_k-1之间关系，那么就可以找到p(x_k)与p(x_k-1)之间的关系。但还需要一个这个随机过程的起点，即x1。否则，你只是知道他们之间的关系，而没有初始值，仍然不能研究这个随机过程。即随机变量不独立，虽有了随机变量间的关系，但初值的选取很重要。

由于不独立性，所以没有办法通过大数定律给P(x1)一个初始值（无法定初始值），但这也得分情况看。

　　1. 有些随机过程，是可以主观给一个初始值，则是可以做随机试验的，比如说：随机游走：X_k=X_k-1+D，这个随机游走，有1/2的概率往前走，也有1/2的概率往后走。那么Xk与Xk-1代表着不同时间段，随机游走的位移。这是一个随机游走的过程，可以主观上设置初始值P(x0=0)=1。（独立性是不能通过逻辑来证明，P(AB)=PA*PB，公式中的PA和PB要通过概率赋值，这就牵扯到了随机实验，是否独立是做随机实验的关键。）

　　2. 有的初始值是不可以做随机实验的。那就只能使用主观概率，猜一个大致得概率。例如，我们找到了随机变量间得关系，但由于选取初始值是主观的，所以不同的主观概率导致不同的结果。

如分子扩散、气温变化等过程本质上是客观过程，不是主观的东西，我们希望尽可能削弱这个主观概率的差异，可以有多种方法解决这样的问题，如贝叶斯滤波。

3. 贝叶斯滤波

贝叶斯滤波的基本思想：引入外部的观测（证据、信息），来对主观的概率进行修正或更新。即主观概率（先验概率（先于实验之前的概率））通过外部观测（即贝叶斯公式），变为相对客观的概率，称这个概率叫做后验概率（实验之后的概率）

主观概率与先验概率的区别：先验概率可以通过看书或学习，大致得到一个概率，这种概率不是猜的，所以它不是主观概率，但是没有做实验，所以是先验概率。在此课程里，不区分主观概率与先验概率，视为等同。

第三讲贝叶斯概率的三大概率

1. 先验概率：没有做实验，通过学习等途径，得到的概率。包括主观概率（靠猜）

概率知识回顾，及用于规范：

3. 似然概率

引例：测温。首先，若无测量工具，就要猜一下今天的温度，即猜一个先验概率P(T=10)=0.8。但还要给出一个完整的先验概率分布，10度的概率是0.8，而11度的概率是0.2，即 ΣP(X=10) + ΣP(X=11) = ΣP(X=x) = 1（对于连续随机变量，其概率密度积分为1。如果给的是离散的先验概率，就要把完整的先验概率分布写出来），其次，利用设备测量，得到后验概率，但测量设备的精准度及测量误差等因素，会导致测量结果的偏差，从而观测结果不一定是真实结果，它也具有不确定性。所以才引入了第三个概率，似然概率。

首先，先验概率分布为，主观感觉温度为10度的概率P(T=10)为0.8，而主观感觉温度为11度的概率P(T=11)为0.2，先验概率分布的总和为1。

其次，观测到的测量结果T_m为10.3度，与先验概率一样，后验概率也是一个概率分布，所以它的概率分布总和为1，即P(T=10|T_m=10.3)+P(T=11|T_m=10.3)=1。但温度计有误差概率，所以引入了似然概率。

第三，似然概率代表着观测/传感器的精度，它衡量的是传感器的误差概率，叫似然概率。温度计的误差有可能是正负0.1，那么P(Tm=10.3|T=10)，真实温度为10的条件下，测量温度为10.3的概率就非常小，因为误差为正负0.1，那么真实值应该在10.2到10.4之间，所以Tm=10.3，而真实的值为10的概率很小。然而如果误差为正负0.5，那么P(Tm=10.3|T=10)就很大了。

似然概率与先验/后验概率不一样，它在不同的真实值情况下的测量结果的概率无需加总为1，因为P(Tm=10.3|T=10) 的本质意思是真实值是10的时候，测量温度为10.3的概率，而P(Tm=10.3|T=11)的本质意思是真实值是11的时候，测量温度为10.3的概率。所以说它是对两个真实值做的不同的测量，他们是两个不同的事件/随机试验，他们的概率没有关系。所以似然概率只是表示传感器的精度，只与传感器本身的误差有关，而与温度的真实值无关，不受其影响。

最后，随机变量的取值与随机变量的分布概率无关。这对P(Tm=10.3)是一个常数至关重要。

P(Tm=10.3) 代表观测值为10.3的概率。教程中，经常认为P(Tm=10.3)与T无关。但其实不然。

我们假设温度计测量的温度为10.3，那么温度为10.3 不是一个已发生的事件吗？所以P(Tm=10.3) 观测值为10.3的概率应该等于1。这样的疑惑是因为把随机变量的取值与随机变量的概率搞混了。如抛硬币，不能把一次正面朝上的结果看成是随机变量的概率，正面朝上只是随机试验的其中一个结果而已，不能把这一次的结果就把它当成是事件的概率，这一次的结果不能改变抛硬币的随机性的性质（也就是说随机事件：正面朝上和反面朝上的概率是0.5），不能把这个结果看成随机变量的概率，所以说随机变量的取值与随机变量的分布率没有关系。随机变量的取值代表一次随机试验的可能的结果，而随机变量的分布率代表你对随机变量的不确定性程度的刻画。

可以利用全概率公式来计算P（Tm=10.3）的概率：P(Tm=10.3) = P(Tm=10.3|T=10)P(T=10) + P(Tm=10.3|T=11)P(T=11)。公式的含义：Tm=10.3，它是结果，T=10，T=11，它是原因。这个全概公式的含义是：当真实值T=10或T=11的时候，测出的温度为Tm=10.3的概率之和，即P(Tm=10.3) 等于T=10的概率乘以传感器的精度+T=11时乘以传感器的精度。

P(Tm=10.3)是一个常数，是因为测试温度的概率P(Tm=10.3)与T 的取值无关，与T的分布律有关。 因为T的分布律时不会受实验的结果所影响，正如：抛硬币，无论是正面还是反面，它都不会影响它的分布率。T的分布率是先验概率，即实验前定的/猜的；而似然概率属于传感器的精度，它属于传感器本身固有的性质，所以说似然概率也是客观存在。如果先验概率定了，且似然概率也定了，那么P(Tm=10.3)的概率也就定了，是一个常数。

既然P(Tm=10.3)是一个常数（设这个常数为η），我们可以把这个贝叶斯公式进行改写：后验=η*似然*先验。

第四讲连续随机变量的贝叶斯公式

1. 似然概率得概念：

问：为什么叫把准确度叫做似然概率？
这是从最大似然估计那里借来的概念，似然的英语为likelihood，意思叫可能性，似然还有一个意思叫相似，相像，似然概率表示当我的观测值为10.3时，他最有可能，最像是哪个原因引起的，我们认为是因为有温度客观存在，才有对温度的观测，温度是因，观测是果，如果后验概率是说已知观测结果问原因的话，似然就是问哪个原因最像是引起观测是10.3这个结果，所以后验概率：T在条件概率前面，Tm在后面，表示已知结果问原因，而似然概率相反，T在后面，表示已知原因问结果。

从A/B班抽一个人出来，观测结果是一个女孩，这个女孩有可能从A班或B班出来，最大似然得意思就是这个女孩最大得可能是从B班出来。而后验概率就是知道这个女孩是从A/B班出来的，去寻找什么原因导致的（由果推因）。

2. 后验与似然的因果关系

一般把状态看成原因，观测看成结果，要先有状态，才能有观测，得先有温度，才可以观测。所以说状态是因，而观测是果。

似然概率是由因去推果（真实值->观测值的概率），真实的温度是10度，看导致观测值是10度的概率，观测值是11度的概率呢？

3. 独立未必没有函数关系

如果两个随机变量存在一定的函数关系，他们是不是一定不独立？
答：不一定。
等价命题：如果两个随机变量相互独立，他们是不是一定没有函数关系？
答：不一定。

独立，无关，没有影响。
独立，但又函数关系：1.必然事件，2. 随机事件中正态分布的 u 和 σ 就是独立，但有函数关系。而且只有在正态分布上有这个得关系，在其他分布上是没有这样得关系的。

证明：https://www.bilibili.com/video/av63950352/

4. 连续随机变量的贝叶斯

贝叶斯公式无法直接运用于连续随机变量，要将连续变量化积分为求和，即化积分的区域（-oo，x）为Δx点的累加求和。

连续随机变量X，Y的某个点的概率为0，即P(X=x)=0 或 P(Y=y)=0，严格来说不是0，算是无穷小，可以将无穷小写成极限的形式，这样就可以将离散的表达形式转化为积分形式：

P(y<Y<y+ε) = ∫ f_Y|X(y|u)dy -写成中值定理的形式-> f_Y|X( ξ₁|u) * ε，其中ξ₁属于（y,y+ε）。

中值定理：如果函数f(x)满足：（1）在闭区间[a,b]上连续；（2）在开区间(a,b)内可导；那么在开区间(a,b)内至少有一点使等式成立。上式转化中ξ₁，ξ₂，ξ₃采用了中值定理。

先验概率：f_x(x)

似然概率：f_Y|X(y|x)

后验概率：先验概率*似然概率 = f_x|y(x|y) = f_y|x(y|x) * f_x(x) / f_y(y）= f(x,y)/f(y)

f(x,y) 的边缘密度为：f_x(x) = ∫ f(x,y)dy, f_y(y) = ∫ f(x,y)dx。

第五讲似然概率与狄拉克函数

例子：测温度，先给一个先验的概率密度，数学期望是10，方差为1的标准正态分布。这样就意味着x=10的概率最大，方差是对你猜的不确定性的把握，方差就是误差的精度，这里的方差是1，所以容错的空间为[9,11]。本例中观测的温度为9度。如果有了先验概率，再得到似然概率，根据贝叶斯公式就可以算出后验概率。

y=x+-0.2 的意思就是观测值在0.2上下波动。利用概率的语言表述就是在X=x的情况下，Y在(x-0.2，x+0.2)区间的概率较大，而在这个区间之外的概率较小。即，在真实值 (X) 等于x的情况下，观测值(y)在0.2上下波动的概率很大。其实说的就是传感器的精度。0.9代表传感器的精度，表明还有10%的可能，误差超出0.2的范围。

这个积分的意思：当真实值（X）取x的时候，你的观测值（Y）在正负0.2之间，和似然概率密度f_y|x(y|x)围成的面积是1。红色的竖线代表f_y|x(y|x)，即似然概率密度。这个似然概率密度 f_y|x(y|x) 具体是什么值，一般来说是不知道的，一般来说传感器的精度是正负0.2或0.3，这样一个精度的范围，通过这个精度0.2/0.3，也只能推导似然概率密度的面积等于1，这一步，不能求出似然概率密度f_y|x(y|x)具体等于多少。 只能通过人为假设似然概率模型，再继续往下算。以下为似然概率的模型：　　

等可能型：它所有的概率都相等。当真实值X=x的情况下，Y在[x-0.2,x+0.2]的区间内取得每一个值都是相等的，而超出部分概率为0。说明Y服从均匀密度，即 f_y|x(y|x)=C 。∫fy|x(y|x) = 1 => C=2.5 ==》 |y-x|<=0.2时，似然概率f_y|x(y|x)=2.5，而当|y-x|>2的时候，似然概率f_y|x(y|x)=0。

阶梯型：它的概率密度也是常数，但它的概率密度是分段的常数，它的精度是分段的，越靠近真实值，它的概率就越大，远离真实值，它的概率就越小。阶梯型更靠近真实情况。

直方图型：这是阶梯型的推广，就是将非线性的概率密度，模拟成多个阶梯型的叠加。化连续为离散型。这样做的好处就是积分比较好积。
直方图滤波：由直方图型可以发展出一个直方图滤波（它是非线性卡尔曼滤波的一种），这个模型的优点在于便于积分，它可以把概率密度分成很多段的常数，常数是很好积分的，特别是概率密度的非线性比较复杂。

正态分布型：以真实值为期望，0.2为方差的正态分布曲线，但这样的话，Y在[x-0.2,x+0.2]之间的积分会小于1。一般来说，方差就取传感器的精度。正态分布好处：
- 正态分布的均值和方差比较好控制，其期望和方差没有什么关联（独立的）。而不像均匀分布，它的期望和方差就不好控制，它的期望定了，方差也就定了，它的E(x)=x，而D(x)=0.4²/12 != 2².
- 中心极限定理，只要误差是独立的，及泊松分布相互独立，有大量的误差累积，它的平均值就会满足这个正态分布。因为我们认为似然概率就是反映传感器的误差，那么传感器的误差的来源相互独立的话，那么他们误差叠加的平均值就会满足正态分布。

续前面正态分布的例子：

假设先验概率满足期望是10，方差为1的正态分布模型：

观测值为9，即u=9，精度为0.2，即σ=0.2，由于似然概率概率为

后验概率=(先验*似然概率)/f(y)= fx|y(x|y) = f(x,y)/f(y) = f_y|x(y|x)*f(x)/f(y) = =1，其中f(y) = η（常数）。

由于归一化,所以：

最后得出后验概率为：

总结：后验概率是通过先验概率与似然概率相乘，然后对后验概率做无穷积分，然后算出η ，最后的到后验概率。

方差有明显的降低，不确定性减小，从而达到滤波的效果。

后验概率的重要定理中u₁与u₂前面的系数相加为1，说明他们的系数为权重，u₁和u₂哪个方差大，占的权重就小（可以通过暴力证明）。所以若σ₁²远大于σ₂²，那么后验概率就近似服从N(u₂,σ₂²)，反之亦然。

2. 狄拉克函数

当σ趋于0的时候，这个正态分布的似然概率就会变成狄拉克函数。

δ(x)的本质：实质上为必然事件的概率密度。如P(X=0)=1，这是个离散的随机变量，是个必然事件，可以把这个离散的随机事件连续化得到P(X<x)。它的概率分布的图像是单位阶跃函数：当x=0的时候，突然抬升为1。而H(x)的导数就是狄拉克函数：δ(x)=H'(x)=dH(x)/dx。

狄拉克函数的性质:

1. x=0的时候，δ(x)=00。
2. δ(x)的在-oo到+oo的积分等于1：∫δ(x)=1.
3. δ(x)的选择性。

选择性的证明：

H(x)在-oo到0上是0。H(x)在0到+oo是1，则1-0=1.

推论：只要0属于[a,b]，那么就满足下列的积分，即选择性。

例子：

对这个后验概率进行积分，得出这个概率分布是一个必然事件，这也体现出这个传感器的精度是无限的，无论你的先验概率是怎么样的，我测到的是10，那么实际温度肯定就是10，这就是δ函数的应用。

第六讲随机过程贝叶斯滤波

前5讲，介绍了单个变量的贝叶斯滤波：通过随机变量x的先验概率和它的观测y，得到它的后验概率：。

而第六讲要介绍随机过程的贝叶斯滤波。在随机过程中，有多个随机变量x，且不独立，有一个初始值x₀(猜的)，有k个观测值y₁,y₂,...y_k对应于每个变量。

2. 求解随机过程的解决方法：

1. 对所有的 x₀... x_k的先验概率都靠猜。
- 缺点：过于依赖观测信息，放弃了预测信息取滤波，这样会导致无视随机过程的表达式，如这两个随机过程：x_k=2x_k-1+Q_k，x_k=x_k-1²+Q_k，这个方法会认为这两个表达式/随机过程是等同的。如果能写出X_k与X_k-1的关系式，哪怕是建模不准，它也是部分取反应随机过程的信息的，所以这个方法不是特别好。
2. 只有x₀的概率是猜的，而X₁....X_k的先验概率是递推的。
- 好处：保留了随机过程的特性，把X_k与X_k-1之间的关系给保留下来了。

3.贝叶斯滤波的基本思想　　

本文不涉及马尔可夫假设和观测独立假设，只利用状态方程和观测方程来推导贝叶斯滤波。

首先了解一下状态方程中，X_k与X_k-1是什么关系。例如：X_k=1/2*gt²+Q，Q可以认为是阻力或摩擦之类的随机扰动，可以认为是一个随机变量。我们可以利用一阶的泰勒展开，从而得到X_k与X_k-1的关系。（泰勒公式：f(x)=f(x0)/0! + f'(x0)/1! * (x-x0)...），我们这只展开了泰勒公式的第一阶，而没有展开更高阶，这样的问题是不大，是因为可以认为把建模的误差放入到Q中，然后把Q的方差调大一点就行。

如果可以能把建模建的更精确固然是好的，但大量的问题是无法写出精确的表达式，大多数情况下，都不知道曲线有什么规律。对于下面的曲线，我们用一个线性模型取建模是否是否合适？贝叶斯滤波的奥秘在于融合，融合先验和观测，先验是你没有做实验之前对这个问题的感性的认识，你主观认为的规律，观测是数据，数据可以给你带来信息，贝叶斯滤波就是这两个的融合。当你根本不知道这个问题是什么，或具体的形式，正如下面的曲线，无法预测随机过程或表达式，那么就贝叶斯滤波只能通过观测来帮你估一个大概的东西。只要写出一个大概的表达式出来，就有一些基本信息，贝叶斯滤波就有用武之地。

注：状态方程（预测方程）反映了X_k与X_k-1之间的关系， X_k=f(X_k-1) + Q，而Q是预测噪声。

观测方程反映了状态时如何引起传感器的读数

状态一般是你想要研究的对象，如速度，温度，角度等，状态方程一般不是我们想要的状态数据，往往是一些跟状态有关的函数数据。例如，测温度，你想测的是温度，传感器的读数也是温度，那么观测方程就可以写成：Yk=Xk+Rk，意思是状态是温度，温度X_k加上了长期误差R_k 就变成了传感器的读数Y_k。

但大多数情况下，状态跟观测往往不一致，比如说状态时位移，观测是角度的传感器，那么状态方程为：Yk=arcsinXk+Rk，也就是说位移和角度一一对应，将位移转化成角度，再将角度加上一个传感器的噪声误差，才会生成传感器的读数Yk。

那么观测方程就可以写成：Yk=h(Xk)+Rk，表示状态如何引起的传感器的读数Yk。

这样就可以把上面的状态方程（随机方程）和观测方程联立起来，表示这个随机过程怎么引起我的观测：

若只通过先验概率，没有观测，是无法滤波的。如下面的例子，无观测方程，状态方程的方差越来越大，最终无法滤波。

那么我们得到联立方程，怎么引入观测来递推呢？

首先设置初始值X₀⁺~ N(0,1)，对这个状态方程做预测（预测步），得到X₁^-~N(0,2²)，若没有观测的话，直接利用X₁^-~ N(0,1)，去做预测，那么方差会越来越大，所以观测从这里开始引入（更新步/后验步），利用传感器所采集到的数据y1=0，对先验的数据进行更新（更新方法就是利用贝叶斯公式进行更新：后验=先验*似然（观测）），得到X₁⁺~ N(0,0.8)。再利用后验概率X₁⁺ 做为下一步的先验概率进行再预测，得到X₂^-~ N(0，3.2)，再根据观测数据Y₂，对X₂^-进行再更新，得到X₂⁺，以此类推。　　

4. 贝叶斯滤波的算法

贝叶斯滤波算法所需要的原料和条件：

1. 首先是状态方程和观测方程. 假设：X₀，Q₁...Q_k，R₁....R_k之间是相互独立的。{Q_k}是预测噪声，是建模的误差，{R_k}是观测噪声，是传感器的误差。X₀是主观猜的初始值。
3. 有一系列的观测值y1.y2...yk
4. 设初始值X₀的概率密度pdf f₀(x)，Q_k的概率密度pdf：f_Qk(x)，R_k的概率密度f_Rk(x)。
重要定理：条件概率例的条件可以做逻辑推导。条件概率的转化，不能省略条件。

具体步骤：

1）. 预测步(求先验概率)：

P(X₁<x)=ΣP(X₁=u)，其中u是变量属于(-oo,x)，这是离散转连续的方式。而通过全概公式：P(X₁=u)=P(X₁=u|X₀=v)*P(X₀=v)。

由于条件概率中的条件是X₀=v，所以可以条件概率的等式两边，同减去f(X₀)和f(v)。然后将原料中的状态方程带入，整理后，得到，其中从原料和条件中，可知Q1与X0是独立的，所以可以将X₀=v从条件概率中取出，得到P(Q1=u-f(v)|X₀=v)=P(Q₁=u-f(v)。并要证明Qk与Xk-1是相互独立的（证明后置）。若他们相互独立，从X0到X1，X1到X2....，Xk-1到Xk的递推逻辑才能成立。

接下来就是求上述的连续随机变量的概率（随机变量的概率公式：P=∫ f(x)*ε，其中ε 趋向于0。将概率化为概率密度乘以ε，然后再取极限。 ），并利用中值定理（P(Q₁=u-f(v)=f_Q1[u-f(v)]*ε，p(X₀=v)=f₀(v)*ε）就得到了下面的公式：

所以P(X₁<x)的公式可得（下式漏掉一个ε，并且u属于(-oo,x)）：

2）. 更新步（求似然概率）：

观测方程和更新步，主要是求似然概率

将观测方程带入，得P(y1-h(x)<Y1-h(x1)<y1-h(x)+ε)。

总结：递推公式推到过程。

f1得后验概率与f2得后验概率是不相等的，所以他们的η也不一样。

最后证明一下Q_k与X_k-1的独立性，X_k与R_k独立

总结：预测步就是求先验概率，而更新步就是求似然概率，并求η。得到先验概率，似然概率和常数η，就可以得到后验概率。

贝叶斯滤波只是思想，虽然看上去算法很完美，但缺点就是无穷积分，在大多数情况下无法解析。

三大无穷积分：

解决方法：

卡尔曼率波就是对Q_k，R_k做正态分布假设。基本思想：由上一时刻的期望和方差，推出这一时刻的期望和方差的先验，用这一时刻的方差的先验再推出这一时刻期望和方差的后验，再加上一个卡尔曼增益。

蒙特卡罗积分就是粒子滤波

直方图滤波：把曲线分成多段的常数，看起来就像直方图一样，这样就好积分了。

第七讲贝叶斯滤波与卡尔曼滤波

1. 贝叶斯滤波：

1) 状态方程和观测方程，以及他们得前提假设：X₀，Q₁...Q_k，R₁...R_k相互独立

2) 贝叶斯滤波得三个公式：

2. 卡尔曼滤波

卡尔曼滤波相对于贝叶斯滤波，多了几个假设：

f(x)与h(x)都是线性函数
f_Q(x)与f_R(x)都是服从正态分布。

1）首先假设X_k-1服从期望为u_k-1，方差为σ_k-1，得正态分布。

2）其次预测步：根据贝叶斯滤波，第一步是写出预测步公式，然后将相应得概率密度带入积分就行了。最后得到 u^-_k和σ^-_k。

傅里叶变换：

前提假设：假设Xk-1服从正态分布， FXk-1也服从正态分布，Qk服从正态分布(0,Q)，Xk-1与Qk之间是相互独立的。

Xk=FXk+Qk，实际上是两个相互独立的的随机变量相加，它的概率密度实际上是这两个独立的随机变量的卷积。卷积有一个性质是：如果h是f与g的卷积的话，对h做傅里叶变换，就等于对f和g做傅里叶变换的乘积。那么我们只需要把FXk-1做傅里叶变换，再把Qk做傅里叶变换，然后把他们乘起来，再做傅里叶逆变换，就会得到Xk的概率密度。虚部 iut 是期望，实部的分子σ²是方差，相对于正态分布的F_uk-1，F²_σk-1+Q。

3）更新步：已知f_k^-(x)服从正态分布，那么f_k⁺(x)后验概率就是常数 η 乘以似然概率 f_R(y_k-h*x)和先验概率f_k^-(x)。　　

利用数学公式去算，得到一个随机变量X_k⁺，但这个随机变量服从一个比较复杂的正态分布。通过变形得到u_k⁺和σ_k^-。

这样就可以得到了后验概率的期望和方差。它们有相同的因子：h*σ_k^-/（h₂σ_k^-+R），这个因子K叫做卡尔曼增益。

这样我们我们就可以对上面的后验概率的期望和方差进行改写：

这样我们就可以得到卡尔曼滤波的五个公式：预测步2个，更新步2个，外加卡尔曼增益。分别是代表从k-1的后验，通过两个独立的正态分布相加，推导出k的先验，再利用k的先验乘以似然概率（两个正态分布相乘再乘以常数η），推导到k的后验（其结果也是正态分布）。卡尔曼滤波之所以简单就是无需算积分，通过对期望和方差的计算，无需知道后验概率是什么形式。

可以发现卡尔曼公式没有后验概率密度，没有积分之类的东西，而只是一些普通的数，进行加加减减。所以具有非常好的性质（这是由正态分布非常好的性质决定的）。正态分布具有封闭性，两个独立的正态分布，通过卷积，相乘，归一化，最终还是正态分布。但对于非线性就没有这么好的性质了，所以非线性的卡尔曼滤波就比较难，它的计算量比较大，因为它已经不再是正态分布，所以你必须计算出后验的概率分布出来。

我们对卡尔曼增益做变形，将σ_k^-除下去。

当R》σ_k^-时，k就趋于0，uk+ ~= uk-，即相信预测结果。这是合理的，因为R是观测的方差，σ_k^-是预测的方差。当R远大于σ_k^-的时候，就说明观测不可信（噪声太大），那么卡尔曼滤波就会偏向预测值。
同理，R《σ_k^-时，我们更相信观测值。

矩阵形式的卡尔曼滤波：与一维卡尔曼滤波很像，其中期望与方差，及F，H都会变为矩阵形式，

其公式也很像，但涉及多元的正态分布概率密度的形式，以及多重的无穷积分，比较复杂。

第八讲从零开始码出卡尔曼滤波

应用：more is different　　

模型1：

# loding modules
import numpy as np
import pandas as pd
from pandas import *
import matplotlib.pyplot as plt

#X(k) = F*X(K-1)+Q
#Y(k) = H*X(K)+R
#第一个问题，生成一段随机信号，并滤波

#生成一段时间t，每隔0.01秒进行采样一次。
t = np.arange(0.1,1,0.01)
L = len(t)

#有了时间，就可以生成真实信号x，以及观测y。
#首先初始化
x=pd.Series(np.zeros(L))
y=pd.Series(np.zeros(L))    #不能将x赋值给y（即y=x），这样只会将x与y用做同一内存空间。
#生成信号，设x=t^2
for i in range(L):
    x[i]=t[i]**2
    y[i]=x[i]+np.random.normal(0,0.1)

plt.plot(t,x)
plt.plot(t,y)

#滤波算法
#预测方程怎么写：
#观测方程：Y(K)=X(K)+R R~N(0,1)，贝叶斯滤波也卡尔曼滤波的难点就是Q和R的方差的取值，它的难点在于没有规律性，只能不停的去试。
#预测方程不好写，画图（plt.plot(t,y)）后，可以看出曲线的一些规律性，可以猜一猜是线性增长，但是大多数问题，信号是杂乱无章的，怎么办。
#模型一，最粗糙的建模：
#X(k)=X(k-1)+Q
#Y(k)=X(k)+R
#猜Q~N(0,1)，弄成正态分布，相对好处理一些。
F1=1
H1=1
Q1=1
R1=1
#初始化X(k)+，要设置一个数组来装这个X(k)+，因为卡尔曼滤波是递推的，X(k)+需要X(k-1)+。
Xplus1=Series(np.zeros(L))
Pplus1=Series(np.zeros(L))
Xminus1=Series(np.zeros(L))
Pminus1=Series(np.zeros(L))
K1=Series(np.zeros(L))

#设置一个初值，假设Xplus1(1)~N(0.01,0.01^2)
Xplus1[0]=0.01
Pplus1[0]= 0.01**2

#卡尔曼滤波算法：它是期望与方差之间的递推关系，用前一时刻plus的期望和方差，推出这一时刻minus的期望和方差。以及这一时刻minus的期望和方差，推出这一时刻plus的期望和方差。再加上一个卡尔曼增益。X(0)是没有minus，且带plus的先验概率的，除了X(0)外所有X(k)的minus为先验概率，而plus为后验概率。
#X(k)minus = F*X(k-1)plus
#P(k)minus = F*P(k-1)plus*F+Q
#K = P(k)minus*H*inv(H*P(k)minus*H'+R)
#X(k)plus = X(k)minus+K*(y(k)-H*X(K)minus)
#P(k)plus = (1-K*H)*P(k)minus
for i in range(1,L):
    ##预测步：
    Xminus1[i]=F1*Xplus1[i-1]
    Pminus1[i]=F1*Pplus1[i-1]*F1'+Q1        #F1'是Matlab求转置。
    ##更新步：
    K1[i] = (Pminus1[i]*H1')/(H1*Pminus1[i]*H1'+R1)
    Xplus1[i] = Xminus1[i]+K1[i]*(y[i]-Xminus1[i]*H1)
    Pplus1[i]=(1-K1[i]*H1)*Pminus1[i]

plt.plot(figsize=(20,12))
plt.plot(t,x,'r')
plt.plot(t,y,'g')
plt.plot(t,Xplus1,'b')

#Q1,R1=0.01,控制噪声的方差。
#Xplus1[0]控制初始值，初始值的大小对滤波影响不大，最后滤波会调整。

模型1

模型2：

###### 模型2 ##########
#模型2相对于模型1相比，观测方程还是一样，但预测方程一项，把X(k)写成用X(k-1)展开的一些项（泰勒展开式）。
#X(K)=X(K-1)*dt+X'(K-1)*dt+X''(K-1)*dt^2*(1/2!)
#Y(K)=X(K)+R R~N(0,1)
#此时状态变量不再是标量了，而是向量：X=[X(k) X'(K) X''(K)]T 列向量
#Y(K)=H*X+R H=[1 0 0] 行向量

#预测方程
#X(K)=X(K-1)+X'(K-1)*dt+X''(K-1)*dt^2*(1/2!)+Q2
#X'(K)=0*X(K-1)+X'(K-1)+X''(K-1)*dt+Q3
#X''(K)=0*X(K-1)+0X'(K-1)+X''(K-1)*dt+Q3 多项式信号多求多阶导数，总会比较平缓，而
#X'''(K)=X''(k-1)+Q3正是描述平衡的随机过程，这种建模相对精细一些，适用范围也较广。它的原理就是无论是什么函数，多求几次导，总会求到0，0就比较平缓了。对导数建随机平稳得模型，要比直接对多项式函数建平稳得模型相对来说精细一点。但有可能信号是指数，那么怎么建都不平缓。一般来说，信号不会增长太快，一般物理规律，信号不会增长太快。
#F= 1 dt 0.5*dt^2    #F是一个矩阵，在计算中要注意行向量和列向量。
#   0  1    dt
#   0  0    1
#H=[1 0 0]
#Q = Q2  0  0        #Q是一个噪声的协方差矩阵，Q协方差是用来描述他们随机变量的噪声的相互关系。由于状态变量X是由3个元素组成，Xk,Xk',Xk'',他们是多维随机变量，有各自的方差。而他们之间的关系，若两个随机变量的相关性为0，那么他们是不相关，若Xk增大了，Xk’是否会有变化，这需要用协方差来表述。
#    0   Q3 0
#    0   0  Q4
dt = t[2]-t[1]
F2 = np.mat([[1,dt,0.5*dt**2],[0,1,dt],[0,0,1]])    #此处要注意矩阵是否病态，若dt特别小，易导致矩阵病态或精度丢失
H2 = np.mat([[1,0,0]])
Q2 = np.mat([[1,0,0],[0,0.01,0],[0,0,0.001]])
R2 = 20            #只有一个观测，所以只有是标量，如果有多个观测，得写成矩阵形式。

#设置初始值
Xplus2 = np.matrix(np.zeros((3,L)))
Xplus2[0,0]=0.1**2
Xplus2[1,0]=0
Xplus2[2,0]=0
Pplus2=np.mat([[0.01,0,0],[0,0.01,0],[0,0,0.0001]])     #一般来说，初始得协方差矩阵都写成是对角矩阵形式。
#K2 = np.mat(np.zeros((3,1)))
for i in range(1,L):
    #预测步
    Xminus2[:,i]=F2*Xplus2[:,i-1]
    Pminus2=F2*Pplus2*F2.T+Q2
    #更新步
    K2=(Pminus2*H2.T)*np.linalg.inv(H2*Pminus2*H2.T+R2)
    Xplus2[:,i]=Xminus2[:,i] +K2*(y[i]-H2*Xminus2[:,i])
    Pplus2=(np.identity(3)-K2*H2)*Pminus2

Xplus2_0=np.zeros(L)
for i in range(L):
    Xplus2_0[i]=Xplus2[0,i]
plt.plot(figsize=(20,12))
plt.plot(t,x,'r')
plt.plot(t,y,'g')
plt.plot(t,Xplus1,'b')
plt.plot(t,Xplus2_0,'m')

# 可以进行在线滤波，实时滤波


### 问题2，两个传感器，进行滤波 ###
# Y1(K) = X(K)+R
# Y2(K) = X(K)+R
# H=[1 1]T （列向量） X=X(K)
# H=[[1 0 0],
         [1 0 0]]

模型2

第九到第十二讲：粒子滤波

若f(Xk-1) h(Xk)都是非线性函数，要用粒子滤波，这是处理非线性最准的方法，但是粒子滤波的缺陷是对内存消耗太大。计算非常的慢，因为需要重采样，如果采的不是正态分布的样本，那么会非常慢。粒子滤波的精度是随着粒子增多而提高的。

第十三讲扩展卡尔曼滤波（EKF）

扩展的卡尔曼滤波处理非线性函数的速度比粒子滤波快，但是这个快是要牺牲精度的。特别是在强非线性的问题中，这种方法很有可能会发散。

扩展卡尔曼滤波的思想是对f(X_k-1)和h(X_k)进行线性化。

扩展卡尔曼滤波的工作原理：

设X_k-1服从期望为X_K-1⁺,方差为P_k-1的正态分布，其实X_k-1⁺和P_k-1⁺是一个数，而不是一个随机变量，这是扩展卡尔曼滤波的源头。 X_k-1是一个正态分布，通过递推，扩展卡尔曼滤波所算出来的X_K⁺也得是正态分布，只有这样算法的递推才能进行下去。因为X_k-1是正态分布，推出X_k，又由X_k去推X_k+1，所以必须保证X_k也是正态分布才可以，否则就没有办法递推了。如果X_k不服从正态分布，就没有办法调用卡尔曼滤波来算了。这么一来，{X_k}序列都要服从正态分布。这样就会引起误差，因为F和H本质上是非线性函数。所以说真实的后验概率它一定不是正态分布，因为它是非线性函数，但滤波的结果是正态分布，误差就在这体现。所以扩展的卡尔曼滤波是一种近似，这种近似是有误差的。

那么扩展卡尔曼滤波是怎么处理的：

1）预测步

首先对预测方程进行线性化，即用泰勒公式展开（只展开一阶），得f(X_k-1)~= f'(x_k-1⁺)+f(x_k-1⁺)-f(x_k-1⁺)*x_k-1⁺。其中f(x_k-1)，f'(x_k-1⁺)，x_k-1⁺他们都是数，所以就线性化了。因为他们是数，所以可以将X前的系数改为A=f'(x_k-1⁺)，B=f(x_k-1⁺)-f(x_k-1)*x_k-1⁺。这样就可以得到f(X_k-1) = AX_k-1 + B，其中系数A和B都是数。这样就完成了预测方程的线性化。

然后将预测方程带入到贝叶斯滤波里，贝叶斯滤波中的预测步为：，将上面得到的 f(v) = Av+B 带入到贝叶斯滤波的预测步中，得到。积分后得到了（利用软件积分得到）期望为Ax_k-1⁺+B，方差为A²P_k-1⁺+Q的正态分布。严格来说f_k^-(x)是正态分布的概率密度，而不是正态分布。
将其进一步简化：