原理

假设有一个系统,服从一阶马尔科夫模型,我们知道它的状态方程,和观测方程如下:
x k = f ( x k − 1 , Q k ) y k = h ( x k , R k ) (1) x_k=f(x_{k-1},Q_k) \tag{1}\\ y_k=h(x_k,R_k) xk=f(xk1,Qk)yk=h(xk,Rk)(1)
设已知 k − 1 k-1 k1时刻的概率密度函数 p ( x k − 1 ∣ y 1 : k − 1 ) p(x_{k-1}|y_{1:k-1}) p(xk1y1:k1)
1,预测,根据状态方程求 f : x k − 1 → x k f:x_{k-1}\to x_k f:xk1xk
p ( x k ∣ y 1 : k − 1 ) = ∫ p ( x k , x k − 1 ∣ y 1 : k − 1 ) d x k − 1 = ∫ p ( x k ∣ x k − 1 , y 1 : k − 1 ) p ( x k − 1 ∣ y 1 : k − 1 ) d x k − 1 = ∫ p ( x k ∣ x k − 1 ) p ( x k − 1 ∣ y 1 : k − 1 ) d x k − 1 \begin{aligned} p(x_k|y_{1:k-1}) &= \int p(x_k,x_{k-1}|y_{1:k-1})dx_{k-1} \\ &= \int p(x_k|x_{k-1},y_{1:k-1})p(x_{k-1}|y_{1:k-1})dx_{k-1} \\ &= \int p(x_k|x_{k-1})p(x_{k-1}|y_{1:k-1})dx_{k-1} \end{aligned} p(xky1:k1)=p(xk,xk1y1:k1)dxk1=p(xkxk1,y1:k1)p(xk1y1:k1)dxk1=p(xkxk1)p(xk1y1:k1)dxk1
采样过程噪声 Q k Q_k Qk再叠加 f ( x k − 1 ) f(x_{k-1}) f(xk1)就可以得到 x k x_k xk的分布;
2,更新,根据观测方程求 h : x k → y k h:x_k\to y_k h:xkyk
p ( x k ∣ y 1 : k ) = p ( y k ∣ x k , y 1 : k − 1 ) p ( x k ∣ y 1 : k − 1 ) p ( y k ∣ y 1 : k − 1 ) = p ( y k ∣ x k ) p ( x k ∣ y 1 : k − 1 ) p ( y k ∣ y 1 : k − 1 ) = p ( y k ∣ x k ) p ( x k ∣ y 1 : k − 1 ) ∫ p ( y k ∣ x k ) p ( x k ∣ y 1 : k − 1 ) d x k \begin{aligned} p(x_k|y_{1:k})&=\frac {p(y_k|x_k,y_{1:k-1})p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})} \\ &=\frac {p(y_k|x_k)p(x_k|y_{1:k-1})}{p(y_k|y_{1:k-1})} \\ &=\frac {p(y_k|x_k)p(x_k|y_{1:k-1})}{\int p(y_k|x_k)p(x_k|y_{1:k-1})dx_k} \\ \end{aligned} p(xky1:k)=p(yky1:k1)p(ykxk,y1:k1)p(xky1:k1)=p(yky1:k1)p(ykxk)p(xky1:k1)=p(ykxk)p(xky1:k1)dxkp(ykxk)p(xky1:k1)
p ( y k ∣ x k ) p(y_k|x_k) p(ykxk)与观测噪声 R k R_k Rk的概率分布有关,利用蒙特卡洛采样计算过上式中的分母积分;

蒙特卡洛采样

蒙特卡洛采样解决概率积分难的问题:从概率分布 p ( x ) p(x) p(x)中采样到样本 x 1 , x 2 , ⋯   , x N x_1,x_2,\cdots,x_N x1,x2,,xN,利用这些样本去估计这个分布的某些函数的期望与方差,如下,
E [ f ( x ) ] = ∫ f ( x ) p ( x ) d x ≈ f ( x 1 ) + f ( x 2 ) + ⋯ + f ( x N ) N E[f(x)]=\int f(x)p(x)dx\approx \frac{f(x_1)+f(x_2)+\cdots+f(x_N)}{N} E[f(x)]=f(x)p(x)dxNf(x1)+f(x2)++f(xN)

