forward_backward.c

EM算法的MATLAB实现,EM算法是HMM中用于随机过程模型参数训练的经典算法

/*

 Forward Backward HMM algorithm


 Usage  
 -------


 [X , gamma , logll , [sum_alpha] , [sum_beta]]  = forward_backward(PI , A , L , [filter]);


 Inputs
 -------

 PI            Initial proabilities (N x 1) : Pr(x_1 = i) , i=1,...,N

 A             State transition probabilities matrix Pr(x_{k} = i| x_{k - 1} = j) such
               sum_{x_k}(A) = 1 => sum(A , 1) = 1

 L             Time indexed Likelihood matrix Pr(z_k | x_k = i) (N x K), i=1,...,N, k=1,...,K. 
               Extracted from B matrix such that B = Pr(z | x) (M x N), sum(B , 1) = 1 and B(z_k , :)' = L(: , k).

 filter        Optional flag. If filter = 0 => Just the alpha probabilities are computed, else (default filter = 1)
               the two passes alpha & beta.


 Ouputs
 -------

 X             MAP sequence estimated from gamma probabilies Pr(x_k|Z_1:K)

 gamma         gamma probabilities (N x K) or alpha probabilities if filter = 0

 logll         Loglikelihood of (PI,A,L)

 If the -Dbetanormalize flag is used to compile, two more outputs are given :

 sum_alpha     Pr(z_k|Z_{1:k-1})  (1 x K)

 sum_beta      Pr(z_k|Z_{k+1:K})  (1 x K)


 To compile
 ----------


mex (-Dbetanormalize) -output forward_backward.dll forward_backward.c

The betanormalize flag forces to scale gamma's probabilities with the sum of beta instead of the alpha's.


mex -f mexopts_intel10amd.bat -output forward_backward.dll forward_backward.c

mex -Dbetanormalize -f mexopts_intel10amd.bat -output forward_backward.dll forward_backward.c


Example 1
----------


N                        = 256;
M                        = 256;
K                        = 300;

PI                       = rand(N , 1);
sumPI                    = sum(PI);
PI                       = PI./sumPI(ones(N , 1) , :);

A                        = rand(N , N);
sA                       = sum(A);
A                        = A./sA(ones(N , 1) , :);

L                        = rand(N , K);

[X , gamma , logll]      = forward_backward(PI , A , L);


Example 2
----------


N                        = 2;
M                        = 4;
K                        = 30;
PI                       = rand(N , 1);
sumPI                    = sum(PI);
PI                       = PI./sumPI(ones(N , 1) , :);

A                        = rand(N , N);
sumA                     = sum(A);
A                        = A./sumA(ones(N , 1) , :);

B                        = rand(M , N);
sumB                     = sum(B);
B                        = B./sumB(ones(M , 1) , :);

[Z , X]                  = sample_dhmm(K , PI , A , B);

L                        = densdis(Z , B); %Likelihood

[X_MAP , gamma , logll]  = forward_backward(PI , A , L);


Example3
----------


d                                      = 2;
M                                      = 256;
K                                      = 300;
W                                      = [1 ; M];
H                                      = 1;
R1                                     = 2;
R2                                     = 3;


PI                                     = [0.5 ; 0.5];
A                                      = [0.9 0.1 ; 0.1 0.9];
X                                      = sample_dmc(K , PI , A);


PI1                                    = iniPI(W);
A1                                     = matA_gauss(W , H , R1);
B1                                     = A1;


PI2                                    = iniPI(W);
A2                                     = matA_gauss(W , H , R2);
B2                                     = A2;


[Z1 , X1]                              = sample_dhmm(K , PI1 , A1 , B1);
[Z2 , X2]                              = sample_dhmm(K , PI2 , A2 , B2);

temp                                   = [Z1 ; Z2];
Z                                      = temp((0:K-1)*d + X);

L1                                     = densdis(Z , B1);
L2                                     = densdis(Z , B2);

[X1 , gamma1 , ll1 , calpha1 , cgamma1] = forward_backward(PI1 , A1 , L1);
[X2 , gamma2 , ll2 , calpha2 , cgamma2] = forward_backward(PI2 , A2 , L2);

L                                      = [calpha1.*cgamma1 ; calpha2.*cgamma2];
X_FB                                   = forward_backward(PI , A , L);


 Author : S閎astien PARIS