|
wolffd@0
|
1 function scg_unstable()
|
|
wolffd@0
|
2
|
|
wolffd@0
|
3 % the objective of this script is to test if the stable conditonal gaussian
|
|
wolffd@0
|
4 % inference can handle the numerical instability problem described on
|
|
wolffd@0
|
5 % page.151 of 'Probabilistic networks and expert system' by Cowell, Dawid, Lauritzen and
|
|
wolffd@0
|
6 % Spiegelhalter, 1999.
|
|
wolffd@0
|
7
|
|
wolffd@0
|
8 A = 1; Y = 2;
|
|
wolffd@0
|
9 n = 2;
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 ns = ones(1, n);
|
|
wolffd@0
|
12 dnodes = [A];
|
|
wolffd@0
|
13 cnodes = Y;
|
|
wolffd@0
|
14 ns = [2 1];
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 dag = zeros(n);
|
|
wolffd@0
|
17 dag(A, Y) = 1;
|
|
wolffd@0
|
18
|
|
wolffd@0
|
19 bnet = mk_bnet(dag, ns, dnodes);
|
|
wolffd@0
|
20
|
|
wolffd@0
|
21 bnet.CPD{A} = tabular_CPD(bnet, A, [0.5 0.5]');
|
|
wolffd@0
|
22 bnet.CPD{Y} = gaussian_CPD(bnet, Y, 'mean', [0 1], 'cov', [1e-5 1e-6]);
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 evidence = cell(1, n);
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 pot_type = 'cg';
|
|
wolffd@0
|
27 potYgivenA = convert_to_pot(bnet.CPD{Y}, pot_type, [A Y], evidence);
|
|
wolffd@0
|
28 potA = convert_to_pot(bnet.CPD{A}, pot_type, A, evidence);
|
|
wolffd@0
|
29 potYandA = multiply_by_pot(potYgivenA, potA);
|
|
wolffd@0
|
30 potA2 = marginalize_pot(potYandA, A);
|
|
wolffd@0
|
31
|
|
wolffd@0
|
32 thresh = 1; % 0dp
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 [g,h,K] = extract_can(potA);
|
|
wolffd@0
|
35 assert(approxeq(g(:)', [-0.693147 -0.693147], thresh))
|
|
wolffd@0
|
36
|
|
wolffd@0
|
37
|
|
wolffd@0
|
38 [g,h,K] = extract_can(potYgivenA);
|
|
wolffd@0
|
39 assert(approxeq(g(:)', [4.83752 -499994], thresh))
|
|
wolffd@0
|
40 assert(approxeq(h(:)', [0 1e6]))
|
|
wolffd@0
|
41 assert(approxeq(K(:)', [1e5 1e6]))
|
|
wolffd@0
|
42
|
|
wolffd@0
|
43 [g,h,K] = extract_can(potYandA);
|
|
wolffd@0
|
44 assert(approxeq(g(:)', [4.14437 -499995], thresh))
|
|
wolffd@0
|
45 assert(approxeq(h(:)', [0 1e6]))
|
|
wolffd@0
|
46 assert(approxeq(K(:)', [1e5 1e6]))
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48
|
|
wolffd@0
|
49 [g,h,K] = extract_can(potA2);
|
|
wolffd@0
|
50 %assert(approxeq(g(:)', [-0.69315 -1]))
|
|
wolffd@0
|
51 g
|
|
wolffd@0
|
52 assert(approxeq(g(:)', [-0.69315 -0.69315]))
|
|
wolffd@0
|
53
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55
|
|
wolffd@0
|
56 if 0
|
|
wolffd@0
|
57 pot_type = 'scg';
|
|
wolffd@0
|
58 spotYgivenA = convert_to_pot(bnet.CPD{Y}, pot_type, [A Y], evidence);
|
|
wolffd@0
|
59 spotA = convert_to_pot(bnet.CPD{A}, pot_type, A, evidence);
|
|
wolffd@0
|
60 spotYandA = direct_combine_pots(spotYgivenA, spotA);
|
|
wolffd@0
|
61 spotA2 = marginalize_pot(spotYandA, A);
|
|
wolffd@0
|
62
|
|
wolffd@0
|
63 spotA=struct(spotA);
|
|
wolffd@0
|
64 spotA2=struct(spotA2);
|
|
wolffd@0
|
65 for i=1:2
|
|
wolffd@0
|
66 assert(approxeq(spotA2.scgpotc{i}.p, spotA.scgpotc{i}.p))
|
|
wolffd@0
|
67 assert(approxeq(spotA2.scgpotc{i}.A, spotA.scgpotc{i}.A))
|
|
wolffd@0
|
68 assert(approxeq(spotA2.scgpotc{i}.B, spotA.scgpotc{i}.B))
|
|
wolffd@0
|
69 assert(approxeq(spotA2.scgpotc{i}.C, spotA.scgpotc{i}.C))
|
|
wolffd@0
|
70 end
|
|
wolffd@0
|
71
|
|
wolffd@0
|
72 end
|
|
wolffd@0
|
73
|
|
wolffd@0
|
74
|
|
wolffd@0
|
75 %%%%%%%%%%%
|
|
wolffd@0
|
76
|
|
wolffd@0
|
77 function [g,h,K] = extract_can(pot)
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79 pot = struct(pot);
|
|
wolffd@0
|
80 D = length(pot.can);
|
|
wolffd@0
|
81 g = zeros(1, D);
|
|
wolffd@0
|
82 h = zeros(1, D);
|
|
wolffd@0
|
83 K = zeros(1, D);
|
|
wolffd@0
|
84 for i=1:D
|
|
wolffd@0
|
85 S = struct(pot.can{i});
|
|
wolffd@0
|
86 g(i) = S.g;
|
|
wolffd@0
|
87 if length(S.h) > 0
|
|
wolffd@0
|
88 h(i) = S.h;
|
|
wolffd@0
|
89 K(i) = S.K;
|
|
wolffd@0
|
90 end
|
|
wolffd@0
|
91 end
|
