diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine')
21 files changed, 620 insertions, 0 deletions
diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries new file mode 100644 index 00000000..e1ab6d00 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Entries @@ -0,0 +1,9 @@ +/enter_evidence.m/1.2/Sat Sep 17 17:00:30 2005// +/find_mpe.m/1.1.1.1/Thu Jun 20 00:18:24 2002// +/fwdback_twoslice.m/1.1/Sat Nov 26 01:24:09 2005// +/hmm_inf_engine.m/1.1.1.1/Thu Nov 14 20:05:36 2002// +/marginal_family.m/1.1.1.1/Thu Nov 14 20:05:36 2002// +/marginal_nodes.m/1.1.1.1/Thu Nov 14 20:03:28 2002// +/update_engine.m/1.1.1.1/Wed May 29 15:59:56 2002// +D/Old//// +D/private//// diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository new file mode 100644 index 00000000..b7392efa --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@hmm_inf_engine diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries new file mode 100644 index 00000000..528b5843 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Entries @@ -0,0 +1,4 @@ +/dhmm_inf_engine.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_family.m/1.1.1.1/Wed May 29 15:59:56 2002// +/marginal_nodes.m/1.1.1.1/Wed May 29 15:59:56 2002// +D diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository new file mode 100644 index 00000000..f82b2bea --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@hmm_inf_engine/Old diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m new file mode 100644 index 00000000..2b0c8810 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/dhmm_inf_engine.m @@ -0,0 +1,34 @@ +function engine = dhmm_inf_engine(bnet, onodes) +% DHMM_INF_ENGINE Inference engine for discrete DBNs which uses the forwards-backwards algorithm. +% engine = dhmm_inf_engine(bnet, onodes) +% +% 'onodes' specifies which nodes are observed; these must be leaves, and can be discrete or continuous. +% The remaining nodes are all hidden, and must be discrete. +% The DBN is converted to an HMM, with a single meganode, but which may have factored obs. + +ss = length(bnet.intra); +hnodes = mysetdiff(1:ss, onodes); +evidence = cell(ss, 2); +ns = bnet.node_sizes; +Q = prod(ns(hnodes)); +tmp = dpot_to_table(compute_joint_pot(bnet, hnodes, evidence)); +engine.startprob = reshape(tmp, Q, 1); +tmp = dpot_to_table(compute_joint_pot(bnet, [hnodes hnodes+ss], evidence)); +engine.transprob = mk_stochastic(reshape(tmp, Q, Q)); +engine.obsprob = cell(1, length(onodes)); +for i=1:length(onodes) + tmp = dpot_to_table(compute_joint_pot(bnet, [hnodes onodes(i)], evidence)); + O = ns(onodes(i)); + engine.obsprob{i} = mk_stochastic(reshape(tmp, Q, O)); +end + +% This is where we will store the results between enter_evidence and marginal_nodes +engine.gamma = []; +engine.xi = []; + +engine.onodes = onodes; +engine.hnodes = hnodes; +engine.maximize = []; + +engine = class(engine, 'dhmm_inf_engine', inf_engine(bnet)); + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m new file mode 100644 index 00000000..681e1591 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_family.m @@ -0,0 +1,31 @@ +function marginal = marginal_family(engine, i, t, add_ev) +% MARGINAL_FAMILY Compute the marginal on the specified family (hmm) +% marginal = marginal_nodes(engine, i, t, add_ev) +% + +if nargin < 3, t = 1; end +if nargin < 4, add_ev = 0; end + +bnet = bnet_from_engine(engine); +ss = length(bnet.intra); +if t==1 + fam = family(bnet.dag, i); + bigpot = engine.one_slice_marginal{t}; + nodes = fam; +else + fam = family(bnet.dag, i+ss); + if any(fam <= ss) % family spans 2 slices + bigpot = engine.two_slice_marginal{t-1}; % t-1 and t + nodes = fam + (t-2)*ss; + else + bigpot = engine.one_slice_marginal{t}; + nodes = fam-ss + (t-1)*ss; + end +end + +marginal = pot_to_marginal(marginalize_pot(bigpot, nodes, engine.maximize)); + +if add_ev + marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m new file mode 100644 index 00000000..4b8d6008 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/Old/marginal_nodes.m @@ -0,0 +1,30 @@ +function marginal = marginal_nodes(engine, nodes, t, add_ev) +% MARGINAL_NODES Compute the marginal on the specified query nodes (hmm) +% marginal = marginal_nodes(engine, nodes, t, add_ev) +% +% 't' specifies the time slice of the earliest node in 'nodes'. +% 'nodes' cannot span more than 2 time slices. +% +% Example: +% Consider a DBN with 2 nodes per slice. +% Then t=2, nodes=[1 3] refers to node 1 in slice 2 and node 1 in slice 3, +% i.e., nodes 3 and 5 in the unrolled network, + +if nargin < 3, t = 1; end +if nargin < 4, add_ev = 0; end + +bnet = bnet_from_engine(engine); +ss = length(bnet.intra); +if all(nodes <= ss) + bigpot = engine.one_slice_marginal{t}; +else + bigpot = engine.two_slice_marginal{t}; +end + +nodes = nodes + (t-1)*ss; +marginal = pot_to_marginal(marginalize_pot(bigpot, nodes, engine.maximize)); + +if add_ev + marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m new file mode 100644 index 00000000..8af97edf --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/enter_evidence.m @@ -0,0 +1,64 @@ +function [engine, loglik] = enter_evidence(engine, evidence, varargin) +% ENTER_EVIDENCE Add the specified evidence to the network (hmm) +% [engine, loglik] = enter_evidence(engine, evidence, ...) +% +% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector) +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% maximize - if 1, does max-product (not yet supported), else sum-product [0] +% filter - if 1, does filtering, else smoothing [0] +% oneslice - 1 means only compute marginals on nodes within a single slice [0] +% +% e.g., engine = enter_evidence(engine, ev, 'maximize', 1) + +maximize = 0; +filter = 0; +oneslice = 0; + +% parse optional params +args = varargin; +nargs = length(args); +if nargs > 0 + for i=1:2:nargs + switch args{i}, + case 'maximize', maximize = args{i+1}; + case 'filter', filter = args{i+1}; + case 'oneslice', oneslice = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +[ss T] = size(evidence); +engine.maximize = maximize; +engine.evidence = evidence; +bnet = bnet_from_engine(engine); +engine.node_sizes = repmat(bnet.node_sizes_slice(:), [1 T]); + +obs_bitv = ~isemptycell(evidence(:)); +bitv = reshape(obs_bitv, ss, T); +for t=1:T + onodes = find(bitv(:,t)); + if ~isequal(onodes, bnet.observed(:)) + error(['dbn was created assuming observed nodes per slice were '... + num2str(bnet.observed(:)') ' but the evidence in slice ' num2str(t) ... + ' has observed nodes ' num2str(onodes(:)')]); + end +end + +obslik = mk_hmm_obs_lik_matrix(engine, evidence); + +%[alpha, beta, gamma, loglik, xi] = fwdback(engine.startprob, engine.transprob, obslik, ... +[alpha, beta, gamma, loglik, xi] = fwdback_twoslice(engine, engine.startprob,... + engine.transprob, obslik, ... + 'maximize', maximize, 'fwd_only', filter, ... + 'compute_xi', ~oneslice); + +engine.one_slice_marginal = gamma; % gamma(:,t) for t=1:T +if ~oneslice + Q = size(gamma,1); + engine.two_slice_marginal = reshape(xi, [Q*Q T-1]); % xi(:,t) for t=1:T-1 +end diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m new file mode 100644 index 00000000..ba2cba74 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/find_mpe.m @@ -0,0 +1,16 @@ +function mpe = find_mpe(engine, evidence) +% FIND_MPE Find the most probable explanation (Viterbi) +% mpe = enter_evidence(engine, evidence, ...) +% +% evidence{i,t} = [] if if X(i,t) is hidden, and otherwise contains its observed value (scalar or column vector) +% + +obslik = mk_hmm_obs_lik_matrix(engine, evidence); +path = viterbi_path(engine.startprob, engine.transprob, obslik); +bnet = bnet_from_engine(engine); +ns = bnet.node_sizes_slice; +ns(bnet.observed) = 1; +ass = ind2subv(ns, path); +mpe = num2cell(ass'); +mpe(bnet.observed,:) = evidence(bnet.observed,:); + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m new file mode 100644 index 00000000..0565e727 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/fwdback_twoslice.m @@ -0,0 +1,198 @@ +function [alpha, beta, gamma, loglik, xi, gamma2] = fwdback_twoslice(engine, init_state_distrib, transmat, obslik, varargin) +% FWDBACK Compute the posterior probs. in an HMM using the forwards backwards algo. +% +% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(init_state_distrib, transmat, obslik, ...) +% +% Notation: +% Y(t) = observation, Q(t) = hidden state, M(t) = mixture variable (for MOG outputs) +% A(t) = discrete input (action) (for POMDP models) +% +% INPUT: +% init_state_distrib(i) = Pr(Q(1) = i) +% transmat(i,j) = Pr(Q(t) = j | Q(t-1)=i) +% or transmat{a}(i,j) = Pr(Q(t) = j | Q(t-1)=i, A(t-1)=a) if there are discrete inputs +% obslik(i,t) = Pr(Y(t)| Q(t)=i) +% (Compute obslik using eval_pdf_xxx on your data sequence first.) +% +% Optional parameters may be passed as 'param_name', param_value pairs. +% Parameter names are shown below; default values in [] - if none, argument is mandatory. +% +% For HMMs with MOG outputs: if you want to compute gamma2, you must specify +% 'obslik2' - obslik(i,j,t) = Pr(Y(t)| Q(t)=i,M(t)=j) [] +% 'mixmat' - mixmat(i,j) = Pr(M(t) = j | Q(t)=i) [] +% +% For HMMs with discrete inputs: +% 'act' - act(t) = action performed at step t +% +% Optional arguments: +% 'fwd_only' - if 1, only do a forwards pass and set beta=[], gamma2=[] [0] +% 'scaled' - if 1, normalize alphas and betas to prevent underflow [1] +% 'maximize' - if 1, use max-product instead of sum-product [0] +% +% OUTPUTS: +% alpha(i,t) = p(Q(t)=i | y(1:t)) (or p(Q(t)=i, y(1:t)) if scaled=0) +% 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) +% gamma(i,t) = p(Q(t)=i | y(1:T)) +% loglik = log p(y(1:T)) +% xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:T)) +% gamma2(j,k,t) = p(Q(t)=j, M(t)=k | y(1:T)) (only for MOG outputs) +% +% If fwd_only = 1, these become +% alpha(i,t) = p(Q(t)=i | y(1:t)) +% beta = [] +% gamma(i,t) = p(Q(t)=i | y(1:t)) +% xi(i,j,t-1) = p(Q(t-1)=i, Q(t)=j | y(1:t)) +% gamma2 = [] +% +% Note: we only compute xi if it is requested as a return argument, since it can be very large. +% Similarly, we only compute gamma2 on request (and if using MOG outputs). +% +% Examples: +% +% [alpha, beta, gamma, loglik] = fwdback(pi, A, multinomial_prob(sequence, B)); +% +% [B, B2] = mixgauss_prob(data, mu, Sigma, mixmat); +% [alpha, beta, gamma, loglik, xi, gamma2] = fwdback(pi, A, B, 'obslik2', B2, 'mixmat', mixmat); + + +if nargout >= 5, compute_xi = 1; else compute_xi = 0; end +if nargout >= 6, compute_gamma2 = 1; else compute_gamma2 = 0; end + +[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); + + +[Q T] = size(obslik); + +if isempty(obslik2) + compute_gamma2 = 0; +end + +if isempty(act) + act = ones(1,T); + transmat = { transmat } ; +end + +scale = ones(1,T); + +% scale(t) = Pr(O(t) | O(1:t-1)) = 1/c(t) as defined by Rabiner (1989). +% Hence prod_t scale(t) = Pr(O(1)) Pr(O(2)|O(1)) Pr(O(3) | O(1:2)) = Pr(O(1), ... ,O(T)) +% or log P = sum_t log scale(t). +% Rabiner suggests multiplying beta(t) by scale(t), but we can instead +% normalise beta(t) - the constants will cancel when we compute gamma. + +loglik = 0; + +alpha = zeros(Q,T); +gamma = zeros(Q,T); +if compute_xi + xi = zeros(Q,Q,T-1); +else + xi = []; +end + + +%%%%%%%%% Forwards %%%%%%%%%% + +t = 1; +alpha(:,1) = init_state_distrib(:) .* obslik(:,t); +if scaled + %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t)); + [alpha(:,t), scale(t)] = normalise(alpha(:,t)); +end +if scaled, assert(approxeq(sum(alpha(:,t)),1)), end +for t=2:T + %trans = transmat(:,:,act(t-1))'; + trans = transmat{act(t-1)}; + if maximize + m = max_mult(trans', alpha(:,t-1)); + %A = repmat(alpha(:,t-1), [1 Q]); + %m = max(trans .* A, [], 1); + else + m = trans' * alpha(:,t-1); + end + alpha(:,t) = m(:) .* obslik(:,t); + if scaled + %[alpha(:,t), scale(t)] = normaliseC(alpha(:,t)); + [alpha(:,t), scale(t)] = normalise(alpha(:,t)); + end + if compute_xi & fwd_only % useful for online EM + %xi(:,:,t-1) = normaliseC((alpha(:,t-1) * obslik(:,t)') .* trans); + xi(:,:,t-1) = normalise((alpha(:,t-1) * obslik(:,t)') .* trans); + end + if scaled, assert(approxeq(sum(alpha(:,t)),1)), end +end +if scaled + if any(scale==0) + loglik = -inf; + else + loglik = sum(log(scale)); + end +else + loglik = log(sum(alpha(:,T))); +end + +if fwd_only + gamma = alpha; + beta = []; + gamma2 = []; + return; +end + + +%%%%%%%%% Backwards %%%%%%%%%% + +beta = zeros(Q,T); +if compute_gamma2 + M = size(mixmat, 2); + gamma2 = zeros(Q,M,T); +else + gamma2 = []; +end + +beta(:,T) = ones(Q,1); +%gamma(:,T) = normaliseC(alpha(:,T) .* beta(:,T)); +gamma(:,T) = normalise(alpha(:,T) .* beta(:,T)); +t=T; +if compute_gamma2 + denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing + gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]); + %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); % wrong! +end +for t=T-1:-1:1 + b = beta(:,t+1) .* obslik(:,t+1); + %trans = transmat(:,:,act(t)); + trans = transmat{act(t)}; + if maximize + B = repmat(b(:)', Q, 1); + beta(:,t) = max(trans .* B, [], 2); + else + beta(:,t) = trans * b; + end + if scaled + %beta(:,t) = normaliseC(beta(:,t)); + beta(:,t) = normalise(beta(:,t)); + end + %gamma(:,t) = normaliseC(alpha(:,t) .* beta(:,t)); + gamma(:,t) = normalise(alpha(:,t) .* beta(:,t)); + if compute_xi + %xi(:,:,t) = normaliseC((trans .* (alpha(:,t) * b'))); + xi(:,:,t) = normalise((trans .* (alpha(:,t) * b'))); + %xi(:,:,t) = (trans .* (alpha(:,t) * b')); + end + if compute_gamma2 + denom = obslik(:,t) + (obslik(:,t)==0); % replace 0s with 1s before dividing + gamma2(:,:,t) = obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M]) ./ repmat(denom, [1 M]); + %gamma2(:,:,t) = normaliseC(obslik2(:,:,t) .* mixmat .* repmat(gamma(:,t), [1 M])); + end +end + + +% We now explain the equation for gamma2 +% Let zt=y(1:t-1,t+1:T) be all observations except y(t) +% gamma2(Q,M,t) = P(Qt,Mt|yt,zt) = P(yt|Qt,Mt,zt) P(Qt,Mt|zt) / P(yt|zt) +% = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|zt) / P(yt|zt) +% Now gamma(Q,t) = P(Qt|yt,zt) = P(yt|Qt) P(Qt|zt) / P(yt|zt) +% hence +% 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) +% = P(yt|Qt,Mt) P(Mt|Qt) P(Qt|yt,zt) / P(yt|Qt) +% \ No newline at end of file diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m new file mode 100644 index 00000000..3de17b40 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/hmm_inf_engine.m @@ -0,0 +1,71 @@ +function engine = hmm_inf_engine(bnet, varargin) +% HMM_INF_ENGINE Inference engine for DBNs which uses the forwards-backwards algorithm. +% engine = hmm_inf_engine(bnet, ...) +% +% The following optional arguments can be specified in the form of name/value pairs: +% [default value in brackets] +% +% maximize - 1 means max-product, 0 means sum-product [0] +% +% The DBN is converted to an HMM with a single meganode, but the observed nodes remain factored. +% This can be faster than jtree if the num. hidden nodes is low, because of lower constant factors. +% +% All hidden nodes must be discrete. +% All observed nodes are assumed to be leaves, i.e., they cannot be parents of anything. +% The parents of each observed leaf are assumed to be a subset of the hidden nodes within the same slice. +% The only exception is if bnet is an AR-HMM, where the parents are assumed to be self in the +% previous slice (continuous), plus all the discrete nodes in the current slice. + +ss = bnet.nnodes_per_slice; + +engine.maximize = 0; +% parse optional params +args = varargin; +nargs = length(args); +if nargs > 0 + for i=1:2:nargs + switch args{i}, + case 'maximize', engine.maximize = args{i+1}; + otherwise, + error(['invalid argument name ' args{i}]); + end + end +end + +% Stuff to do with speeding up marginal_family +[int, engine.persist, engine.transient] = compute_interface_nodes(bnet.intra, bnet.inter); +engine.persist_bitv = zeros(1, ss); +engine.persist_bitv(engine.persist) = 1; + + +ns = bnet.node_sizes(:); +ns(bnet.observed) = 1; +ns(bnet.observed+ss) = 1; +engine.eff_node_sizes = ns; + +for o=bnet.observed(:)' + %if bnet.equiv_class(o,1) ~= bnet.equiv_class(o,2) + % error(['observed node ' num2str(o) ' is not tied']) + %end + cs = children(bnet.dag, o); + if ~isempty(cs) + error(['observed node ' num2str(o) ' is not allowed children']) + end +end + +[engine.startprob, engine.transprob, engine.obsprob] = dbn_to_hmm(bnet); + +% This is where we will store the results between enter_evidence and marginal_nodes +engine.one_slice_marginal = []; +engine.two_slice_marginal = []; + +ss = length(bnet.intra); +engine.evidence = []; +engine.node_sizes = []; + +% avoid the need to do bnet_from_engine, which is slow +engine.slice_size = ss; +engine.parents = bnet.parents; + +engine = class(engine, 'hmm_inf_engine', inf_engine(bnet)); + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m new file mode 100644 index 00000000..56b9fb6c --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_family.m @@ -0,0 +1,35 @@ +function marginal = marginal_family(engine, i, t, add_ev) +% MARGINAL_FAMILY Compute the marginal on the specified family (hmm) +% marginal = marginal_family(engine, i, t, add_ev) + +if nargin < 3, t = 1; end +if nargin < 4, add_ev = 0; end + +ns = engine.eff_node_sizes(:); +ss = engine.slice_size; + +if t==1 | ~engine.persist_bitv(i) + bigT = engine.one_slice_marginal(:,t); + ps = engine.parents{i}; + dom = [ps i] + (t-1)*ss; + bigdom = 1:ss; + bigsz = ns(bigdom); + bigdom = bigdom + (t-1)*ss; +else % some parents are in previous slice + bigT = engine.two_slice_marginal(:,t-1); % t-1 and t + ps = engine.parents{i+ss}; + dom = [ps i+ss] + (t-2)*ss; + bigdom = 1:(2*ss); % domain of xi(:,:,t) + bigsz = ns(bigdom); + bigdom = bigdom + (t-2)*ss; +end +marginal.domain = dom; + +marginal.T = marg_table(bigT, bigdom, bigsz, dom, engine.maximize); +marginal.mu = []; +marginal.Sigma = []; + +if add_ev + marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m new file mode 100644 index 00000000..0a2bec4f --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/marginal_nodes.m @@ -0,0 +1,29 @@ +function marginal = marginal_nodes(engine, nodes, t, add_ev) +% MARGINAL_NODES Compute the marginal on the specified query nodes (hmm) +% marginal = marginal_nodes(engine, nodes, t, add_ev) +% +% 'nodes' must be a single node. +% t is the time slice. + +if nargin < 3, t = 1; end +if nargin < 4, add_ev = 0; end + +assert(length(nodes)==1) +ss = engine.slice_size; + +i = nodes(1); +bigT = engine.one_slice_marginal(:,t); +dom = i + (t-1)*ss; + +ns = engine.eff_node_sizes(:); +bigdom = 1:ss; +marginal.T = marg_table(bigT, bigdom + (t-1)*ss, ns(bigdom), dom, engine.maximize); + +marginal.domain = dom; +marginal.mu = []; +marginal.Sigma = []; + +if add_ev + marginal = add_ev_to_dmarginal(marginal, engine.evidence, engine.node_sizes); +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries new file mode 100644 index 00000000..a35185ea --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Entries @@ -0,0 +1,3 @@ +/mk_hmm_obs_lik_matrix.m/1.1.1.1/Sun May 4 21:42:26 2003// +/mk_hmm_obs_lik_vec.m/1.1.1.1/Thu Jan 23 18:50:10 2003// +D diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository new file mode 100644 index 00000000..20dfc6fd --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Repository @@ -0,0 +1 @@ +FullBNT/BNT/inference/dynamic/@hmm_inf_engine/private diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root new file mode 100644 index 00000000..f3bd14a6 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/CVS/Root @@ -0,0 +1 @@ +:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m new file mode 100644 index 00000000..441f3c0a --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_matrix.m @@ -0,0 +1,30 @@ +function obslik = mk_hmm_obs_lik_matrix(engine, evidence) + +T = size(evidence,2); +Q = length(engine.startprob); +obslik = ones(Q, T); +bnet = bnet_from_engine(engine); +% P(o1,o2| Q1,Q2) = P(o1|Q1,Q2) * P(o2|Q1,Q2) +onodes = bnet.observed; +for i=1:length(onodes) + data = cell2num(evidence(onodes(i),:)); + if bnet.auto_regressive(onodes(i)) + params = engine.obsprob{i}; + mu = params.big_mu; + Sigma = params.big_Sigma, + W = params.big_W; + mu0 = params.big_mu0; + Sigma0 = params.big_Sigma0; + %obslik_i = mk_arhmm_obs_lik(data, mu, Sigma, W, mu0, Sigma0 + obslik_i = clg_prob(data(:,1:T-1), data(:,2:T), mu, Sigma, W); + obslik_i = [mixgauss_prob(data(:,1), mu0, Sigma0) obslik_i]; + elseif myismember(onodes(i), bnet.dnodes) + %obslik_i = eval_pdf_cond_multinomial(data, engine.obsprob{i}.big_CPT); + obslik_i = multinomial_prob(data, engine.obsprob{i}.big_CPT); + else + %obslik_i = eval_pdf_cond_gauss(data, engine.obsprob{i}.big_mu, engine.obsprob{i}.big_Sigma); + obslik_i = mixgauss_prob(data, engine.obsprob{i}.big_mu, engine.obsprob{i}.big_Sigma); + end + obslik = obslik .* obslik_i; +end + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m new file mode 100644 index 00000000..16d30aec --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/private/mk_hmm_obs_lik_vec.m @@ -0,0 +1,52 @@ +function obslik = mk_hmm_obs_lik_vec(engine, evidence) + +% P(o1,o2| h) = P(o1|h) * P(o2|h) where h = Q1,Q2,... + +bnet = bnet_from_engine(engine); +ss = length(bnet.intra); +onodes = bnet.observed; +hnodes = mysetdiff(1:ss, onodes); +ns = bnet.node_sizes(:); +ns(onodes) = 1; + +Q = length(engine.startprob); +obslik = ones(Q, 1); + +for i=1:length(onodes) + o = onodes(i); + %data = cell2num(evidence(o,1)); + data = evidence{o,1}; + if myismember(o, bnet.dnodes) + obslik_i = eval_pdf_cond_multinomial(data, engine.obsprob{i}.CPT); + else + if bnet.auto_regressive(o) + error('can''t handle AR nodes') + end + %% calling mk_ghmm_obs_lik, which calls gaussian_prob, is slow, so we inline it + %% and use the pre-computed inverse matrix + %obslik_i = mk_ghmm_obs_lik(data, engine.obsprob{i}.mu, engine.obsprob{i}.Sigma); + x = data(:); + m = engine.obsprob{i}.mu; + Qi = size(m, 2); + obslik_i = size(Qi, 1); + invC = engine.obsprob{i}.inv_Sigma; + denom = engine.obsprob{i}.denom; + for j=1:Qi + numer = exp(-0.5 * (x-m(:,j))' * invC(:,:,j) * (x-m(:,j))); + obslik_i(j) = numer / denom(j); + end + end + % convert P(o|ps) into P(o|h) by multiplying onto a (h,o) potential of all 1s + ps = bnet.parents{o}; + dom = [ps o]; + obspot_i = dpot(dom, ns(dom), obslik_i); + dom = [hnodes o]; + obspot = dpot(dom, ns(dom)); + obspot = multiply_by_pot(obspot, obspot_i); + % compute p(oi|h) * p(oj|h) + S = struct(obspot); + obslik = obslik .* S.T(:); +end + + + diff --git a/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m new file mode 100644 index 00000000..e6cd1f79 --- /dev/null +++ b/sourcecodes/bnt-master/BNT/inference/dynamic/@hmm_inf_engine/update_engine.m @@ -0,0 +1,8 @@ +function engine = update_engine(engine, newCPDs) +% UPDATE_ENGINE Update the engine to take into account the new parameters (hmm) +% engine = update_engine(engine, newCPDs) + +%engine.inf_engine.bnet.CPD = newCPDs; +engine.inf_engine = update_engine(engine.inf_engine, newCPDs); +[engine.startprob, engine.transprob, engine.obsprob] = dbn_to_hmm(bnet_from_engine(engine)); + |