comparison toolboxes/FullBNT-1.0.7/bnt/learning/learn_params_dbn_em.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [bnet, LL, engine] = learn_params_dbn_em(engine, evidence, varargin)
+% LEARN_PARAMS_DBN Set the parameters in a DBN to their ML/MAP values using batch EM.
+% [bnet, LLtrace, engine] = learn_params_dbn_em(engine, data, ...)
+%
+% data{l}{i,t} = value of node i in slice t of time-series l, or [] if hidden.
+%   Suppose you have L time series, each of length T, in an O*T*L array D, 
+%   where O is the num of observed scalar nodes, and N is the total num nodes per slice.
+%   Then you can create data as follows, where onodes is the index of the observable nodes:
+%      data = cell(1,L);
+%      for l=1:L
+%        data{l} = cell(N, T);
+%        data{l}(onodes,:) = num2cell(D(:,:,l));
+%      end
+% Of course it is possible for different sets of nodes to be observed in
+% each slice/ sequence, and for each sequence to be a different length.
+%
+% LLtrace is the learning curve: the vector of log-likelihood scores at each iteration.
+%
+% Optional arguments [default]
+%
+% max_iter - specifies the maximum number of iterations [100]
+% thresh - specifies the thresold for stopping EM [1e-3]
+%   We stop when |f(t) - f(t-1)| / avg < threshold,
+%   where avg = (|f(t)| + |f(t-1)|)/2 and f is log lik.
+% verbose - display loglik at each iteration [1]
+% anneal - 1 means do deterministic annealing (only for entropic priors) [0]
+% anneal_rate - geometric cooling rate [0.8]
+% init_temp - initial annealing temperature [10]
+% final_temp - final annealing temperature [1e-3]
+%
+
+max_iter = 100;
+thresh = 1e-3;
+anneal = 0;
+anneal_rate = 0.8;
+init_temp = 10;
+final_temp = 1e-3;
+verbose = 1;
+
+for i=1:2:length(varargin)
+  switch varargin{i}
+   case 'max_iter', max_iter = varargin{i+1}; 
+   case 'thresh', thresh = varargin{i+1}; 
+   case 'anneal', anneal = varargin{i+1}; 
+   case 'anneal_rate', anneal_rate = varargin{i+1}; 
+   case 'init_temp', init_temp = varargin{i+1}; 
+   case 'final_temp', final_temp = varargin{i+1}; 
+   otherwise, error(['unrecognized argument' varargin{i}])
+  end
+end
+
+% take 1 EM step at each temperature value, then when temp=0, run to convergence
+% When using an entropic prior, Z = 1-T, so 
+% T=2 => Z=-1 (max entropy)
+% T=1 => Z=0 (max likelihood)
+% T=0 => Z=1 (min entropy / max structure)
+num_iter = 1;
+LL = [];
+if anneal
+  temperature = init_temp;
+  while temperature > final_temp
+    [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
+    if verbose
+      fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f, temp=%8.4f\n', ...
+	      num_iter, loglik, logpost, temperature);
+    end
+    num_iter = num_iter + 1;
+    LL = [LL loglik];
+    temperature = temperature * anneal_rate;
+  end
+  temperature = 0;
+  previous_loglik = loglik;
+  previous_logpost = logpost;
+else
+  temperature = 0;
+  previous_loglik = -inf;
+  previous_logpost = -inf;
+end
+
+converged = 0;
+while ~converged & (num_iter <= max_iter)
+  [engine, loglik, logpost] = EM_step(engine, evidence, temperature);
+  if verbose
+    %fprintf('EM iteration %d, loglik = %8.4f, logpost = %8.4f\n', ...
+    %	    num_iter, loglik, logpost);
+    fprintf('EM iteration %d, loglik = %8.4f\n', num_iter, loglik);
+  end
+  num_iter = num_iter + 1;
+  [converged, decreased] = em_converged(loglik, previous_loglik, thresh);
+  %[converged, decreased] = em_converged(logpost, previous_logpost, thresh);
+  previous_loglik = loglik;
+  previous_logpost = logpost;
+  LL = [LL loglik];
+end
+
+bnet = bnet_from_engine(engine);
+
+%%%%%%%%%
+
+function [engine, loglik, logpost] = EM_step(engine, cases, temp)
+
+bnet = bnet_from_engine(engine); % engine contains the old params that are used for the E step
+ss = length(bnet.intra);
+CPDs = bnet.CPD; % these are the new params that get maximized
+num_CPDs = length(CPDs);
+
+% log P(theta|D) = (log P(D|theta) + log P(theta)) - log(P(D))
+% where log P(D|theta) = sum_cases log P(case|theta)
+% and log P(theta) = sum_CPDs log P(CPD) - only count once even if tied!
+% logpost = log P(theta,D) (un-normalized)
+% This should be negative, and increase at every step.
+
+adjustable = zeros(1,num_CPDs);
+logprior = zeros(1, num_CPDs);
+for e=1:num_CPDs
+  adjustable(e) = adjustable_CPD(CPDs{e});
+end
+adj = find(adjustable);
+
+for e=adj(:)'
+  logprior(e) = log_prior(CPDs{e});
+  CPDs{e} = reset_ess(CPDs{e});
+end
+
+loglik = 0;
+for l=1:length(cases)
+  evidence = cases{l};
+  if ~iscell(evidence)
+    error('training data must be a cell array of cell arrays')
+  end
+  [engine, ll] = enter_evidence(engine, evidence);
+  assert(~isnan(ll))
+  loglik = loglik + ll;
+  T = size(evidence, 2);
+  
+  % We unroll ns etc because in update_ess, we refer to nodes by their unrolled number
+  % so that they extract evidence from the right place.
+  % (The CPD should really store its own version of ns and cnodes...)
+  ns = repmat(bnet.node_sizes_slice(:), [1 T]);
+  cnodes = unroll_set(bnet.cnodes_slice, ss, T);
+ 
+  %hidden_bitv = repmat(bnet.hidden_bitv(1:ss), [1 T]);
+  hidden_bitv = zeros(ss, T);
+  hidden_bitv(isemptycell(evidence))=1;
+  % hidden_bitv(i) = 1 means node i is hidden.
+  % We pass this in, rather than using isemptycell(evidence(dom)), because
+  % isemptycell is very slow.
+  
+  t = 1;
+  for i=1:ss
+    e = bnet.equiv_class(i,1);
+    if adjustable(e)
+      fmarg = marginal_family(engine, i, t);
+      CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:)); 
+    end
+  end
+  
+  for i=1:ss   
+    e = bnet.equiv_class(i,2);
+    if adjustable(e)
+      for t=2:T
+	fmarg = marginal_family(engine, i, t);
+	CPDs{e} = update_ess(CPDs{e}, fmarg, evidence, ns(:), cnodes(:), hidden_bitv(:));
+      end
+    end
+  end
+end
+
+logpost = loglik + sum(logprior(:));
+
+for e=adj(:)'
+  CPDs{e} = maximize_params(CPDs{e}, temp);
+end
+
+engine = update_engine(engine, CPDs);
+
+
+
+