隐马尔可夫模型（HMM）的MATLAB实现——Baum-Welch算法

Baum-Welch算法用来对隐马尔可夫模型的参数进行学习，Baum-Welch算法是EM算法的一种特例，属于非监督学习算法，下边第一部分程序的迭代次数为1的算法代码，其中调用了计算中所需要的两个变量Gamma和Xi的计算代码，一并在此给出，第二部分为迭代n次的代码，代码中已经给出了示例。

1，迭代次数为1的Baum-Welch算法代码

function [A_1,B_1,Pi_1] = BaumWelchAlgo(A,B,Pi,O)
% 函数功能：实现迭代次数为1的Baum-Welch算法
%
% 参考文献：1，李航《统计学习方法》
%          2，Bilmes L,et al. A Gentle Tutorial of EM Algorithm and its Application to
%               Parameter Estimation for Gaussian Mixture and Hidden Markov Models
%          3，Rabiner L. A Tutorial on Hidden Markov Models and Selected
%               Applications in Speech Recognition
%
% 思路：
% pi_i = Gamma_1(i);
% a_ij = sum_t-from-1-to-(T-1)(Xi_t(i,j)) / sum_t-from-1-to-T(Gamma_t(i));
% b_j(k) = sum_t-from-1-to-T_ot-=-vk(Gamma_t(i)) / sum_t-from-1-to-T(Gamma_t(j));
% 其中：
% sum_i-from-1-to-N(pi_i) = 1
% sum_j-from-1-to-N(a_ij) = 1
% sum_k-from-1-to-M(b_j(k)) = 1
% 
% 输入：
%     第一次选取的模型参数A,B,Pi和观测序列O。
% 输出：
%     经过一次迭代Baum-Welch算法之后得到的模型参数A_1,B_1,Pi_1
%     
% 示例：
% 输入
% A =
%     0.5000    0.2000    0.3000
%     0.3000    0.5000    0.2000
%     0.2000    0.3000    0.5000
% B =
%     0.5000    0.5000
%     0.4000    0.6000
%     0.7000    0.3000
% Pi =
%     0.2000
%     0.4000
%     0.4000
% O =
%      1
%      2
%      1
%      2
% 输出：
% >>[a,b,c] = BaumWelchAlgo(A,B,Pi,O)
% 
% a =
%     0.5151    0.2212    0.2637
%     0.2987    0.5231    0.1783
%     0.2210    0.3674    0.4116
% b =
%     0.4325    0.5675
%     0.4257    0.5743
%     0.6389    0.3611
% c =
%     0.1872
%     0.3241
%     0.4887
%
A_size = size(A);
B_size = size(B);
O_size = size(O);
N = A_size(1,1);%状态集个数
M = A_size(1,2);
Y = B_size(1,2);%B的列数，即观测集个数
K = O_size(1,1);
%
% ----调用GammaXiAlgo函数，计算Gamma矩阵和Xi矩阵----
[Gamma,Xi] = GammaXiAlgo(A,B,Pi,O);
%
% ---------------计算转移概率矩阵A_1---------------
% 转移概率a_ij = 在观测O下由状态i转移到状态j的期望值 / 在观测O下由状态i转移的期望值
%             = expected number of transitions from state Si to state Sj /
%             excepted number of transitions from state Si
A_1 = zeros();
f = sum(Gamma,2);%计算Gamma的每一行之和，表示在观测时间序列的所有时刻都处于某一个状态的概率和
for i = 1:M
    for j = 1:N
        A_1(i,j) = sum(sum(Xi(i,j,:))) / (f(i,1) - Gamma(i,K));
    end
end
%
% ---------------计算观测概率矩阵B_1---------------
% 观测概率b_j(k) = 满足条件(ot = vk)的观测O下由状态i转移的期望值 / 观测O下由状态i转移的期望值
%               = excepted number of times in state j and observing symbol
%               vk / excepted number of times in state j
B_1 = zeros();
G = zeros();
c = 0;
for y = 1:Y
    for j = 1:N
        for t = 1:K
            if O(t,1) == y % 用于比较输出的观测值是否由对应的状态值输出，即条件ot = vk;
                c = c + 1;
                G(c,1) = Gamma(j,t);
            end
        end
         B_1(j,y) = sum(G) / f(j,1);
         c = 0;
         G = zeros();
    end
end
%
% ---------------计算初始概率pi_1---------------
Pi_1 = zeros();
for i = 1:M
    Pi_1(i,1) = Gamma(i,1);
end

2，变量Gamma和Xi的计算代码

