详解HMM模型原理 及 实现（之四：matlab实现曲线分类）

本文详解paper “A Tutorial on HMMs and Selected Applications in Speech Recognition"并进行matlab实现（尽量用其他编程语言通用的实现）
**为方便阅读代码，注释部分用python的”#"

实现HMM模型用于简单的曲线分类
先定义HMM变量
曲线在T个时刻取值
observation定义为曲线可以取的值，这里先把曲线做归一化，然后从0到1之间以0.1为间隔切为10个区间，依次标号1到10，那么observation就有10个，为1到10的值
train_data定义为训练的曲线集，集中到一个矩阵，每行为一个observation序列，行数为train_data的数量
N为state数，这里设为5

   maxIter = 100;  #迭代次数
   N = 5; #number of states
   M = 10; #number of observations
   T = size(train_data, 2);  #training data列数应一致

初始化HMM参数，可以用均匀分布，也可以用随机数



这里pai和A用均匀分布，B用随机分布

   pai = 1/N * ones(N,1);
   A = 1/N * ones(N,N);  #state transition matrix
   B = rand(N,M); #row: states, col: observations

初始化计算过程中的变量alpha, beta, gamma, xi
xi的size本来应该是state * state * T
但注意在计算aij时只用到了xi的1到T-1时刻求和，因此减掉一个维度T，直接定义为1到T-1时刻求和的xi, 变量定义为xi_sumT_1

   alpha = zeros([N, T]); #row: states, col: time, for calculating problem 1
   beta = zeros([N, T]); #row: states, col: time, for calculating problem 1
   gamma = zeros([N, T]); #row: states, col: time, for calculating problem 2
   #xi = zeros([N*N, T-1]); #row: state*state, col:time-1, for calculating problem 3 
   xi_sumT_1 = zeros(N, N); #直接对t时刻求和 row:states, col:states, for calculating problem 3
   trainNum = size(train_data, 1);

另外一些辅助变量

   iter = 1;  #迭代次数计数器
   converged = 0; #0:未收敛 1:收敛
   previous_loglik = -inf; #记录前一次likelihood
   converge_thres = 1e-4; #判断收敛的阈值

HMM设计的3个基本问题：

1. 给定一个HMM模型，一个观测信号序列，估计出序列的概率
   给出observation序列, 模型lambda=(A,B,pi), 计算P(O | lambda)
   即模型与Observation有多匹配
2. 决定最好的模型states
   给出observation序列，模型lambda=(A,B,pi), 选择state序列Q=q1q2…qT
   即探索隐性状态
3. 调整模型参数来最好地解释观测到的信号
   调整lamba=(A,B,pi)使P(O | lambda)最大
   即用observation进行训练

这里要训练的HMM模型是用来识别曲线的，每个曲线类型都要训练一个模型(problem3)，而后输入测试曲线，计算在每个模型下的概率(problem1)，概率最大的为匹配类型，达到分类效果
每个类型的模型训练都是相同的方法，因此具体写一个类型的train_data

problem3中需要用到的变量为xi

而计算xi需要用到状态转移矩阵A，state下observation概率矩阵b, problem1中的alpha和beta
state下observation概率矩阵b的计算：
根据state下observation概率矩阵B
给定一个训练曲线，这个曲线在1到T时刻有值，取值范围在observation内
那么根据每个时刻t处的observation, 和指定的state，可以在B矩阵中查找到概率，可以知道b的size为state * T, 例如state为j, 时刻t时的值为bj(Ot)

for t=1:T
  b(:,t) = B(:, train_data