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