obsinit.m

VARHMMBOX, version 1.1, Iead Rezek, Oxford University, MAR 2002 Matlab toolbox for Hidden Markov Mo

function [hmm] = obsinit (X,hmm,options)

% function [hmm] = obsinit (X,hmm,options)
%
% Initialise observation model in HMM
% 
% X         N x p data matrix
% hmm       hmm data structure
% options. 
%         covtype	'full' or 'diag' covariance matrices
%         gamma		weighting of each of N data points 
%         Bins          Number of bins for Dirichlet observation model
%         prrasc        prior range scale (scales prior variances);

obsmodelist={'Gauss','GaussCom','Poisson'};


K=hmm.K;
[T,ndim]=size(X);
if length(X)~=T,
  X=X';
  [T,ndim]=size(X);
end;

if ~isfield(hmm,'obsmodel')
  hmm.obsmodel='Gauss';
else
  if ~ismember(hmm.obsmodel,obsmodelist);
    error('Error: Undefined observation model');
  end;
end

defaultoptions=struct('covtype',[],'gamma',[]);
defaultoptions.covtype='full';
defaultoptions.gamma=ones(T,1);

if nargin<3
  options=defaultoptions;
else
  if ~isfield(options,'covtype'),
    options.covtype=defaultoptions.covtype;
  elseif ~ismember(options.covtype,{'full','diag'}),
    options.covtype=defaultoptions.covtype;
  end;
  if ~isfield(options,'gamma') 
    options.gamma=ones(T,1); 
  end
end;

hmm=initpriors(X,hmm,options);
hmm=initparams(X,hmm,options);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hmm] = initpriors(X,hmm,options)


[T,ndim]=size(X);

if ~isfield(options,'prrasc'),
  rangescale=1;
else
  rangescale=options.prrasc;
end;

% define priors
switch hmm.obsmodel
 case {'GaussCom','Gauss'},
  % priors
  midpoint=mean(X)';			% educated guess with scaling
  midscale=median(X);
  drange=range(X)'.*rangescale;

  defstatepriors.Norm_Mu=midscale;
  %    defstatepriors.Norm_Mu=midpoint;
  defstatepriors.Norm_Cov=diag(drange.^2);
  defstatepriors.Norm_Prec=inv(defstatepriors.Norm_Cov);
  defstatepriors.Wish_B=diag(drange);
  defstatepriors.Wish_alpha=ndim+1;
  defstatepriors.Wish_k=ndim;
 
 case 'Poisson',
  defstatepriors=struct('Gamma_alpha',[],'Gamma_beta',[]);
  defstatepriors.Gamma_alpha=1;
  defstatepriors.Gamma_beta=1;	
 otherwise
  error('Unknown observation model');
end;

% assigning default priors for observation models
if ~isfield(hmm,'state') | ~isfield(hmm.state,'priors'),
  for j=1:hmm.K,
    hmm.state(j).priors=defstatepriors;
  end;
else
  % priors not specified are set to default
  statepriorlist=fieldnames(defstatepriors)
  fldname=fieldnames(hmm.priors.state);
  misfldname=find(~ismember(statepriorlist,fldname));
  for j=1:hmm.K,
    for i=1:length(misfldname),
      priorval=getfield(defstatepriors,statepriorlist{i});
      hmm.state.priors=setfield(hmm.state,j,'priors',statepriorlist{i},priorval);
    end;
  end;
end;      

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [hmm] = initparams(X,hmm,options)

[T,ndim]=size(X);

covtype=options.covtype;
gamma=options.gamma;

% Initialising the posteriors
switch hmm.obsmodel
 case 'Gauss',
  mix=gmm(ndim,hmm.K,covtype);

  netlaboptions=foptions;
  netlaboptions(14) = 5; % Just use 5 iterations of k-means initialisation
  mix = gmminit(mix, X, netlaboptions);
  
  netlaboptions = zeros(1, 18);
  netlaboptions(1)  = 0;                % Prints out error values.
  
  % Termination criteria
  netlaboptions(3) = 0.000001;          % tolerance in likelihood
  netlaboptions(14) = 100;              % Max. Number of iterations.
  
  % Reset cov matrix if singular values become too small
  netlaboptions(5)=1;              
  [mix, netlaboptions, errlog] = wgmmem(mix, X, gamma,netlaboptions);
  hmm.gmmll=netlaboptions(8);     % Log likelihood of gmm model
  
  for k=1:mix.ncentres;
    hmm.state(k).Mu=mix.centres(k,:);
    switch covtype
     case 'full',
      hmm.state(k).Cov=squeeze(mix.covars(:,:,k));
      hmm.init_val(k).Cov = hmm.state(k).Cov; % In case we need to re-init
     case 'diag',
      hmm.state(k).Cov=diag(mix.covars(k,:));
      hmm.init_val(k).Cov = hmm.state(k).Cov; % In case we need to re-init
     otherwise,      
      disp('Unknown type of covariance matrix');
    end
  end
  hmm.train.init='gmm';
  hmm.mix=mix;
 case 'Poisson'
  mix=gmm(1,hmm.K,'diag');

  netlaboptions=foptions;
  netlaboptions(14) = 5; % Just use 5 iterations of k-means initialisation
  mix = gmminit(mix, X(:,2), netlaboptions);
  
  netlaboptions = zeros(1, 18);
  netlaboptions(1)  = 0;                % Prints out error values.
  
  % Termination criteria
  netlaboptions(3) = 0.000001;          % tolerance in likelihood
  netlaboptions(14) = 100;              % Max. Number of iterations.
  centres=zeros(1,hmm.K);
  
  for k=1:hmm.K,			% sorting: lowest rate first
    centres(k)=mix.centres(k,:);
  end
%  centres=sort(centres);
  for k=1:hmm.K,
    hmm.state(k).lambda=centres(k);		% Poisson rate
  end
 otherwise
  disp('Unknown observation model');
end
hmm.options=options;