Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [alpha, beta, gamma, loglik, xi, gamma2] = fwdback_twoslice(engine, init_state_distrib, transmat, obslik, varargin)
wolffd@0 2 % FWDBACK Compute the posterior probs. in an HMM using the forwards backwards algo.
wolffd@0 3 %
wolffd@0 4 % [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(init_state_distrib, transmat, obslik, ...)
wolffd@0 5 %
wolffd@0 6 % Notation:
wolffd@0 7 % Y(t) = observation, Q(t) = hidden state, M(t) = mixture variable (for MOG outputs)
wolffd@0 8 % A(t) = discrete input (action) (for POMDP models)
wolffd@0 9 %
wolffd@0 10 % INPUT:
wolffd@0 11 % init_state_distrib(i) = Pr(Q(1) = i)
wolffd@0 12 % transmat(i,j) = Pr(Q(t) = j | Q(t-1)=i)
wolffd@0 13 % or transmat{a}(i,j) = Pr(Q(t) = j | Q(t-1)=i, A(t-1)=a) if there are discrete inputs
wolffd@0 14 % obslik(i,t) = Pr(Y(t)| Q(t)=i)
wolffd@0 15 % (Compute obslik using eval_pdf_xxx on your data sequence first.)
wolffd@0 16 %
wolffd@0 17 % Optional parameters may be passed as 'param_name', param_value pairs.
wolffd@0 18 % Parameter names are shown below; default values in [] - if none, argument is mandatory.
wolffd@0 19 %
wolffd@0 20 % For HMMs with MOG outputs: if you want to compute gamma2, you must specify
wolffd@0 21 % 'obslik2' - obslik(i,j,t) = Pr(Y(t)| Q(t)=i,M(t)=j) []
wolffd@0 22 % 'mixmat' - mixmat(i,j) = Pr(M(t) = j | Q(t)=i) []
wolffd@0 23 %
wolffd@0 24 % For HMMs with discrete inputs:
wolffd@0 25 % 'act' - act(t) = action performed at step t
wolffd@0 26 %
wolffd@0 27 % Optional arguments:
wolffd@0 28 % 'fwd_only' - if 1, only do a forwards pass and set beta=[], gamma2=[] [0]
wolffd@0 29 % 'scaled' - if 1, normalize alphas and betas to prevent underflow [1]
wolffd@0 30 % 'maximize' - if 1, use max-product instead of sum-product [0]
wolffd@0 31 %
wolffd@0 32 % OUTPUTS:
wolffd@0 33 % alpha(i,t) = p(Q(t)=i | y(1:t)) (or p(Q(t)=i, y(1:t)) if scaled=0)
wolffd@0 34 % beta(i,t) = p(y(t+1:T) | Q(t)=i)*p(y(t+1:T)|y(1:t)) (or p(y(t+1:T) | Q(t)=i) if scaled=0)
wolffd@0 35 % gamma(i,t) = p(Q(t)=i | y(1:T))
wolffd@0 36 % loglik = log p(y(1:T))
wolffd@0 37 % xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:T))
wolffd@0 38 % gamma2(j,k,t) = p(Q(t)=j, M(t)=k | y(1:T)) (only for MOG outputs)
wolffd@0 39 %
wolffd@0 40 % If fwd_only = 1, these become
wolffd@0 41 % alpha(i,t) = p(Q(t)=i | y(1:t))
wolffd@0 42 % beta = []
wolffd@0 43 % gamma(i,t) = p(Q(t)=i | y(1:t))
wolffd@0 44 % xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:t))
wolffd@0 45 % gamma2 = []
wolffd@0 46 %
wolffd@0 47 % Note: we only compute xi if it is requested as a return argument, since it can be very large.
wolffd@0 48 % Similarly, we only compute gamma2 on request (and if using MOG outputs).
wolffd@0 49 %
wolffd@0 50 % Examples:
wolffd@0 51 %
wolffd@0 52 % [alpha, beta, gamma, loglik] = fwdback(pi, A, multinomial_prob(sequence, B));
wolffd@0 53 %
wolffd@0 54 % [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat);
wolffd@0 55 % [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(pi, A, B, 'obslik2', B2, 'mixmat', mixmat);
wolffd@0 56
wolffd@0 57
wolffd@0 58 if nargout >= 5, compute_xi = 1; else compute_xi = 0; end
wolffd@0 59 if nargout >= 6, compute_gamma2 = 1; else compute_gamma2 = 0; end
wolffd@0 60
wolffd@0 61 [obslik2, mixmat, fwd_only, scaled, act, maximize, compute_xi, compute_gamma2] = process_options(varargin, 'obslik2', [], 'mixmat', [], 'fwd_only', 0, 'scaled', 1, 'act', [], 'maximize', 0, 'compute_xi', compute_xi, 'compute_gamma2', compute_gamma2);
wolffd@0 62
wolffd@0 63
wolffd@0 64 [Q T] = size(obslik);
wolffd@0 65
wolffd@0 66 if isempty(obslik2)
wolffd@0 67 compute_gamma2 = 0;
wolffd@0 68 end
wolffd@0 69
wolffd@0 70 if isempty(act)
wolffd@0 71 act = ones(1,T);
wolffd@0 72 transmat = { transmat } ;
wolffd@0 73 end
wolffd@0 74
wolffd@0 75 scale = ones(1,T);
wolffd@0 76
wolffd@0 77 % scale(t) = Pr(O(t) | O(1:t-1)) = 1/c(t) as defined by Rabiner (1989).
wolffd@0 78 % Hence prod_t scale(t) = Pr(O(1)) Pr(O(2)|O(1)) Pr(O(3) | O(1:2)) = Pr(O(1), ... ,O(T))
wolffd@0 79 % or log P = sum_t log scale(t).
wolffd@0 80 % Rabiner suggests multiplying beta(t) by scale(t), but we can instead
wolffd@0 81 % normalise beta(t) - the constants will cancel when we compute gamma.
wolffd@0 82
wolffd@0 83 loglik = 0;
wolffd@0 84
wolffd@0 85 alpha = zeros(Q,T);
wolffd@0 86 gamma = zeros(Q,T);
wolffd@0 87 if compute_xi
wolffd@0 88 xi = zeros(Q,Q,T-1);
wolffd@0 89 else
wolffd@0 90 xi = [];
wolffd@0 91 end
wolffd@0 92
wolffd@0 93
wolffd@0 94 %%%%%%%%% Forwards %%%%%%%%%%
wolffd@0 95
wolffd@0 96 t = 1;
wolffd@0 97 alpha(:,1) = init_state_distrib(:) .* obslik(:,t);
wolffd@0 98 if scaled
wolffd@0 99 %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
wolffd@0 100 [alpha(:,t), scale(t)] = normalise(alpha(:,t));
wolffd@0 101 end
wolffd@0 102 if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
wolffd@0 103 for t=2:T
wolffd@0 104 %trans = transmat(:,:,act(t-1))';
wolffd@0 105 trans = transmat{act(t-1)};
wolffd@0 106 if maximize
wolffd@0 107 m = max_mult(trans', alpha(:,t-1));
wolffd@0 108 %A = repmat(alpha(:,t-1), [1 Q]);
wolffd@0 109 %m = max(trans .* A, [], 1);
wolffd@0 110 else
wolffd@0 111 m = trans' * alpha(:,t-1);
wolffd@0 112 end
wolffd@0 113 alpha(:,t) = m(:) .* obslik(:,t);
wolffd@0 114 if scaled
wolffd@0 115 %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t));
wolffd@0 116 [alpha(:,t), scale(t)] = normalise(alpha(:,t));
wolffd@0 117 end
wolffd@0 118 if compute_xi & fwd_only % useful for online EM
wolffd@0 119 %xi(:,:,t-1) = normaliseC((alpha(:,t-1) * obslik(:,t)') .* trans);
wolffd@0 120 xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* trans);
wolffd@0 121 end
wolffd@0 122 if scaled, assert(approxeq(sum(alpha(:,t)),1)), end
wolffd@0 123 end
wolffd@0 124 if scaled
wolffd@0 125 if any(scale==0)
wolffd@0 126 loglik = -inf;
wolffd@0 127 else
wolffd@0 128 loglik = sum(log(scale));
wolffd@0 129 end
wolffd@0 130 else
wolffd@0 131 loglik = log(sum(alpha(:,T)));
wolffd@0 132 end
wolffd@0 133
wolffd@0 134 if fwd_only
wolffd@0 135 gamma = alpha;
wolffd@0 136 beta = [];
wolffd@0 137 gamma2 = [];
wolffd@0 138 return;
wolffd@0 139 end
wolffd@0 140
wolffd@0 141
wolffd@0 142 %%%%%%%%% Backwards %%%%%%%%%%
wolffd@0 143
wolffd@0 144 beta = zeros(Q,T);
wolffd@0 145 if compute_gamma2
wolffd@0 146 M = size(mixmat, 2);
wolffd@0 147 gamma2 = zeros(Q,M,T);
wolffd@0 148 else
wolffd@0 149 gamma2 = [];
wolffd@0 150 end
wolffd@0 151
wolffd@0 152 beta(:,T) = ones(Q,1);
wolffd@0 153 %gamma(:,T) = normaliseC(alpha(:,T) .* beta(:,T));
wolffd@0 154 gamma(:,T) = normalise(alpha(:,T) .* beta(:,T));
wolffd@0 155 t=T;
wolffd@0 156 if compute_gamma2
wolffd@0 157 denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
wolffd@0 158 gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
wolffd@0 159 %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); % wrong!
wolffd@0 160 end
wolffd@0 161 for t=T-1:-1:1
wolffd@0 162 b = beta(:,t+1) .* obslik(:,t+1);
wolffd@0 163 %trans = transmat(:,:,act(t));
wolffd@0 164 trans = transmat{act(t)};
wolffd@0 165 if maximize
wolffd@0 166 B = repmat(b(:)', Q, 1);
wolffd@0 167 beta(:,t) = max(trans .* B, [], 2);
wolffd@0 168 else
wolffd@0 169 beta(:,t) = trans * b;
wolffd@0 170 end
wolffd@0 171 if scaled
wolffd@0 172 %beta(:,t) = normaliseC(beta(:,t));
wolffd@0 173 beta(:,t) = normalise(beta(:,t));
wolffd@0 174 end
wolffd@0 175 %gamma(:,t) = normaliseC(alpha(:,t) .* beta(:,t));
wolffd@0 176 gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
wolffd@0 177 if compute_xi
wolffd@0 178 %xi(:,:,t) = normaliseC((trans .* (alpha(:,t) * b')));
wolffd@0 179 xi(:,:,t) = normalise((trans .* (alpha(:,t) * b')));
wolffd@0 180 %xi(:,:,t) = (trans .* (alpha(:,t) * b'));
wolffd@0 181 end
wolffd@0 182 if compute_gamma2
wolffd@0 183 denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing
wolffd@0 184 gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]);
wolffd@0 185 %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]));
wolffd@0 186 end
wolffd@0 187 end
wolffd@0 188
wolffd@0 189
wolffd@0 190 % We now explain the equation for gamma2
wolffd@0 191 % Let zt=y(1:t-1,t+1:T) be all observations except y(t)
wolffd@0 192 % gamma2(Q,M,t) = P(Qt,Mt|yt,zt) = P(yt|Qt,Mt,zt) P(Qt,Mt|zt) / P(yt|zt)
wolffd@0 193 % = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|zt) / P(yt|zt)
wolffd@0 194 % Now gamma(Q,t) = P(Qt|yt,zt) = P(yt|Qt) P(Qt|zt) / P(yt|zt)
wolffd@0 195 % hence
wolffd@0 196 % P(Qt,Mt|yt,zt) = P(yt|Qt,Mt) P(Mt|Qt) [P(Qt|yt,zt) P(yt|zt) / P(yt|Qt)] / P(yt|zt)
wolffd@0 197 % = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|yt,zt) / P(yt|Qt)
wolffd@0 198 %