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