function [Gamma,Xi] = GammaXiAlgo(A,B,Pi,O)
% 函数功能：计算Baum-Welch算法公式中的gamma概率和xi概率
%
% 参考文献：李航《统计学习方法》
%
% 思路：
% gamma_t(i) = P(it = qi|O,lambda) = P(it = qi,O|lambda) / P(O|lambda);
% xi_t(i) = P(i_t = qi,i_t+1 = qj|O,lambda) = P(i_t = qi,i_t+1 = qj,O|lambda) / P(O|lambda)
% gamma_t(i)表示在时刻t处于状态qi的概率
% xi_t(i,j)表示在时刻t处于状态qi，且在时刻t+1处于状态qj的概率
%     其中t的范围是观测时间段长度，即序列O的长度K
%     qi的范围是状态集合中状态的个数，即转移矩阵A的维数N
% 
% 输入：隐马尔可夫模型lambda(A,B,Pi)，观测序列O；
%      其中A为转移概率，B为观测概率（行数表示状态集合中的元的个数，列数表示观测集中的元的个数）；
%      Pi为初始概率，初始概率Pi和观测序列O须为列矩阵。
% 输出：矩阵Gamma：存放不同时刻t处于不同状态qi的gamma概率值，为N*K阶
%       矩阵Xi：存放不同时刻t处于不同状态qi且t+1处于不同状态qj的概率值，为N*N*K阶
%             其中N*N维矩阵中，列对应t时刻，行对应t+1时刻，第K维表示的是时间序列
%             示例：
%             y(:,:,2) =
%             0.1474    0.0523    0.1157
%             0.1268    0.1874    0.1106
%             0.0464    0.0617    0.1518
%             那么，y(:,:,2)中的2表示此时处在时刻t=2；
%             0.1474表示在时刻2处于状态1且在时刻3处于状态1的概率；
%             0.1268表示在时刻2处于状态2且在时刻3处于状态1的概率；
%             0.0523表示在时刻2处于状态1且在时刻3处于状态2的概率；
%             以此类推。
%
% 示例：
% 输入：
% A =
%     0.5000    0.2000    0.3000
%     0.3000    0.5000    0.2000
%     0.2000    0.3000    0.5000
% B =
%     0.5000    0.5000
%     0.4000    0.6000
%     0.7000    0.3000
% Pi =
%     0.2000
%     0.4000
%     0.4000
% O2 =
%      1
%      2
%      1
%      2
% 输出：
% >>[x,y] = GammaXiAlgo(A,B,Pi,O2)
% x =
%     0.1872    0.3154    0.3205    0.3508
%     0.3241    0.4248    0.3014    0.4192
%     0.4887    0.2599    0.3781    0.2301
% y(:,:,1) =
%     0.1024    0.0462    0.0386
%     0.0983    0.1847    0.0412
%     0.1147    0.1939    0.1801
% y(:,:,2) =
%     0.1474    0.0523    0.1157
%     0.1268    0.1874    0.1106
%     0.0464    0.0617    0.1518
% y(:,:,3) =
%     0.1742    0.0836    0.0627
%     0.0886    0.1773    0.0355
%     0.0879    0.1583    0.1319
%
P1 = zeros();
A_size = size(A);
O_size = size(O);
N = A_size(1,1);
M = A_size(1,2);
K = O_size(1,1);
%
% 调用前向后向算法函数ForwardBackwardAlgo_1，计算得到alpha矩阵P1，beta矩阵P2和观测序列O的概率p1
[P1,P2,p1] = ForwardBackwardAlgo(A,B,Pi,O);
%
% t表示观测序列时间长度数，i表示状态数；
% -------------计算单个状态概率gamma = Gamma--------------
Gamma = zeros();
for t = 1:K
    for i = 1:M
        Gamma(i,t) = P1(i,t) * P2(i,t) / p1;
    end
end
%
% --------------计算两个状态概率xi = Xi-----------------
Xi = zeros();
for t = 1:K-1
    for i = 1:M
        for j = 1:N
           Xi(i,j,t) = P1(i,t) * A(i,j) * B(j,O(t+1,1)) * P2(j,t+1) / p1;
        end
    end
end

3，迭代次数为n的Baum-Welch算法

function [A_n,B_n,Pi_n] = BaumWelchAlgo_n(A,B,Pi,O,n)
% 函数功能：n次迭代的BaumWelchAlgo_n算法
for r = 1:n
     [A,B,Pi] = BaumWelchAlgo(A,B,Pi,O);
end
A_n = A;
B_n = B;
Pi_n = Pi;