enter_evidence.m

来自「Bayes网络工具箱」· M 代码 · 共 115 行

function [engine, loglik] = enter_evidence(engine, evidence, filter)
% ENTER_EVIDENCE Add the specified evidence to the network (hmm)
% [engine, loglik] = enter_evidence(engine, evidence, filter)
%
% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector)
% If filter = 1, we do filtering, otherwise smoothing (default).

if nargin < 3, filter = 0; end

n = length(engine.bnet.intra);
onodes = engine.onodes;
hnodes = mysetdiff(1:n, onodes);
vals = cat(1, evidence{onodes,:});
T = size(evidence, 2);
O = length(onodes); % num. observables per slice
vals = reshape(vals, [O T]);
% collapse the observables into a single value per slice
ns = engine.bnet.node_sizes(:);
os = ns(onodes);
data = subv2ind(os, vals'); % each row of vals is a vector of indices

obslik = mk_dhmm_obs_lik(data, engine.obsmat, engine.obsmat1);
[gamma, xi, loglik] = forwards_backwards(engine.prior, engine.transmat, obslik, filter);
  
%engine.loglik = loglik;
%engine.onodes = onodes;

% Wrap the posterior inside a potential, so it can be marginalized easily
engine.one_slice_marginal = cell(1,T);
engine.two_slice_marginal = cell(1,T);
ns(onodes) = 1;
ns(onodes+n) = 1;
ss = length(engine.bnet.intra);
for t=1:T
  dom = (1:n);
  engine.one_slice_marginal{t} = dpot(dom+(t-1)*ss, ns(dom), myreshape(gamma(:,t), ns(dom)));
end
for t=1:T-1
  dom = (1:(2*n));
  engine.two_slice_marginal{t} = dpot(dom+(t-1)*ss, ns(dom), myreshape(xi(:,:,t), ns(dom)));
end
     

%%%%%%%%%%%%%%

function B = mk_dhmm_obs_lik(data, obsmat, obsmat1)
% MK_DHMM_OBS_LIK  Make the observation likelihood vector for a discrete HMM.
% B = mk_dhmm_obs_lik(data, obsmat, obsmat1)
%
% Inputs:
% data(t) = y(t) = observation at time t
% obsmat(i,o) = Pr(Y(t)=o | Q(t)=i)
% obsmat1(i,o) = Pr(Y(1)=o | Q(1)=i). Defaults to obsmat if omitted.
%
% Output:
% B(i,t) = Pr(y(t) | Q(t)=i)

[Q O] = size(obsmat);
T = length(data);
B = zeros(Q,T);

t = 1;
B(:,t) = obsmat1(:, data(t));
for t=2:T
  B(:,t) = obsmat(:, data(t));
end               


%%%%%%%%%%%

function [gamma, xi, loglik] = forwards_backwards(prior, transmat, obslik, forwards_only)
%
% Outputs:
% gamma(i,t) = Pr(X(t)=i | O(1:tau))
% xi(i,j,t)  = Pr(X(t)=i, X(t+1)=j | O(1:tau)) t <= T-1
%
% where tau = t if forwards_only = 1, and tau = T otherwise.


T = size(obslik, 2);
Q = length(prior);

loglik = 0;
scale = ones(1,T);
alpha = zeros(Q,T); % alpha(i,t) = Pr(X(t)=i | O(1:t))
xi = zeros(Q,Q,T-1);

t = 1;
alpha(:,1) = prior(:) .* obslik(:,t);
[alpha(:,t), scale(t)] = normalise(alpha(:,t));
for t=2:T
  [alpha(:,t), scale(t)] = normalise((transmat' * alpha(:,t-1)) .* obslik(:,t));
  xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* transmat);
end
loglik = sum(log(scale));

if forwards_only
  gamma = alpha;
  return;
end

beta = zeros(Q,T); % beta(i,t)  = Pr(O(t+1:T) | X(t)=i)
gamma = zeros(Q,T);
beta(:,T) = ones(Q,1);
gamma(:,T) = normalise(alpha(:,T) .* beta(:,T));
t=T;
for t=T-1:-1:1
  b = beta(:,t+1) .* obslik(:,t+1); 
  beta(:,t) = normalise((transmat * b));
  gamma(:,t) = normalise(alpha(:,t) .* beta(:,t));
  xi(:,:,t) = normalise((transmat .* (alpha(:,t) * b')));
end