vittest.m

来自「HMMBOX, version 3.2, William Penny, Impe」· M 代码 · 共 48 行

clear all
Ntest=1000;
for t=1:Ntest,
  T=750;
  K=2;
  covscale=11;
  simhmm=struct('Pi',rand(1,K),'P',rand(K),'K',K);
  simhmm.Pi=simhmm.Pi./sum(simhmm.Pi);
  simhmm.P=rdiv(simhmm.P,sum(simhmm.P,2))';
  simhmm.obsmodel='Gauss';
  for i=1:K,
    simhmm.state(i).Mu=rand(1,1);
    C=rand(2);
    simhmm.state(i).Cov=covscale*C'*C;
  end
  [simdata] = hmmsim (simhmm,T);
  hmm=struct('K',K);
  X=simdata.Xseries;
  hmm=hmminit(X,hmm,'full');
  hmm.train.cyc=30;
  hmm.obsmodel='Gauss';
  hmm.train.obsupdate=ones(1,hmm.K);    % update observation models ?
  hmm.train.init=1;         % Yes, we've already done initialisation
  
  hmm=hmmtrain(X,T,hmm);

  [simdata] = hmmsim (simhmm,T);
  X=simdata.Xseries;
  hmm.tain.obsupdate=zeros(1,hmm.K);
  hmm.train.pupdate=0;
  hmm.train.cyc=1;
  
  hmm=hmmtrain(X,T,hmm);
  [block]=hmmdecode(X,T,hmm);

  vit=block.q_star;
  [tmp,gamma]=max(hmm.gamma,[],2);
  data=[vit', gamma, simdata.Xclass];
  if data(1,1)~=data(1,3),
    data(:,3)=K+1-data(:,3);
  end
  subplot(211),plot(1:T,data(:,1),1:T,data(:,3));
  axis([1 T .9 K+.1])
  subplot(212),plot(1:T,data(:,2),1:T,data(:,3));
  axis([1 T .9 K+.1])
  H(t,1)=sum(data(:,1)~=data(:,3));
  H(t,2)=sum(data(:,2)~=data(:,3));
end;