comparison toolboxes/FullBNT-1.0.7/bnt/examples/dynamic/fhmm_infer.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [loglik, gamma] = fhmm_infer(inter, CPTs_slice1, CPTs, obsmat, node_sizes)
+% FHMM_INFER Exact inference for a factorial HMM.
+% [loglik, gamma] = fhmm_infer(inter, CPTs_slice1, CPTs, obsmat, node_sizes)
+%
+% Inputs:
+% inter - the inter-slice adjacency matrix
+% CPTs_slice1{s}(j) = Pr(Q(s,1) = j) where Q(s,t) = hidden node s in slice t
+% CPT{s}(i1, i2, ..., j) = Pr(Q(s,t) = j | Pa(s,t-1) = i1, i2, ...),
+% obsmat(i,t) = Pr(y(t) | Q(t)=i)
+% node_sizes is a vector with the cardinality of the hidden nodes
+%
+% Outputs:
+% gamma(i,t) = Pr(X(t)=i | O(1:T)) as in an HMM,
+% except that i is interpreted as an M digit, base-K number (if there are M chains each of cardinality K).
+%
+%
+% For M chains each of cardinality K, the frontiers  (i.e., cliques)
+% contain M+1 nodes, and it takes M steps to advance the frontier by one time step,
+% so the run time is O(T M K^(M+1)).
+% An HMM takes O(T S^2) where S is the size of the state space.
+% Collapsing the FHMM to an HMM results in S = K^M.
+% For details, see
+%   "The Factored Frontier Algorithm for Approximate Inference in DBNs",
+%    Kevin Murphy and Yair Weiss, submitted to NIPS 2000.
+%
+% The frontier algorithm makes the following topological assumptions:
+% 
+%  - All nodes are persistent (connect to the next slice)
+%  - No connections within a timeslice
+%  - There is a single observation variable, which depends on all the hidden nodes
+%  - Each node can have several parents in the previous time slice (generalizes a FHMM slightly)
+%
+
+% The forwards pass of the frontier algorithm can be explained with the following example.
+% Suppose we have 3 hidden nodes per slice, A, B, C.
+% The goal is to compute alpha(j, t) = Pr( (A_t,B_t,C_t)=j | Y(1:t))
+% We move alpha from t to t+1 one node at a time, as follows.
+% We define the following quantities:
+% s([a1 b1 c1], 1) = Prob(A(t)=a1, B(t)=b1, C(t)=c1 | Y(1:t)) = alpha(j, t)
+% s([a2 b1 c1], 2) = Prob(A(t+1)=a2, B(t)=b1, C(t)=c1 | Y(1:t))
+% s([a2 b2 c1], 3) = Prob(A(t+1)=a2, B(t+1)=b2, C(t)=c1 | Y(1:t))
+% s([a2 b2 c2], 4) = Prob(A(t+1)=a2, B(t+1)=b2, C(t+1)=c2 | Y(1:t))
+% s([a2 b2 c2], 5) = Prob(A(t+1)=a2, B(t+1)=b2, C(t+1)=c2 | Y(1:t+1)) = alpha(j, t+1)
+%
+% These can be computed recursively as follows:
+%
+% s([a2 b1 c1], 2) = sum_{a1} P(a2|a1) s([a1 b1 c1], 1)
+% s([a2 b2 c1], 3) = sum_{b1} P(b2|b1) s([a2 b1 c1], 2)
+% s([a2 b2 c2], 4) = sum_{c1} P(c2|c1) s([a2 b2 c1], 1)
+% s([a2 b2 c2], 5) = normalise( s([a2 b2 c2], 4) .* P(Y(t+1)|a2,b2,c2)
+
+
+[kk,ll,mm] = make_frontier_indices(inter, node_sizes); % can pass in as args
+
+scaled = 1;
+
+M = length(node_sizes);
+S = prod(node_sizes);
+T = size(obsmat, 2);
+
+alpha = zeros(S, T);
+beta = zeros(S, T);
+gamma = zeros(S, T);
+scale = zeros(1,T);
+tiny = exp(-700);
+
+
+alpha(:,1) = make_prior_from_CPTs(CPTs_slice1, node_sizes);
+alpha(:,1) = alpha(:,1) .* obsmat(:, 1);
+
+if scaled
+  s = sum(alpha(:,1));
+  if s==0, s = s + tiny; end
+  scale(1) = 1/s;
+else
+  scale(1) = 1;
+end
+alpha(:,1) = alpha(:,1) * scale(1);
+
+%a = zeros(S, M+1);
+%b = zeros(S, M+1);
+anew = zeros(S,1);
+aold = zeros(S,1);
+bnew = zeros(S,1);
+bold = zeros(S,1);
+
+for t=2:T
+  %a(:,1) = alpha(:,t-1);
+  aold =  alpha(:,t-1);
+  
+  c = 1;
+  for i=1:M
+    ns = node_sizes(i);
+    cpt = CPTs{i};
+    for j=1:S
+      s = 0;
+      for xx=1:ns
+	%k = kk(xx,j,i);
+	%l = ll(xx,j,i);
+	k = kk(c);
+	l = ll(c);
+	c = c + 1;
+	% s = s + a(k,i) * CPTs{i}(l);
+	s = s + aold(k) * cpt(l);
+      end
+      %a(j,i+1) = s;
+      anew(j) = s;
+    end
+    aold = anew;
+  end
+  
+  %alpha(:,t) = a(:,M+1) .* obsmat(:, obs(t));
+  alpha(:,t) = anew .* obsmat(:, t);
+
+  if scaled
+    s = sum(alpha(:,t));
+    if s==0, s = s + tiny; end
+    scale(t) = 1/s;
+  else
+    scale(t) = 1;
+  end
+  alpha(:,t) = alpha(:,t) * scale(t);
+
+end
+
+
+beta(:,T) = ones(S,1) * scale(T);
+for t=T-1:-1:1
+  %b(:,1) = beta(:,t+1) .* obsmat(:, obs(t+1));
+  bold = beta(:,t+1) .* obsmat(:, t+1);
+
+  c = 1;
+  for i=1:M
+    ns = node_sizes(i);
+    cpt = CPTs{i};
+    for j=1:S
+      s = 0;
+      for xx=1:ns
+	%k = kk(xx,j,i);
+	%m = mm(xx,j,i);
+	k = kk(c);
+	m = mm(c);
+	c = c + 1;
+	% s = s + b(k,i) * CPTs{i}(m);
+	s = s + bold(k) * cpt(m);
+      end
+      %b(j,i+1) = s;
+      bnew(j) = s;
+    end
+    bold = bnew;
+  end
+  % beta(:,t) = b(:,M+1) * scale(t);
+  beta(:,t) = bnew * scale(t);
+end
+
+
+if scaled
+  loglik = -sum(log(scale)); % scale(i) is finite
+else
+  lik = alpha(:,1)' * beta(:,1);
+  loglik = log(lik+tiny);
+end
+
+for t=1:T
+  gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
+end
+
+%%%%%%%%%%%
+
+function [kk,ll,mm] = make_frontier_indices(inter, node_sizes)
+%
+% Precompute indices for use in the frontier algorithm.
+% These only depend on the topology, not the parameters or data.
+% Hence we can compute them outside of fhmm_infer.
+% This saves a lot of run-time computation.
+
+M = length(node_sizes);
+S = prod(node_sizes);
+
+mns = max(node_sizes);
+kk = zeros(mns, S, M);
+ll = zeros(mns, S, M);
+mm = zeros(mns, S, M);
+
+for i=1:M
+  for j=1:S
+    u = ind2subv(node_sizes, j);
+    x = u(i);
+    for xx=1:node_sizes(i)
+      uu = u;
+      uu(i) = xx;
+      k = subv2ind(node_sizes, uu);
+      kk(xx,j,i) = k;
+      ps = find(inter(:,i)==1);
+      ps = ps(:)';
+      l = subv2ind(node_sizes([ps i]), [uu(ps) x]); % sum over parent
+      ll(xx,j,i) = l;
+      m = subv2ind(node_sizes([ps i]), [u(ps) xx]); % sum over child
+      mm(xx,j,i) = m;
+    end
+  end
+end
+
+%%%%%%%%%
+
+function prior=make_prior_from_CPTs(indiv_priors, node_sizes)
+%
+% composite_prior=make_prior(individual_priors, node_sizes)
+% Make the prior for the first node in a Markov chain
+% from the priors on each node in the equivalent DBN.
+% prior{i}(j) = Pr(X_i=j), where X_i is the i'th node in slice 1.
+% composite_prior(i) = Pr(slice1 = i).
+
+n = length(indiv_priors);
+S = prod(node_sizes);
+prior = zeros(S,1);
+for i=1:S
+  vi = ind2subv(node_sizes, i);
+  p = 1;
+  for k=1:n
+    p = p * indiv_priors{k}(vi(k));
+  end
+  prior(i) = p;
+end
+
+
+
+%%%%%%%%%%%
+
+function [loglik, alpha, beta] = FHMM_slow(inter, CPTs_slice1, CPTs, obsmat, node_sizes, data)
+% 
+% Same as the above, except we don't use the optimization of computing the indices outside the loop.
+
+
+scaled = 1;
+
+M = length(node_sizes);
+S = prod(node_sizes);
+[numex T] = size(data);
+
+obs = data;
+
+alpha = zeros(S, T);
+beta = zeros(S, T);
+a = zeros(S, M+1);
+b = zeros(S, M+1);
+scale = zeros(1,T);
+
+alpha(:,1) = make_prior_from_CPTs(CPTs_slice1, node_sizes);
+alpha(:,1) = alpha(:,1) .* obsmat(:, obs(1));
+if scaled
+  s = sum(alpha(:,1));
+  if s==0, s = s + tiny; end
+  scale(1) = 1/s;
+else
+  scale(1) = 1;
+end
+alpha(:,1) = alpha(:,1) * scale(1);
+
+for t=2:T
+  fprintf(1, 't %d\n', t);
+  a(:,1) = alpha(:,t-1);
+  for i=1:M
+    for j=1:S
+      u = ind2subv(node_sizes, j);
+      xnew = u(i);
+      s = 0;
+      for xold=1:node_sizes(i)
+	uold = u;
+	uold(i) = xold;
+	k = subv2ind(node_sizes, uold);
+	ps = find(inter(:,i)==1);
+	ps = ps(:)';
+	l = subv2ind(node_sizes([ps i]), [uold(ps) xnew]);
+	s = s + a(k,i) * CPTs{i}(l);
+      end
+      a(j,i+1) = s;
+    end
+  end
+  alpha(:,t) = a(:,M+1) .* obsmat(:, obs(t));
+
+  if scaled
+    s = sum(alpha(:,t));
+    if s==0, s = s + tiny; end
+    scale(t) = 1/s;
+  else
+    scale(t) = 1;
+  end
+  alpha(:,t) = alpha(:,t) * scale(t);
+
+end
+
+
+beta(:,T) = ones(S,1) * scale(T);
+for t=T-1:-1:1
+  fprintf(1, 't %d\n', t);
+  b(:,1) = beta(:,t+1) .* obsmat(:, obs(t+1));
+  for i=1:M
+    for j=1:S
+      u = ind2subv(node_sizes, j);
+      xold = u(i);
+      s = 0;
+      for xnew=1:node_sizes(i)
+	unew = u;
+	unew(i) = xnew;
+	k = subv2ind(node_sizes, unew);
+	ps = find(inter(:,i)==1);
+	ps = ps(:)';
+	l = subv2ind(node_sizes([ps i]), [u(ps) xnew]);
+	s = s + b(k,i) * CPTs{i}(l);
+      end
+      b(j,i+1) = s;
+    end
+  end
+  beta(:,t) = b(:,M+1) * scale(t);
+end
+
+
+if scaled
+  loglik = -sum(log(scale)); % scale(i) is finite
+else
+  lik = alpha(:,1)' * beta(:,1);
+  loglik = log(lik+tiny);
+end