learn_dhmm_entropic1.m

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

function [hmm, LL] = learn_dhmm_entropic(data, hmm, varargin)
% LEARN_DHMM_ENTROPIC Find the MAP params of an HMM with discrete outputs with an entropic prior using EM
%
% [hmm, LL] = learn_dhmm_entropic(data, hmm, ...)
%
% This has the same interface as learn_dhmm_simple.
%
% Extra optional params
% 'trimtrans' - trim uninformative outgoing transitions? [0]
% 'trimobs' - trim uninformative observations? [0]
% 'trimstates' - trim low occupancy states? [0]
% 'anneal' - do deterministic annealing? [0]

% Based on "Structure learning in conditional probability models via an entropic  prior
% and parameter extinction", M. Brand, Neural Computation 11 (1999): 1155--1182
% For the annealed case, see "Pattern discovery via entropy minimization",
% M. Brand, AI & Statistics 1999. Equation numbers refer to this paper.
        
max_iter = 30;
thresh = 1e-2;
verbose = 1;
dirichlet = 0;
trimtrans = 0;
trimobs = 0;
trimstates = 0;
anneal = 0;

if nargin >= 3
  args = varargin;
  for i=1:2:length(args)
    switch args{i},
     case 'max_iter', max_iter = args{i+1};
     case 'thresh', thresh = args{i+1};
     case 'verbose', verbose = args{i+1};
     case 'dirichlet', dirichlet = args{i+1};
     case 'trimtrans', trimtrans = args{i+1};
     case 'trimobs', trimobs = args{i+1};
     case 'trimstates', trimstates = args{i+1};
     case 'anneal', anneal = args{i+1};
    end
  end
end      

previous_loglik = -inf;
loglik = 0;
converged = 0;
num_iter = 1;
LL = [];

if ~iscell(data)
  data = num2cell(data, 2); % each row gets its own cell
end
numex = length(data);

startprob = hmm.startprob;
endprob = hmm.endprob;
transmat = hmm.transmat;
obsmat = hmm.obsmat;


if anneal
  % schedule taken from Ueda and Nakano, "Determinsitic Annealing EM algorithm",
  % Neural Networks 11 (1998): 271-282, p276
  b = [];
  temp = [];
  i = 1;
  b(i)=0.1;
  temp(i)=1/b(i);
  while b(i) < 1
    i = i + 1;
    b(i)=b(i-1)*1.2;
    temp(i)=1/b(i);
  end
  temp_schedule = temp;
end

Q = hmm.nstates;
O = hmm.nobs;

% record what has already been trimmed
trimmed_trans = zeros(1,Q);
trimmed_obs = zeros(1,Q);
trimmed_states = zeros(1,Q);

while (num_iter <= max_iter) & ~converged
  % Z = 1 is the min entropy case, Z = 0 is ML, Z = -1 is max ent
  % Z << 0 is the high temperature case
  if anneal
    if num_iter <= length(temp_schedule)
      temp = temp_schedule(num_iter);
    else
      temp = temp_schedule(end);
    end
    T0 = 0;
    if temp <= 1.0
      Z = 1;
    else
      Z = T0 - temp;
    end
  else
    Z = 1;
  end

  % E step
  [loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ...
      compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet);

  converged = em_converged(loglik, previous_loglik, thresh);
  if converged
    Z = 1; % do the last step with min entropy
  end

  if verbose, fprintf(1, 'iteration %d, loglik = %7.4f, Z=%5.3f\n', num_iter, loglik, Z); end
  num_iter =  num_iter + 1;

  previous_loglik = loglik;
  LL = [LL loglik];

  % M step
  startprob = normalise(exp_num_visits1);
  endprob = normalise(exp_num_visitsT);
  
  %transmat = mk_stochastic(exp_num_trans);
  for i=1:Q
    ndx = find(transmat(i,:)==0);
    assert(all(exp_num_trans(i,ndx)==0))
    transmat(i,:) = entropic_map(exp_num_trans(i,:), Z);
    assert(all(transmat(i,ndx)==0))
    
    % only trim if we are in the min entropy setting
    % If Z << 0, we would trim everything!
    if trimtrans & ~trimmed_trans(i) & (Z==1) 
      % grad(j) = d log lik / d theta(i ->j)
      % transmat(i,j) = 0 => exp_num_trans(i,j) = 0
      % so we can safely replace 0s by 1s in the denominator
      denom = transmat(i,:) + (transmat(i,:)==0);
      grad = exp_num_trans(i,:) ./ denom;
      trim = find(transmat(i,:) <= exp(-(1/Z)*grad)); % eqn 32
      if ~isempty(trim)
	transmat(i,trim) = 0;
	trimmed_trans(i) = 1;
	disp(['trimming transitions ' num2str(i) ' -> ' num2str(trim)])
      end
    end
  end

  %obsmat = mk_stochastic(exp_num_emit);
  for i=1:Q
    obsmat(i,:) = entropic_map(exp_num_emit(i,:), Z);
    if trimobs & ~trimmed_obs(i) & (Z==1)
      denom = obsmat(i,:) + (obsmat(i,:)==0);
      grad = exp_num_emit(i,:) ./ denom;
      trim = find(obsmat(i,:) <= exp(-(1/Z)*grad)); % eqn 32
      if ~isempty(trim)
	obsmat(i,trim) = 0;
	trimmed_obs(i) = 1;
	disp(['trimming observations ' num2str(i) ' -> ' num2str(trim)])
      end
    end
  end

  if trimstates & (Z==1)
    prob_occ = sum(exp_num_emit, 2);
    trim = find((prob_occ < 1e-10) & ~trimmed_states);
    if ~isempty(trim)
      disp(['trimming states ' num2str(trim)])
      trimmed_states(trim) = 1;
      for i=trim(:)'
	transmat(:,i) = 0;
	transmat(i,:) = 0;
	obsmat(i,:) = 0;
      end
    end
  end

end

% compute log lik with the final param values
[loglik, exp_num_trans, exp_num_visits1, exp_num_emit, exp_num_visitsT] = ...
    compute_ess_dhmm(startprob, transmat, obsmat, data, dirichlet);
if verbose, fprintf(1, 'iteration %d, loglik = %7.4f, Z=%5.3f\n', num_iter, loglik, Z); end
LL = [LL loglik];

hmm.startprob = startprob;
hmm.endprob = endprob;
hmm.transmat = transmat;
hmm.obsmat = obsmat;