重要性采样

重要性采样解决分布难预知的问题(分布转换技巧)
E [ f ( x k ) ] = ∫ f ( x k ) p ( x k ∣ y 1 : k ) d x k = ∫ f ( x k ) p ( x k ∣ y 1 : k ) q ( x k ∣ y 1 : k ) q ( x k ∣ y 1 : k ) d x k = ∫ f ( x k ) p ( y 1 : k ∣ x k ) p ( x k ) p ( y 1 : k ) q ( x k ∣ y 1 : k ) q ( x k ∣ y 1 : k ) d x k = ∫ f ( x k ) W k ( x k ) p ( y 1 : k ) q ( x k ∣ y 1 : k ) d x k = 1 p ( y 1 : k ) ∫ f ( x k ) W k ( x k ) q ( x k ∣ y 1 : k ) d x k = ∫ f ( x k ) W k ( x k ) q ( x k ∣ y 1 : k ) d x k ∫ p ( y 1 : k ∣ x k ) p ( x k ) d x k = ∫ f ( x k ) W k ( x k ) q ( x k ∣ y 1 : k ) d x k ∫ W k ( x k ) q ( x k ∣ y 1 : k ) d x k = E q ( x k ∣ y 1 : k ) [ W k ( x k ) f ( x k ) ] E q ( x k ∣ y 1 : k ) [ W k ( x k ) ] ≈ 1 N ∑ i = 1 N W k ( x k ( i ) ) f ( x k ( i ) ) 1 N ∑ i = 1 N W k ( x k ( i ) ) = ∑ i = 1 N W k ( x k ( i ) ) ∑ i = 1 N W k ( x k ( i ) ) f ( x k ( i ) ) = ∑ i = 1 N W ~ k ( x k ( i ) ) f ( x k ( i ) ) \begin{aligned} E[f(x_k)]&=\int f(x_k)p(x_k|y_{1:k})dx_k \\ &=\int f(x_k)\frac{p(x_k|y_{1:k})}{q(x_k|y_{1:k})}q(x_k|y_{1:k})dx_k \\ &=\int f(x_k)\frac{p(y_{1:k}|x_k)p(x_k)}{p(y_{1:k})q(x_k|y_{1:k})}q(x_k|y_{1:k})dx_k \\ &=\int f(x_k)\frac{W_k(x_k)}{p(y_{1:k})}q(x_k|y_{1:k})dx_k \\ &=\frac{1}{p(y_{1:k})}\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k \\ &=\frac{\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k}{\int p(y_{1:k}|x_k)p(x_k)dx_k} \\ &=\frac{\int f(x_k)W_k(x_k)q(x_k|y_{1:k})dx_k}{\int W_k(x_k)q(x_k|y_{1:k})dx_k} \\ &=\frac{E_{q(x_k|y_{1:k})}[W_k(x_k)f(x_k)]}{E_{q(x_k|y_{1:k})}[W_k(x_k)]} \\ &\approx \frac{\frac{1}{N}\sum_{i=1}^{N} W_k(x_k^{(i)})f(x_k^{(i)})}{\frac{1}{N}\sum_{i=1}^{N} W_k(x_k^{(i)})} \\ &=\sum_{i=1}^{N} \frac{W_k(x_k^{(i)})}{\sum_{i=1}^{N} W_k(x_k^{(i)})}f(x_k^{(i)}) \\ &=\sum_{i=1}^{N} \tilde{W}_k(x_k^{(i)})f(x_k^{(i)}) \end{aligned} E[f(xk)]=f(xk)p(xky1:k)dxk=f(xk)q(xky1:k)p(xky1:k)q(xky1:k)dxk=f(xk)p(y1:k)q(xky1:k)p(y1:kxk)p(xk)q(xky1:k)dxk=f(xk)p(y1:k)Wk(xk)q(xky1:k)dxk=p(y1:k)1f(xk)Wk(xk)q(xky1:k)dxk=p(y1:kxk)p(xk)dxkf(xk)Wk(xk)q(xky1:k)dxk=Wk(xk)q(xky1:k)dxkf(xk)Wk(xk)q(xky1:k)dxk=Eq(xky1:k)[Wk(xk)]Eq(xky1:k)[Wk(xk)f(xk)]N1i=1NWk(xk(i))N1i=1NWk(xk(i))f(xk(i))=i=1Ni=1NWk(xk(i))Wk(xk(i))f(xk(i))=i=1NW~k(xk(i))f(xk(i))
权重的递推公式 W k ( i ) W_k^{(i)} Wk(i)
W k ( i ) = W k ( x k ( i ) ) = p ( y 1 : k ∣ x k ) p ( x k ) q ( x k ∣ y 1 : k ) ∝ p ( x k ∣ y 1 : k ) q ( x k ∣ y 1 : k ) ≡ p ( x 0 : k ∣ y 1 : k ) q ( x 0 : k ∣ y 1 : k ) = p ( y k ∣ x 0 : k , y 1 : k − 1 ) p ( x 0 : k ∣ y 1 : k − 1 ) / p ( y k ∣ y 1 : k − 1 ) q ( x 0 : k − 1 ∣ y 1 : k − 1 ) q ( x k ∣ x 0 : k − 1 , y 1 : k ) = p ( y k ∣ x 0 : k , y 1 : k − 1 ) p ( x k ∣ x 0 : k − 1 , y 1 : k − 1 ) p ( x 0 : k − 1 ∣ y 1 : k − 1 ) / p ( y k ∣ y 1 : k − 1 ) q ( x 0 : k − 1 ∣ y 1 : k − 1 ) q ( x k ∣ x 0 : k − 1 , y 1 : k ) = p ( y k ∣ x k ) p ( x k ∣ x k − 1 ) p ( x 0 : k − 1 ∣ y 1 : k − 1 ) / p ( y k ∣ y 1 : k − 1 ) q ( x 0 : k − 1 ∣ y 1 : k − 1 ) q ( x k ∣ x 0 : k − 1 , y 1 : k ) ∝ p ( y k ∣ x k ) p ( x k ∣ x k − 1 ) p ( x 0 : k − 1 ∣ y 1 : k − 1 ) q ( x 0 : k − 1 ∣ y 1 : k − 1 ) q ( x k ∣ x 0 : k − 1 , y 1 : k ) = W k − 1 ( i ) p ( y k ∣ x k ( i ) ) p ( x k ( i ) ∣ x k − 1 ( i ) ) q ( x k ( i ) ∣ x 0 : k − 1 ( i ) , y 1 : k ) \begin{aligned} W_k^{(i)}&=W_k(x_k^{(i)}) \\ &=\frac{p(y_{1:k}|x_k)p(x_k)}{q(x_k|y_{1:k})} \\ &\propto \frac{p(x_k|y_{1:k})}{q(x_k|y_{1:k})} \\ &\equiv \frac{p(x_{0:k}|y_{1:k})}{q(x_{0:k}|y_{1:k})} \\ &=\frac{p(y_k|x_{0:k},y_{1:k-1})p(x_{0:k}|y_{1:k-1})/p(y_k|y_{1:k-1})}{q(x_{0:k-1}|y_{1:k-1})q(x_k|x_{0:k-1},y_{1:k})} \\ &=\frac{p(y_k|x_{0:k},y_{1:k-1})p(x_k|x_{0:k-1},y_{1:k-1})p(x_{0:k-1}|y_{1:k-1})/p(y_k|y_{1:k-1})}{q(x_{0:k-1}|y_{1:k-1})q(x_k|x_{0:k-1},y_{1:k})} \\ &=\frac{p(y_k|x_k)p(x_k|x_{k-1})p(x_{0:k-1}|y_{1:k-1})/p(y_k|y_{1:k-1})}{q(x_{0:k-1}|y_{1:k-1})q(x_k|x_{0:k-1},y_{1:k})} \\ &\propto \frac{p(y_k|x_k)p(x_k|x_{k-1})p(x_{0:k-1}|y_{1:k-1})}{q(x_{0:k-1}|y_{1:k-1})q(x_k|x_{0:k-1},y_{1:k})} \\ &=W_{k-1}^{(i)}\frac{p(y_k|x_k^{(i)})p(x_k^{(i)}|x_{k-1}^{(i)})}{q(x_k^{(i)}|x_{0:k-1}^{(i)},y_{1:k})} \\ \end{aligned} Wk(i)=Wk(xk(i))=q(xky1:k)p(y1:kxk)p(xk)q(xky1:k)p(xky1:k)q(x0:ky1:k)p(x0:ky1:k)=q(x0:k1y1:k1)q(xkx0:k1,y1:k)p(ykx0:k,y1:k1)p(x0:ky1:k1)/p(yky1:k1)=q(x0:k1y1:k1)q(xkx0:k1,y1:k)p(ykx0:k,y1:k1)p(xkx0:k1,y1:k1)p(x0:k1y1:k1)/p(yky1:k1)=q(x0:k1y1:k1)q(xkx0:k1,y1:k)p(ykxk)p(xkxk1)p(x0:k1y1:k1)/p(yky1:k1)q(x0:k1y1:k1)q(xkx0:k1,y1:k)p(ykxk)p(xkxk1)p(x0:k1y1:k1)=Wk1(i)q(xk(i)x0:k1(i),y1:k)p(ykxk(i))p(xk(i)xk1(i))
前提假设重要性分布满足 q ( x k ∣ x 0 : k − 1 , y 1 : k ) = q ( x k ∣ x k − 1 , y k ) q(x_k|x_{0:k-1},y_{1:k})=q(x_k|x_{k-1},y_k) q(xkx0:k1,y1:k)=q(xkxk1,yk)

