hmmposterior.m

Matlab 马尔科夫计算工具包

function [posterior_prob]=hmmposterior(data,hmm);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [posterior_prob]=hmmposterior(data,hmm)
%
% computes the posterior probablity of a data set, given
% a HMM prob. density:
% hmm.state.Mu:              state centres
% hmm.state.Cov:           state covariances
% hmm.prior:                  state prior probabilities; (dafault: flat)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
if nargin<2,
  error('Insufficient number of arguments!');
end;

[N,dim]=size(data);
if N<dim,
  data=data';
  [N,dim]=size(data);
end;

if ~isfield(hmm,'state'),
  error('No states specified!');
end;

if ~isfield(hmm.state,'Mu'),
  error('No states centers specified');
elseif length(hmm.state(1).Mu)~=dim,
    error('Dimensions of centres do not match data!');
end;

if ~isfield(hmm.state,'Cov'),
  error('No states centers specified');
elseif size(hmm.state(1).Cov,1)~=size(hmm.state(1).Cov,2),
  error('Covariance matrix must be square!');
elseif  size(hmm.state(1).Cov,1)~=dim,
  error('Dimensions of covariance matrix do not match data!');
end;

if ~isfield(hmm,'prior'),
  disp('Using flat prior');
  hmm.prior=ones(1,length(hmm.state));
  hmm.prior=hmm.prior./length(hmm.prior);
else
  if length(hmm.prior)~=length(hmm.state),
    error('Number of priors must be identical to number of states');
  end;
end;

K=length(hmm.state);
posterior_prob=zeros(N,K);

for i=1:K,
  posterior_prob(:,i)=gausspdf(data,hmm.state(i).Mu,hmm.state(i).Cov);
  posterior_prob(:,i)=posterior_prob(:,i).*hmm.prior(i);
end;
S=sum(sort(posterior_prob,2),2);
ndx=find(S==0);
posterior_prob(ndx,:)=K.*ones(length(ndx),K);
S=sum(sort(posterior_prob,2),2);
posterior_prob=posterior_prob./(S*ones(1,K));

return;
  


function [px] = gausspdf (x,mu,sigma)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   [px] = gausspdf (x,mu,sigma)
%
%   computes m-D probability density from x values given
%   mean mu and standard deviation sigma
%              1                                     T       -1
%   p(x)= ------------------------- exp (-0.5  (x-mu)   Cov    (x-mu)  )
%         (2*pi)^(d/2) |Cov|^0.5
%
%  e.g: [X,Y] = meshdom(-2:.2:2, -2:.2:2);
%       mu=[0;0];sigma=[1 0;0 1];
%       for i=1:21, for j=1:21, 
%             p(i,j)=gaussmd([ X(i,j) Y(i,j)],mu,sigma); 
%       end; end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

d=size(sigma,1);
if (size(x,2)~=d)  x=x'; end;
if (size(mu,2)~=d)  mu=mu'; end;

N=size(x,1);
ndim=size(x,2);


%px=zeros(N,1);

z=(x-ones(N,1)*mu);
IS=inv(sigma);
DS=det(sigma);

px=exp(-0.5.*sum((z*IS).*z,2))./sqrt((2*pi)^ndim*DS);
return;