dempoissonhmm.m

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

% A demonstration of the HMM software using a Poisson observation
% model on count data
clear all

load dempoisson
refval=50;
scaledcountdat=normalise(countdat,'a',refval);
scaledcountdat=round(scaledcountdat+.5);

Xtrain=[scaledcountdat];
Xtrain=[ones(length(Xtrain),1) Xtrain];
T=length(Xtrain);

%plot(cumsum(countdat),countdat);
%title('Original data'),xlabel('Time (msec)'),ylabel('RR Interval (msec)');
disp(' The data consist of a sequence of R-R intervals obtained from an');
disp(' ECG recording. It can be considered a sequence of counts, the');
disp(' counts being the number of sample intervals between each');
disp(' heart-beat. In this particular case, the subject is asked to');
disp(' take 6 deep breaths at 30secs into the recoding (lasting until');
disp(' 90 seconds into the recording');
disp(' ');
disp('Press a key to continue');
pause


disp(' ');
disp('We will take random samples to initialise the HMM.');
disp(' ');
disp('Press a key to continue');
pause

% Train up HMM on this data
hmm=struct('K',6);

disp(' ');
hmm=hmmhsinit(Xtrain,hmm);
hmm.obsmodel='Poisson';
obsoptions=struct('prrasc',1);
hmm=obsinit(Xtrain,hmm,obsoptions);

disp('Mean rates of HMM initialisation');
hmm.state(1).lambda/refval*60
hmm.state(2).lambda/refval*60

% Train up HMM on observation sequence data using Baum-Welch
% This uses the forward-backward method as a sub-routine
disp('We will now train the HMM using Baum/Welch');
disp(' ');
disp('Press a key to continue');
pause
disp('Estimated HMM');

hmm.train.cyc=50;
hmm.train.tol=-inf;
hmm.train.obsupdate=ones(1,hmm.K);    % update observation models ?
hmm.train.init=1;         % Yes, we've already done initialisation

hmm=hmmtrain(Xtrain,T,hmm);
disp('Rates');
hmm.state(1).lambda/refval*60
hmm.state(2).lambda/refval*60
disp('Initial State Probabilities, Pi');
hmm.Pi
disp('State Transition Matrix, P');
hmm.P

[block]=hmmdecode(Xtrain,T,hmm);


% Find most likely hidden state sequence using Viterbi method
[block]=hmmdecode(Xtrain,T,hmm);


% sorting the labels s.t. lowest starts first
lambda=zeros(1,hmm.K);
for i=1:hmm.K lambda(i)=hmm.state(i).lambda; end;
[tmp,ndx]=sort(lambda);

block.sq_star=zeros(size(block.q_star));
l=block.q_star(1);
for i=1:hmm.K
  ndx2=find(block.q_star==ndx(i));
  block.sq_star(ndx2(:))=i*ones(1,length(ndx2));
end;

figure
subplot(211),plot(cumsum(countdat),block.sq_star);
title('Viterbi decoding');
av=axis;av(3:4)=[0 hmm.K]+0.5;axis(av);
subplot(212),plot(1:length(hrv),hrv*60);
title('Interpolated HR signal');