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