|
wolffd@0
|
1 % Lawn sprinker example from Russell and Norvig p454
|
|
wolffd@0
|
2 % See www.cs.berkeley.edu/~murphyk/Bayes/usage.html for details.
|
|
wolffd@0
|
3
|
|
wolffd@0
|
4 N = 4;
|
|
wolffd@0
|
5 dag = zeros(N,N);
|
|
wolffd@0
|
6 C = 1; S = 2; R = 3; W = 4;
|
|
wolffd@0
|
7 dag(C,[R S]) = 1;
|
|
wolffd@0
|
8 dag(R,W) = 1;
|
|
wolffd@0
|
9 dag(S,W)=1;
|
|
wolffd@0
|
10
|
|
wolffd@0
|
11 false = 1; true = 2;
|
|
wolffd@0
|
12 ns = 2*ones(1,N); % binary nodes
|
|
wolffd@0
|
13
|
|
wolffd@0
|
14 bnet = mk_bnet(dag, ns);
|
|
wolffd@0
|
15 bnet.CPD{C} = tabular_CPD(bnet, C, [0.5 0.5]);
|
|
wolffd@0
|
16 bnet.CPD{R} = tabular_CPD(bnet, R, [0.8 0.2 0.2 0.8]);
|
|
wolffd@0
|
17 bnet.CPD{S} = tabular_CPD(bnet, S, [0.5 0.9 0.5 0.1]);
|
|
wolffd@0
|
18 bnet.CPD{W} = tabular_CPD(bnet, W, [1 0.1 0.1 0.01 0 0.9 0.9 0.99]);
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 CPT = cell(1,N);
|
|
wolffd@0
|
21 for i=1:N
|
|
wolffd@0
|
22 s=struct(bnet.CPD{i}); % violate object privacy
|
|
wolffd@0
|
23 CPT{i}=s.CPT;
|
|
wolffd@0
|
24 end
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 % Generate training data
|
|
wolffd@0
|
27 nsamples = 50;
|
|
wolffd@0
|
28 samples = cell(N, nsamples);
|
|
wolffd@0
|
29 for i=1:nsamples
|
|
wolffd@0
|
30 samples(:,i) = sample_bnet(bnet);
|
|
wolffd@0
|
31 end
|
|
wolffd@0
|
32 data = cell2num(samples);
|
|
wolffd@0
|
33
|
|
wolffd@0
|
34 % Make a tabula rasa
|
|
wolffd@0
|
35 bnet2 = mk_bnet(dag, ns);
|
|
wolffd@0
|
36 seed = 0;
|
|
wolffd@0
|
37 rand('state', seed);
|
|
wolffd@0
|
38 bnet2.CPD{C} = tabular_CPD(bnet2, C, 'clamped', 1, 'CPT', [0.5 0.5], ...
|
|
wolffd@0
|
39 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
|
|
wolffd@0
|
40 bnet2.CPD{R} = tabular_CPD(bnet2, R, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
|
|
wolffd@0
|
41 bnet2.CPD{S} = tabular_CPD(bnet2, S, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
|
|
wolffd@0
|
42 bnet2.CPD{W} = tabular_CPD(bnet2, W, 'prior_type', 'dirichlet', 'dirichlet_weight', 0);
|
|
wolffd@0
|
43
|
|
wolffd@0
|
44
|
|
wolffd@0
|
45 % Find MLEs from fully observed data
|
|
wolffd@0
|
46 bnet4 = learn_params(bnet2, samples);
|
|
wolffd@0
|
47
|
|
wolffd@0
|
48 % Bayesian updating with 0 prior is equivalent to ML estimation
|
|
wolffd@0
|
49 bnet5 = bayes_update_params(bnet2, samples);
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 CPT4 = cell(1,N);
|
|
wolffd@0
|
52 for i=1:N
|
|
wolffd@0
|
53 s=struct(bnet4.CPD{i}); % violate object privacy
|
|
wolffd@0
|
54 CPT4{i}=s.CPT;
|
|
wolffd@0
|
55 end
|
|
wolffd@0
|
56
|
|
wolffd@0
|
57 CPT5 = cell(1,N);
|
|
wolffd@0
|
58 for i=1:N
|
|
wolffd@0
|
59 s=struct(bnet5.CPD{i}); % violate object privacy
|
|
wolffd@0
|
60 CPT5{i}=s.CPT;
|
|
wolffd@0
|
61 assert(approxeq(CPT5{i}, CPT4{i}))
|
|
wolffd@0
|
62 end
|
|
wolffd@0
|
63
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 if 1
|
|
wolffd@0
|
66 % Find MLEs from partially observed data
|
|
wolffd@0
|
67
|
|
wolffd@0
|
68 % hide 50% of the nodes
|
|
wolffd@0
|
69 samplesH = samples;
|
|
wolffd@0
|
70 hide = rand(N, nsamples) > 0.5;
|
|
wolffd@0
|
71 [I,J]=find(hide);
|
|
wolffd@0
|
72 for k=1:length(I)
|
|
wolffd@0
|
73 samplesH{I(k), J(k)} = [];
|
|
wolffd@0
|
74 end
|
|
wolffd@0
|
75
|
|
wolffd@0
|
76 engine = jtree_inf_engine(bnet2);
|
|
wolffd@0
|
77 max_iter = 5;
|
|
wolffd@0
|
78 [bnet6, LL] = learn_params_em(engine, samplesH, max_iter);
|
|
wolffd@0
|
79
|
|
wolffd@0
|
80 CPT6 = cell(1,N);
|
|
wolffd@0
|
81 for i=1:N
|
|
wolffd@0
|
82 s=struct(bnet6.CPD{i}); % violate object privacy
|
|
wolffd@0
|
83 CPT6{i}=s.CPT;
|
|
wolffd@0
|
84 end
|
|
wolffd@0
|
85
|
|
wolffd@0
|
86 end
|
