annotate toolboxes/FullBNT-1.0.7/bnt/CPDs/@gaussian_CPD/CPD_to_lambda_msg.m @ 0:cc4b1211e677 tip

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