mk_arhmm_obs_lik.m

上载文件为Matlab环境下的高斯以马尔科夫模型例程

function B = mk_arhmm_obs_lik(data, mu, Sigma, W, mu0, Sigma0)
% MK_ARHMM_OBS_LIK Make the observation likelihood vector for an autoregressive(1) HMM.
% B = mk_arhmm_obs_lik(data, mu, Sigma, W, mu0, Sigma0)
%
% P(Y(t)=y | Y(t-1)=u, Q(t)=j) = N(y; mu(:,j) + W(:,:,j) * u, Sigma(:,:,j))
% P(Y(1)=y | Q(1)=j) = N(y; mu0(:,j), Sigma0(:,:,j))
% y(t) = data(:,t) = observation vector at time t
%
% Output:
% B(i,t) = Pr(y(t) | Q(t)=i, y(t-1))

[d Q] = size(mu);
[d T] = size(data);
B = zeros(Q,T);

for j=1:Q
  X = data(:, 2:T);
  M = repmat(mu(:,j), 1, T-1) + W(:,:,j) * data(:, 1:T-1);
  C = Sigma(:,:,j);
  Cinv = inv(C);
  denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
  % Given X(i,t) and ignoring M, [X'*C](t,j) = sum_i X(t,i) C(i,j)
  % [X'*C*X'](t,j) = (sum_i X(t,i) C(i,j)) .* X(t,j)
  % mahal(t,1) = sum_i sum_j X(t,i) C(i,j) X(t,j)
  mahal = sum(((X-M)'*Cinv).*(X-M)',2); 
  B(j,2:T) = exp(-0.5*mahal') / denom;
end

[d Q0] = size(mu0);
for j=1:Q0
  x = data(:, 1);
  m = mu0(:,j);
  C = Sigma0(:,:,j);
  Cinv = inv(C);
  denom = (2*pi)^(d/2)*sqrt(abs(det(C)));
  mahal = sum(((x-m)'*Cinv).*(x-m)',2); 
  B(j,1) = exp(-0.5*mahal) / denom;
end