Mercurial > hg > camir-aes2014
view toolboxes/FullBNT-1.0.7/bnt/inference/static/@var_elim_inf_engine/find_mpe.m @ 0:e9a9cd732c1e tip
first hg version after svn
author | wolffd |
---|---|
date | Tue, 10 Feb 2015 15:05:51 +0000 |
parents | |
children |
line wrap: on
line source
function mpe = find_mpe(engine, new_evidence, max_over) % FIND_MPE Find the most probable explanation of the data (assignment to the hidden nodes) % function mpe = find_mpe(engine, evidence, order) % % PURPOSE: % CALC_MPE Computes the most probable explanation to the network nodes % given the evidence. % % [mpe, ll] = calc_mpe(engine, new_evidence, max_over) % % INPUT: % bnet - the bayesian network % new_evidence - optional, if specified - evidence to be incorporated [cell(1,n)] % max_over - optional, if specified determines the variable elimination order [1:n] % % OUTPUT: % mpe - the MPE assignmet for the net variables (or [] if no satisfying assignment) % ll - log assignment probability. % % Notes: % 1. Adapted from '@var_elim_inf_engine\marginal_nodes' for MPE by Ron Zohar, 8/7/01 % 2. Only discrete potentials are supported at this time. % 3. Complexity: O(nw*) where n is the number of nodes and w* is the induced tree width. % 4. Implementation based on: % - R. Dechter, "Bucket Elimination: A Unifying Framework for Probabilistic Inference", % UA1 96, pp. 211-219. bnet = bnet_from_engine(engine); ns = bnet.node_sizes; n = length(bnet.dag); evidence = cell(1,n); if (nargin<2) new_evidence = evidence; end onodes = find(~isemptycell(new_evidence)); % observed nodes hnodes = find(isemptycell(new_evidence)); % hidden nodes pot_type = determine_pot_type(bnet, onodes); if pot_type ~= 'd' error('only disrete potentials supported at this time') end for i=1:n fam = family(bnet.dag, i); CPT{i} = convert_to_pot(bnet.CPD{bnet.equiv_class(i)}, pot_type, fam(:), evidence); end % handle observed nodes: set impossible cases' probability to zero % rather than prun matrix (this makes backtracking easier) for ii=onodes lIdx = 1:ns(ii); lIdx = setdiff(lIdx, new_evidence{ii}); sCPT=struct(CPT{ii}); % violate object privacy sargs = ''; for jj=1:(length(sCPT.domain)-1) sargs = [sargs, ':,']; end for jj=lIdx eval(['sCPT.T(', sargs, num2str(jj), ')=0;']); end CPT{ii}=dpot(sCPT.domain, sCPT.sizes, sCPT.T); end B = cell(1,n); for b=1:n B{b} = mk_initial_pot(pot_type, [], [], [], []); end if (nargin<3) max_over = (1:n); end order = max_over; % no attempt to optimize this % Initialize the buckets with the CPDs assigned to them for i=1:n b = bucket_num(domain_pot(CPT{i}), order); B{b} = multiply_pots(B{b}, CPT{i}); end % Do backward phase max_over = max_over(length(max_over):-1:1); % reverse maximize = 1; for i=max_over(1:end-1) % max-ing over variable i which occurs in bucket j j = bucket_num(i, order); rest = mysetdiff(domain_pot(B{j}), i); %temp = marginalize_pot_max(B{j}, rest); temp = marginalize_pot(B{j}, rest, maximize); b = bucket_num(domain_pot(temp), order); % fprintf('maxing over bucket %d (var %d), putting result into bucket %d\n', j, i, b); sB=struct(B{b}); % violate object privacy if ~isempty(sB.domain) B{b} = multiply_pots(B{b}, temp); else B{b} = temp; end end result = B{1}; marginal = pot_to_marginal(result); [prob, mpe] = max(marginal.T); % handle impossible cases if ~(prob>0) mpe = []; ll = -inf; %warning('evidence has zero probability') return end ll = log(prob); % Do forward phase for ii=2:n marginal = pot_to_marginal(B{ii}); mpeidx = []; for jj=order(1:length(mpe)) %assert(ismember(jj, marginal.domain)) %%% bug temp = find_equiv_posns(jj, marginal.domain); mpeidx = [mpeidx, temp] ; if isempty(temp) mpeidx = [mpeidx, Inf] ; end end [mpeidxsorted sortedtompe] = sort(mpeidx) ; % maximize the matrix obtained from assigning values from previous buckets. % this is done by building a string and using eval. kk=1; sargs = '('; for jj=1:length(marginal.domain) if (jj~=1) sargs = [sargs, ',']; end if (mpeidxsorted(kk)==jj) sargs = [sargs, num2str(mpe(sortedtompe(kk)))]; if (kk<length(mpe)) kk = kk+1 ; end else sargs = [sargs, ':']; end end sargs = [sargs, ')'] ; eval(['[val, loc] = max(marginal.T', sargs, ');']) mpe = [mpe loc]; end [I,J] = sort(order); mpe = mpe(J); mpe = num2cell(mpe); %%%%%%%%% function b = bucket_num(domain, order) b = max(find_equiv_posns(domain, order));