伪代码
[ { x k ( i ) , w k ( i ) } i = 1 N ] = S I S [ { x k − 1 ( i ) , w k − 1 ( i ) } i = 1 N , y 1 : k ] [\{x_k^{(i)},w_k^{(i)}\}_{i=1}^N]=SIS[\{x_{k-1}^{(i)},w_{k-1}^{(i)}\}_{i=1}^N,y_{1:k}] [{xk(i),wk(i)}i=1N]=SIS[{xk1(i),wk1(i)}i=1N,y1:k]
F o r i = 1 : N For\quad i=1:N Fori=1:N
(1)采样: x k ( i ) ∼ q ( x k ( i ) ∣ x k − 1 ( i ) , y k ) x_k^{(i)}\sim q(x_k^{(i)}|x_{k-1}^{(i)},y_k) xk(i)q(xk(i)xk1(i),yk)
(2)根据 w k ( i ) ∝ w k − 1 ( i ) p ( y k ∣ x k ( i ) ) p ( x k ( i ) ∣ x k − 1 ( i ) ) q ( x k ( i ) ∣ x k − 1 ( i ) , y k ) w_{k}^{(i)}\propto w_{k-1}^{(i)}\frac{p(y_k|x_k^{(i)})p(x_k^{(i)}|x_{k-1}^{(i)})}{q(x_k^{(i)}|x_{k-1}^{(i)},y_k)} wk(i)wk1(i)q(xk(i)xk1(i),yk)p(ykxk(i))p(xk(i)xk1(i))递推计算各个粒子的权重;

重采样

重采样解决粒子权重退化的问题

Logo

有“AI”的1024 = 2048,欢迎大家加入2048 AI社区

更多推荐