[Scikit-learn] Dynamic Bayesian Network - Kalman Filter

看上去不錯的網站：http://iacs-courses.seas.harvard.edu/courses/am207/blog/lecture-18.html

SciPy Cookbook：http://scipy-cookbook.readthedocs.io/items/KalmanFiltering.html

 

---

 

良心視頻：卡爾曼濾波器的原理以及在matlab中的實現

講解思路貌似是在已知迭代結果的基礎上作講解，不是很透徹。

1. 用矩陣表示

2. 本質就是：二維高斯的協方差與sampling效果

3. 不肯定性在狀態之間的傳遞

4. 矩陣表示觀察數據

5. Kalman係數

6. 噪聲協方差矩陣的更新

7. Matlab實現

 

思考： 

與數學領域 openBUGS 的估參的關係是什麼？[Bayes] openBUGS: this is not the annoying bugs in programming

一個是對逐漸增多數據的實時預測；一個是對整體數據的迴歸擬合。

 

代碼示例：純python代碼

# Kalman filter example demo in Python

# A Python implementation of the example given in pages 11-15 of "An # Introduction to the Kalman Filter" by Greg Welch and Gary Bishop, # University of North Carolina at Chapel Hill, Department of Computer # Science, TR 95-041, # http://www.cs.unc.edu/~welch/kalman/kalmanIntro.html

# by Andrew D. Straw

import numpy as np import matplotlib.pyplot as plt plt.rcParams['figure.figsize'] = (10, 8) # intial parameters
n_iter = 50 sz = (n_iter,) # size of array
x  = -0.37727  # truth value (typo in example at top of p. 13 calls this z)
z  = np.random.normal(x,0.1,size=sz) # observations (normal about x, sigma=0.1)
# 已得到一組隨機數

Q = 1e-5 # process variance

# allocate space for arrays
xhat =np.zeros(sz)         # a posteri estimate of x
P    =np.zeros(sz)         # a posteri error estimate
xhatminus =np.zeros(sz)    # a priori estimate of x
Pminus    =np.zeros(sz)    # a priori error estimate
K         =np.zeros(sz)    # gain or blending factor
 R = 0.1**2                 # estimate of measurement variance, change to see effect

# intial guesses
xhat[0] = 0.0 P[0] = 1.0

# 開始迭代
for k in range(1, n_iter): # time update
    xhatminus[k] = xhat[k-1] Pminus[k] = P[k-1]+Q # measurement update
    K[k]    = Pminus[k]/( Pminus[k]+R ) xhat[k] = xhatminus[k]+K[k]*(z[k]-xhatminus[k]) P[k] = (1-K[k])*Pminus[k] plt.figure() plt.plot(z,'k+',label='noisy measurements') plt.plot(xhat,'b-',label='a posteri estimate') plt.axhline(x,color='g',label='truth value') plt.legend() plt.title('Estimate vs. iteration step', fontweight='bold') plt.xlabel('Iteration') plt.ylabel('Voltage') plt.figure() valid_iter = range(1,n_iter) # Pminus not valid at step 0
plt.plot(valid_iter,Pminus[valid_iter],label='a priori error estimate') plt.title('Estimated $\it{\mathbf{a \ priori}}$ error vs. iteration step', fontweight='bold') plt.xlabel('Iteration') plt.ylabel('$(Voltage)^2$') plt.setp(plt.gca(),'ylim',[0,.01]) plt.show()

Result: 

 

Goto: [OpenCV] Samples 14: kalman filter

其實，真正的Kalman Filter用得是以下理論，上述例子只是教小學生的入門讀物。

 

Goto: https://www.youtube.com/watch?v=UVNeulkWWUM by XU Yida

關鍵須要理解： http://www.cnblogs.com/rubbninja/p/6220284.html

 

【重點】證實過程的理解關鍵是：

由於是線性濾波器，自己又具有一個alpha迭代的過程，那麼先找出joint distribution，

而後，根據高斯的性質直接得出條件機率，便是Update Rule，這樣正好對應於濾波器的alpha迭代過程的形式。

這個條件機率就是關於x_t的，也就是最新的狀態的機率分佈，那麼指望也就是miu，就是最新的x_t。

大概就是這麼個思路，筆記在本本上，具體請看視頻。符號比較多，但大致就是如上脈絡。