隐马尔可夫模型（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.4887A_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 SiA_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));    endend% ---------------计算观测概率矩阵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 jB_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();    endend% ---------------计算初始概率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.1319P1 = 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;    endend% --------------计算两个状态概率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    endend

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);endA_n = A;B_n = B;Pi_n = Pi;

下一篇介绍HMM解码程序-维特比算法。