fg1.m

Bayesian网络工具箱.

% make an unrolled HMM, convert to factor graph, and check that 
% loopy propagation on the fgraph gives the exact answers.

seed = 1;
rand('state', seed);
randn('state', seed);

T = 3;
H = 1:T;
O = H+T;
N = 2*T;
G = zeros(N);
for t=1:T
  G(H(t),O(t)) = 1;
end
for t=1:T-1
  G(H(t),H(t+1)) = 1;
end
ns = ones(1,N);
ns(H) = 2; % binary hidden
ns(O) = 3; % ternary observed
bnet = mk_bnet(G, ns, 'discrete', 1:N);
% no param tieing
for i=1:N
  bnet.CPD{i} = tabular_CPD(bnet, i);
end

fgraph = bnet_to_fgraph(bnet);
big_bnet = fgraph_to_bnet(fgraph);
% converting factor graph back does not recover the structure of the original bnet

max_iter = 2*T;

engine = {};
engine{1} = jtree_inf_engine(bnet);
engine{2} = belprop_inf_engine(bnet, 'max_iter', max_iter);
engine{3} = belprop_fg_inf_engine(fgraph, 'max_iter', max_iter);
engine{4} = jtree_inf_engine(big_bnet);
nengines = length(engine);

big_engine = 4;
fgraph_engine = 3;


data = sample_bnet(bnet, 1);
evidence = cell(1,N);
evidence(O) = data(O);

bigN = length(big_bnet.dag);
big_evidence = cell(1, bigN);
big_evidence(O) = data(O);
big_evidence(N+1:end) = {1}; % factors are observed to be 1

ll = zeros(1, nengines);
for i=1:nengines
  if i==big_engine
    tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, big_evidence); toc
  else
    tic; [engine{i}, ll(i)] = enter_evidence(engine{i}, evidence); toc
  end
end

% compare all engines to engine{1}

% the log likelihood values may be bogus...
for i=2:nengines
  %assert(approxeq(ll(1), ll(i)));
end

m = cell(nengines, T);
for i=1:T
  for e=1:nengines
    m{e,i} = marginal_nodes(engine{e}, H(i));
  end
  for e=2:nengines
    assert(approxeq(m{e,i}.T, m{1,i}.T));
  end
end


ll = zeros(1, nengines);
for e=1:nengines
  if e==big_engine
    [mpe{e}, ll(e)] = calc_mpe(engine{e}, big_evidence);
    mpe{e} = mpe{e}(1:N); % chop off dummy nodes
  elseif e==fgraph_engine
    mpe{e} = calc_mpe(engine{e}, evidence);
  else
    [mpe{e}, ll(e)] = calc_mpe(engine{e}, evidence);
  end
end
ll
celldisp(mpe)

% fgraph can't compute loglikelihood for software reasons
% jtree on the big_bnet gives the wrong ll
for e=2:nengines
  %assert(approxeq(ll(1), ll(e)));
  assert(approxeq(mpe{1}, mpe{e}))
end