wolffd@0: function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
wolffd@0: % UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node
wolffd@0: % function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
wolffd@0: 
wolffd@0: %if nargin < 6
wolffd@0: %  hidden_bitv = zeros(1, max(fmarginal.domain));
wolffd@0: %  hidden_bitv(find(isempty(evidence)))=1;
wolffd@0: %end
wolffd@0: 
wolffd@0: dom = fmarginal.domain;
wolffd@0: self = dom(end);
wolffd@0: ps = dom(1:end-1);
wolffd@0: hidden_self = hidden_bitv(self);
wolffd@0: cps = myintersect(ps, cnodes);
wolffd@0: dps = mysetdiff(ps, cps);
wolffd@0: hidden_cps = all(hidden_bitv(cps));
wolffd@0: hidden_dps = all(hidden_bitv(dps));
wolffd@0: 
wolffd@0: CPD.nsamples = CPD.nsamples + 1;            
wolffd@0: [ss cpsz dpsz] = size(CPD.weights); % ss = self size
wolffd@0: 
wolffd@0: % Let X be the cts parent (if any), Y be the cts child (self).
wolffd@0: 
wolffd@0: if ~hidden_self & (isempty(cps) | ~hidden_cps) & hidden_dps % all cts nodes are observed, all discrete nodes are hidden
wolffd@0:   % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0
wolffd@0:   % Since discrete parents are hidden, we do not need to add evidence to w.
wolffd@0:   w = fmarginal.T(:);
wolffd@0:   CPD.Wsum = CPD.Wsum + w;
wolffd@0:   y = evidence{self};
wolffd@0:   Cyy = y*y';
wolffd@0:   if ~CPD.useC
wolffd@0:      W = repmat(w(:)',ss,1); % W(y,i) = w(i)
wolffd@0:      W2 = repmat(reshape(W, [ss 1 dpsz]), [1 ss 1]); % W2(x,y,i) = w(i)
wolffd@0:      CPD.WYsum = CPD.WYsum +  W .* repmat(y(:), 1, dpsz);
wolffd@0:      CPD.WYYsum = CPD.WYYsum + W2  .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]);
wolffd@0:   else
wolffd@0:      W = w(:)';
wolffd@0:      W2 = reshape(W, [1 1 dpsz]);
wolffd@0:      CPD.WYsum = CPD.WYsum +  rep_mult(W, y(:), size(CPD.WYsum)); 
wolffd@0:      CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum));
wolffd@0:   end
wolffd@0:   if cpsz > 0 % X exists
wolffd@0:     x = cat(1, evidence{cps}); x = x(:);
wolffd@0:     Cxx = x*x';
wolffd@0:     Cxy = x*y';
wolffd@0:     if ~CPD.useC
wolffd@0:        CPD.WXsum = CPD.WXsum + W .* repmat(x(:), 1, dpsz);
wolffd@0:        CPD.WXXsum = CPD.WXXsum + W2 .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]);
wolffd@0:        CPD.WXYsum = CPD.WXYsum + W2 .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]);
wolffd@0:     else
wolffd@0:        CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum));
wolffd@0:        CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum));
wolffd@0:        CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum));
wolffd@0:     end
wolffd@0:   end
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: % general (non-vectorized) case
wolffd@0: fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow!
wolffd@0: 
wolffd@0: if dpsz == 1 % no discrete parents
wolffd@0:   w = 1;
wolffd@0: else
wolffd@0:   w = fullm.T(:);
wolffd@0: end
wolffd@0: 
wolffd@0: CPD.Wsum = CPD.Wsum + w;
wolffd@0: xi = 1:cpsz;
wolffd@0: yi = (cpsz+1):(cpsz+ss);
wolffd@0: for i=1:dpsz
wolffd@0:   muY = fullm.mu(yi, i);
wolffd@0:   SYY = fullm.Sigma(yi, yi, i);
wolffd@0:   CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
wolffd@0:   CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
wolffd@0:   if cpsz > 0
wolffd@0:     muX = fullm.mu(xi, i);
wolffd@0:     SXX = fullm.Sigma(xi, xi, i);
wolffd@0:     SXY = fullm.Sigma(xi, yi, i);
wolffd@0:     CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
wolffd@0:     CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
wolffd@0:     CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
wolffd@0:   end
wolffd@0: end                
wolffd@0: