从零开始的机器学习--概率图模型--隐马尔科夫模型

本文参考机器学习周志华



main.m

clear;clc;
P = [0.2,0.4,0.4];
n = 3;%隐状态个数
A = [0.5,0.2,0.3;0.3,0.5,0.2;0.2,0.3,0.5];
B = [0.5,0.5;0.4,0.6;0.7,0.3];
k = 2;%观测种类数
% B = repmat(ones(1,k)./k,n,1);
% A = repmat(ones(1,n)./n,n,1);
% O = [1,2,1,2,1,1,1];
% Label = [1,2,3,1,2,3,1];
O = [1,2,1];
% O = [1,2,1,2,1,2,1,1,1,2,2,2];
% O = [2,1,2,1,2,1,2,2,2,1,1,1];
% O = [2,1,2];
t = length(O);
Label = [1,2,3];
%% 计算乾后项概率
% 
% alpha = zeros(t,n);
%  for i = 1:t
%     if i == 1
%         alpha(i,:) = P.*B(:,O(i))';
%     else
%         alpha(i,:) = alpha(i-1,:)*A.*B(:,O(i))';
%     end
%  end
% % p1 = sum(alpha,2);
% %计算后向概率
% beta = zeros(t,n);
% for i = t:-1:1
%     if i == t
%         beta(i,:) = 1;
%     else
%         beta(i,:) =  beta(i+1,:).*B(:,O(i+1))'*A';
%     end
% end
% % p = sum(alpha.*beta,2);

%% 迭代EM算法

% ITER = 10;
% t = length(O);
% for Times = 1:ITER
%     alpha = zeros(t,n);
%     for i = 1:t
%         if i == 1
%             alpha(i,:) = P.*B(:,O(i))';
%         else
%             alpha(i,:) = alpha(i-1,:)*A.*B(:,O(i))';
%         end
%     end
%     beta = zeros(t,n);
%     for i = t:-1:1
%         if i == t
%             beta(i,:) = 1;
%         else
%             beta(i,:) =  beta(i+1,:).*B(:,O(i+1))'*A';
%         end
%     end
%     PO = sum(alpha(1,:).*beta(1,:));
%     ksi = zeros(n);
%     gamma = zeros(n,1);
%     a_num = zeros(n);
%     a_den = zeros(n,1);
%     b_num = zeros(n,k);
%     for se = 1:t-1
%         temp1 = beta(O(se+1),:).*B(:,O(se+1))';
%         temp2 = alpha(O(se),:).*A';
%         ksi = (temp1.*temp2')./PO;
%         clear temp1 temp2
%         gamma = sum(ksi,2);
%         if se == 1
%             P = gamma';
%         end
%         a_num = a_num+ksi;
%         a_den = a_den+gamma;
%         b_num(:,O(se)) = b_num(:,O(se))+gamma;
%     end
%     A = a_num./a_den;
%     B = b_num./a_den;
% end

%% 维比特解码
viterbipath = zeros(n,t);
viterbiprob = zeros(n,t);
% viterbipath(:,1) = 1:n;
viterbiprob(:,1) = P'.*B(:,O(1));
% for i = 2:t
%     for state = 1:n
% %         tempp = viterbiprob(:,i-1)'.*A(state,:)*B(state,O(i));
%         tempp = viterbiprob(:,i-1)*A;
%         [val,row] = max(tempp);
%         delta = B(:,O(i+1)).*val
%         [~,col] = find(tempp==max(tempp));
%         viterbipath(state,i) = col;
%         viterbiprob(state,i) = max(tempp);
%     end
% end
for i = 2:t
    tempp = (viterbiprob(:,i-1)'.*A')';
    [val,row] = max(tempp);
    viterbiprob(:,i) = B(:,O(i)).*val';
    viterbipath(:,i) = row';
end
bestpathindex = find(viterbiprob(:,end)==max(viterbiprob(:,end)));
bestpath = [viterbipath(bestpathindex,:),bestpathindex];
bestpath = bestpath(:,2:end);