[Math] Hidden Markov Model

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

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

良心视频：卡尔曼滤波器的原理以及在matlab中的实现

讲解思路貌似是在已知迭代结果的基础上做讲解，不是很透彻。

1. 用矩阵表示

[Math] Hidden Markov Model

2. 本质就是：二维高斯的协方差与sampling效果

[Math] Hidden Markov Model

3. 不确定性在状态之间的传递

[Math] Hidden Markov Model

4. 矩阵表示观察数据

[Math] Hidden Markov Model

5. Kalman系数

[Math] Hidden Markov Model

6. 噪声协方差矩阵的更新

[Math] Hidden Markov Model

7. Matlab实现

[Math] Hidden Markov Model

思考：

与数学领域 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()