|
wolffd@0
|
1 function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
|
|
wolffd@0
|
2 % CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)
|
|
wolffd@0
|
3 % lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p, evidence)
|
|
wolffd@0
|
4 % Pearl p183 eq 4.52
|
|
wolffd@0
|
5
|
|
wolffd@0
|
6 switch msg_type
|
|
wolffd@0
|
7 case 'd',
|
|
wolffd@0
|
8 error('gaussian_CPD can''t create discrete msgs')
|
|
wolffd@0
|
9 case 'g',
|
|
wolffd@0
|
10 cps = ps(CPD.cps);
|
|
wolffd@0
|
11 cpsizes = CPD.sizes(CPD.cps);
|
|
wolffd@0
|
12 self_size = CPD.sizes(end);
|
|
wolffd@0
|
13 i = find_equiv_posns(p, cps); % p is n's i'th cts parent
|
|
wolffd@0
|
14 psz = cpsizes(i);
|
|
wolffd@0
|
15 if all(msg{n}.lambda.precision == 0) % no info to send on
|
|
wolffd@0
|
16 lam_msg.precision = zeros(psz, psz);
|
|
wolffd@0
|
17 lam_msg.info_state = zeros(psz, 1);
|
|
wolffd@0
|
18 return;
|
|
wolffd@0
|
19 end
|
|
wolffd@0
|
20 [m, Q, W] = gaussian_CPD_params_given_dps(CPD, [ps n], evidence);
|
|
wolffd@0
|
21 Bmu = m;
|
|
wolffd@0
|
22 BSigma = Q;
|
|
wolffd@0
|
23 for k=1:length(cps) % only get pi msgs from cts parents
|
|
wolffd@0
|
24 pk = cps(k);
|
|
wolffd@0
|
25 if pk ~= p
|
|
wolffd@0
|
26 %bk = block(k, cpsizes);
|
|
wolffd@0
|
27 bk = CPD.cps_block_ndx{k};
|
|
wolffd@0
|
28 Bk = W(:, bk);
|
|
wolffd@0
|
29 m = msg{n}.pi_from_parent{k};
|
|
wolffd@0
|
30 BSigma = BSigma + Bk * m.Sigma * Bk';
|
|
wolffd@0
|
31 Bmu = Bmu + Bk * m.mu;
|
|
wolffd@0
|
32 end
|
|
wolffd@0
|
33 end
|
|
wolffd@0
|
34 % BSigma = Q + sum_{k \neq i} B_k Sigma_k B_k'
|
|
wolffd@0
|
35 %bi = block(i, cpsizes);
|
|
wolffd@0
|
36 bi = CPD.cps_block_ndx{i};
|
|
wolffd@0
|
37 Bi = W(:,bi);
|
|
wolffd@0
|
38 P = msg{n}.lambda.precision;
|
|
wolffd@0
|
39 if (rcond(P) > 1e-3) | isinf(P)
|
|
wolffd@0
|
40 if isinf(P) % Y is observed
|
|
wolffd@0
|
41 Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
|
|
wolffd@0
|
42 mu_lambda = msg{n}.lambda.mu; % observed_value;
|
|
wolffd@0
|
43 else
|
|
wolffd@0
|
44 Sigma_lambda = inv(P);
|
|
wolffd@0
|
45 mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
|
|
wolffd@0
|
46 end
|
|
wolffd@0
|
47 C = inv(Sigma_lambda + BSigma);
|
|
wolffd@0
|
48 lam_msg.precision = Bi' * C * Bi;
|
|
wolffd@0
|
49 lam_msg.info_state = Bi' * C * (mu_lambda - Bmu);
|
|
wolffd@0
|
50 else
|
|
wolffd@0
|
51 % method that uses matrix inversion lemma to avoid inverting P
|
|
wolffd@0
|
52 A = inv(P + inv(BSigma));
|
|
wolffd@0
|
53 C = P - P*A*P;
|
|
wolffd@0
|
54 lam_msg.precision = Bi' * C * Bi;
|
|
wolffd@0
|
55 D = eye(self_size) - P*A;
|
|
wolffd@0
|
56 z = msg{n}.lambda.info_state;
|
|
wolffd@0
|
57 lam_msg.info_state = Bi' * (D*z - D*P*Bmu);
|
|
wolffd@0
|
58 end
|
|
wolffd@0
|
59 end
|
