comparison toolboxes/FullBNT-1.0.7/bnt/examples/static/SCG/scg_unstable.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function scg_unstable()
+
+% the objective of this script is to test if the stable conditonal gaussian
+% inference can handle the numerical instability problem described on
+% page.151 of 'Probabilistic networks and expert system' by Cowell, Dawid, Lauritzen and
+% Spiegelhalter, 1999.
+
+A = 1; Y = 2;
+n = 2;
+
+ns = ones(1, n);
+dnodes = [A];
+cnodes = Y;
+ns = [2 1];
+
+dag = zeros(n);
+dag(A, Y) = 1;
+
+bnet = mk_bnet(dag, ns, dnodes);
+
+bnet.CPD{A} = tabular_CPD(bnet, A, [0.5 0.5]'); 
+bnet.CPD{Y} = gaussian_CPD(bnet, Y, 'mean', [0 1], 'cov', [1e-5 1e-6]);
+
+evidence = cell(1, n);
+
+pot_type = 'cg';
+potYgivenA = convert_to_pot(bnet.CPD{Y}, pot_type, [A Y], evidence);
+potA = convert_to_pot(bnet.CPD{A}, pot_type, A, evidence);
+potYandA = multiply_by_pot(potYgivenA, potA);
+potA2 = marginalize_pot(potYandA, A);
+
+thresh = 1; % 0dp
+
+[g,h,K] = extract_can(potA);
+assert(approxeq(g(:)', [-0.693147 -0.693147], thresh))
+
+
+[g,h,K] = extract_can(potYgivenA);
+assert(approxeq(g(:)', [4.83752 -499994], thresh))
+assert(approxeq(h(:)', [0 1e6]))
+assert(approxeq(K(:)', [1e5 1e6]))
+
+[g,h,K] = extract_can(potYandA);
+assert(approxeq(g(:)', [4.14437 -499995], thresh))
+assert(approxeq(h(:)', [0 1e6]))
+assert(approxeq(K(:)', [1e5 1e6]))
+
+
+[g,h,K] = extract_can(potA2);
+%assert(approxeq(g(:)', [-0.69315 -1]))
+g
+assert(approxeq(g(:)', [-0.69315 -0.69315]))
+
+
+
+if 0
+pot_type = 'scg';
+spotYgivenA = convert_to_pot(bnet.CPD{Y}, pot_type, [A Y], evidence);
+spotA = convert_to_pot(bnet.CPD{A}, pot_type, A, evidence);
+spotYandA = direct_combine_pots(spotYgivenA, spotA); 
+spotA2 = marginalize_pot(spotYandA, A);
+
+spotA=struct(spotA);
+spotA2=struct(spotA2);
+for i=1:2
+  assert(approxeq(spotA2.scgpotc{i}.p, spotA.scgpotc{i}.p))
+  assert(approxeq(spotA2.scgpotc{i}.A, spotA.scgpotc{i}.A))
+  assert(approxeq(spotA2.scgpotc{i}.B, spotA.scgpotc{i}.B))
+  assert(approxeq(spotA2.scgpotc{i}.C, spotA.scgpotc{i}.C))
+end
+
+end
+
+
+%%%%%%%%%%%
+
+function [g,h,K] = extract_can(pot)
+
+pot = struct(pot);
+D = length(pot.can);
+g = zeros(1, D);
+h = zeros(1, D);
+K = zeros(1, D);
+for i=1:D
+  S = struct(pot.can{i});
+  g(i) = S.g;
+  if length(S.h) > 0
+    h(i) = S.h;
+    K(i) = S.K;
+  end
+end