BNW using Octave instead of Matlab.

This version of BNW should perform the same as the original version. The only difference is that it uses Octave instead of Matlab when running BayesNet Toolbox during parameter learning.

I am calling this BNW_1.02. It can be accessed at:
compbio.uthsc.edu/BNW_1.02

9 files changed, 666 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m
new file mode 100644
index 00000000..5a6d398a
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CPD_to_lambda_msg.m
@@ -0,0 +1,64 @@
+function lam_msg = CPD_to_lambda_msg(CPD, msg_type, n, ps, msg, p)
+% CPD_TO_LAMBDA_MSG Compute lambda message (gaussian)
+% lam_msg = compute_lambda_msg(CPD, msg_type, n, ps, msg, p)
+% Pearl p183 eq 4.52
+
+switch msg_type
+ case 'd',
+  error('gaussian_CPD can''t create discrete msgs')
+ case 'g',
+  self_size = CPD.sizes(end);
+  if all(msg{n}.lambda.precision == 0) % no info to send on
+    lam_msg.precision = zeros(self_size);
+    lam_msg.info_state = zeros(self_size, 1);
+    return;
+  end
+  cpsizes = CPD.sizes(CPD.cps);
+  dpval = 1;
+  Q = CPD.cov(:,:,dpval);
+  Sigmai = Q;
+  wmu = zeros(self_size, 1);
+  for k=1:length(ps)
+    pk = ps(k);
+    if pk ~= p
+      bk = block(k, cpsizes);
+      Bk = CPD.weights(:, bk, dpval);
+      m = msg{n}.pi_from_parent{k};
+      Sigmai = Sigmai + Bk * m.Sigma * Bk';
+      wmu = wmu + Bk * m.mu; % m.mu = u(k)
+    end
+  end
+  % Sigmai = Q + sum_{k \neq i} B_k Sigma_k B_k'
+  i = find_equiv_posns(p, ps);
+  bi = block(i, cpsizes);
+  Bi = CPD.weights(:,bi, dpval);
+  
+  if 0
+  P = msg{n}.lambda.precision;
+  if isinf(P) % inv(P)=Sigma_lambda=0
+    precision_temp = inv(Sigmai);
+    lam_msg.precision = Bi' * precision_temp * Bi;
+    lam_msg.info_state = precision_temp * (msg{n}.lambda.mu - wmu);
+  else
+    A = inv(P + inv(Sigmai));
+    precision_temp = P + P*A*P;
+    lam_msg.precision = Bi' * precision_temp * Bi;
+    self_size = length(P);
+    C = eye(self_size) + P*A;
+    z = msg{n}.lambda.info_state;
+    lam_msg.info_state = C*z - C*P*wmu;
+  end
+  end
+  
+  if isinf(msg{n}.lambda.precision)
+    Sigma_lambda = zeros(self_size, self_size); % infinite precision => 0 variance
+    mu_lambda = msg{n}.lambda.mu; % observed_value;
+  else
+    Sigma_lambda = inv(msg{n}.lambda.precision);
+    mu_lambda = Sigma_lambda * msg{n}.lambda.info_state;
+  end
+  precision_temp = inv(Sigma_lambda + Sigmai);
+  lam_msg.precision = Bi' * precision_temp * Bi;
+  lam_msg.info_state = Bi' * precision_temp * (mu_lambda - wmu);
+end
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries
new file mode 100644
index 00000000..ea2f5a4c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Entries
@@ -0,0 +1,7 @@
+/CPD_to_lambda_msg.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/gaussian_CPD.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/log_prob_node.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/maximize_params.m/1.1.1.1/Thu Jan 30 22:38:16 2003//
+/update_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//
+/update_tied_ess.m/1.1.1.1/Wed May 29 15:59:52 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository
new file mode 100644
index 00000000..c89b5b86
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/CPDs/@gaussian_CPD/Old
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m
new file mode 100644
index 00000000..6f7138fc
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/gaussian_CPD.m
@@ -0,0 +1,184 @@
+function CPD = gaussian_CPD(varargin)
+% GAUSSIAN_CPD Make a conditional linear Gaussian distrib.
+%
+% To define this CPD precisely, call the continuous (cts) parents (if any) X,
+% the discrete parents (if any) Q, and this node Y. Then the distribution on Y is:
+% - no parents: Y ~ N(mu, Sigma)
+% - cts parents : Y|X=x ~ N(mu + W x, Sigma)
+% - discrete parents: Y|Q=i ~ N(mu(i), Sigma(i))
+% - cts and discrete parents: Y|X=x,Q=i ~ N(mu(i) + W(i) x, Sigma(i))
+%
+% CPD = gaussian_CPD(bnet, node, ...) will create a CPD with random parameters,
+% where node is the number of a node in this equivalence class.
+%
+% The list below gives optional arguments [default value in brackets].
+% (Let ns(i) be the size of node i, X = ns(X), Y = ns(Y) and Q = prod(ns(Q)).)
+%
+% mean       - mu(:,i) is the mean given Q=i [ randn(Y,Q) ]
+% cov        - Sigma(:,:,i) is the covariance given Q=i [ repmat(eye(Y,Y), [1 1 Q]) ]
+% weights    - W(:,:,i) is the regression matrix given Q=i [ randn(Y,X,Q) ]
+% cov_type   - if 'diag', Sigma(:,:,i) is diagonal [ 'full' ]
+% tied_cov   - if 1, we constrain Sigma(:,:,i) to be the same for all i [0]
+% clamp_mean - if 1, we do not adjust mu(:,i) during learning [0]
+% clamp_cov  - if 1, we do not adjust Sigma(:,:,i) during learning [0]
+% clamp_weights - if 1, we do not adjust W(:,:,i) during learning [0]
+% cov_prior_weight - weight given to I prior for estimating Sigma [0.01]
+%
+% e.g., CPD = gaussian_CPD(bnet, i, 'mean', [0; 0], 'clamp_mean', 'yes')
+%
+% For backwards compatibility with BNT2, you can also specify the parameters in the following order
+%   CPD = gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weight)
+%
+% Sometimes it is useful to create an "isolated" CPD, without needing to pass in a bnet.
+% In this case, you must specify the discrete and cts parents (dps, cps) and the family sizes, followed
+% by the optional arguments above:
+%   CPD = gaussian_CPD('self', i, 'dps', dps, 'cps', cps, 'sz', fam_size, ...)
+
+
+if nargin==0
+  % This occurs if we are trying to load an object from a file.
+  CPD = init_fields;
+  clamp = 0;
+  CPD = class(CPD, 'gaussian_CPD', generic_CPD(clamp));
+  return;
+elseif isa(varargin{1}, 'gaussian_CPD')
+  % This might occur if we are copying an object.
+  CPD = varargin{1};
+  return;
+end
+CPD = init_fields;
+ 
+CPD = class(CPD, 'gaussian_CPD', generic_CPD(0));
+
+
+% parse mandatory arguments
+if ~isstr(varargin{1}) % pass in bnet
+  bnet = varargin{1};
+  self = varargin{2};
+  args = varargin(3:end);
+  ns = bnet.node_sizes;
+  ps = parents(bnet.dag, self);
+  dps = myintersect(ps, bnet.dnodes);
+  cps = myintersect(ps, bnet.cnodes);
+  fam_sz = ns([ps self]);
+else
+  disp('parsing new style')
+  for i=1:2:length(varargin)
+    switch varargin{i},
+     case 'self', self = varargin{i+1}; 
+     case 'dps',  dps = varargin{i+1};
+     case 'cps',  cps = varargin{i+1};
+     case 'sz',   fam_sz = varargin{i+1};
+    end
+  end
+  ps = myunion(dps, cps);
+  args = varargin;
+end
+
+CPD.self = self;
+CPD.sizes = fam_sz;
+
+% Figure out which (if any) of the parents are discrete, and which cts, and how big they are
+% dps = discrete parents, cps = cts parents
+CPD.cps = find_equiv_posns(cps, ps); % cts parent index
+CPD.dps = find_equiv_posns(dps, ps);
+ss = fam_sz(end);
+psz = fam_sz(1:end-1);
+dpsz = prod(psz(CPD.dps));
+cpsz = sum(psz(CPD.cps));
+
+% set default params
+CPD.mean = randn(ss, dpsz);
+CPD.cov = 100*repmat(eye(ss), [1 1 dpsz]);    
+CPD.weights = randn(ss, cpsz, dpsz);
+CPD.cov_type = 'full';
+CPD.tied_cov = 0;
+CPD.clamped_mean = 0;
+CPD.clamped_cov = 0;
+CPD.clamped_weights = 0;
+CPD.cov_prior_weight = 0.01;
+
+nargs = length(args);
+if nargs > 0
+  if ~isstr(args{1})
+    % gaussian_CPD(bnet, self, mu, Sigma, W, cov_type, tied_cov, clamp_mean, clamp_cov, clamp_weights)
+    if nargs >= 1 & ~isempty(args{1}), CPD.mean = args{1}; end
+    if nargs >= 2 & ~isempty(args{2}), CPD.cov = args{2}; end
+    if nargs >= 3 & ~isempty(args{3}), CPD.weights = args{3}; end
+    if nargs >= 4 & ~isempty(args{4}), CPD.cov_type = args{4}; end
+    if nargs >= 5 & ~isempty(args{5}) & strcmp(args{5}, 'tied'), CPD.tied_cov = 1; end
+    if nargs >= 6 & ~isempty(args{6}), CPD.clamped_mean = 1; end
+    if nargs >= 7 & ~isempty(args{7}), CPD.clamped_cov = 1; end
+    if nargs >= 8 & ~isempty(args{8}), CPD.clamped_weights = 1; end
+  else
+    CPD = set_fields(CPD, args{:});
+  end
+end
+
+% Make sure the matrices have 1 dimension per discrete parent.
+% Bug fix due to Xuejing Sun 3/6/01
+CPD.mean = myreshape(CPD.mean, [ss ns(dps)]);
+CPD.cov = myreshape(CPD.cov, [ss ss ns(dps)]);
+CPD.weights = myreshape(CPD.weights, [ss cpsz ns(dps)]);
+  
+CPD.init_cov = CPD.cov;  % we reset to this if things go wrong during learning
+
+% expected sufficient statistics 
+CPD.Wsum = zeros(dpsz,1);
+CPD.WYsum = zeros(ss, dpsz);
+CPD.WXsum = zeros(cpsz, dpsz);
+CPD.WYYsum = zeros(ss, ss, dpsz);
+CPD.WXXsum = zeros(cpsz, cpsz, dpsz);
+CPD.WXYsum = zeros(cpsz, ss, dpsz);
+
+% For BIC
+CPD.nsamples = 0;
+switch CPD.cov_type
+  case 'full',
+    ncov_params = ss*(ss-1)/2; % since symmetric (and positive definite)
+  case 'diag',
+    ncov_params = ss;
+  otherwise
+    error(['unrecognized cov_type ' cov_type]);
+end
+% params = weights + mean + cov
+if CPD.tied_cov
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + ncov_params;
+else
+  CPD.nparams = ss*cpsz*dpsz + ss*dpsz + dpsz*ncov_params;
+end
+
+
+
+clamped = CPD.clamped_mean & CPD.clamped_cov & CPD.clamped_weights;
+CPD = set_clamped(CPD, clamped);
+
+%%%%%%%%%%%
+
+function CPD = init_fields()
+% This ensures we define the fields in the same order 
+% no matter whether we load an object from a file,
+% or create it from scratch. (Matlab requires this.)
+
+CPD.self = [];
+CPD.sizes = [];
+CPD.cps = [];
+CPD.dps = [];
+CPD.mean = [];
+CPD.cov = [];
+CPD.weights = [];
+CPD.clamped_mean = [];
+CPD.clamped_cov = [];
+CPD.clamped_weights = [];
+CPD.init_cov = [];
+CPD.cov_type = [];
+CPD.tied_cov = [];
+CPD.Wsum = [];
+CPD.WYsum = [];
+CPD.WXsum = [];
+CPD.WYYsum = [];
+CPD.WXXsum = [];
+CPD.WXYsum = [];
+CPD.nsamples = [];
+CPD.nparams = [];            
+CPD.cov_prior_weight = [];
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m
new file mode 100644
index 00000000..3fa398c8
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/log_prob_node.m
@@ -0,0 +1,59 @@
+function L = log_prob_node(CPD, self_ev, pev)
+% LOG_PROB_NODE Compute prod_m log P(x(i,m)| x(pi_i,m), theta_i) for node i (gaussian)
+% L = log_prob_node(CPD, self_ev, pev)
+%
+% self_ev(m) is the evidence on this node in case m.
+% pev(i,m) is the evidence on the i'th parent in case m (if there are any parents).
+% (These may also be cell arrays.)
+
+if iscell(self_ev), usecell = 1; else usecell = 0; end
+
+use_log = 1;
+ncases = length(self_ev);
+nparents = length(CPD.sizes)-1;
+assert(ncases == size(pev, 2));
+
+if ncases == 0
+  L = 0;
+  return;
+end
+
+if length(CPD.dps)==0 % no discrete parents, so we can vectorize
+  i = 1;
+  if usecell
+    Y = cell2num(self_ev);
+  else
+    Y = self_ev;
+  end
+  if length(CPD.cps) == 0 
+    L = gaussian_prob(Y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
+  else
+    if usecell
+      X = cell2num(pev);
+    else
+      X = pev;
+    end
+    L = gaussian_prob(Y, CPD.mean(:,i) + CPD.weights(:,:,i)*X, CPD.cov(:,:,i), use_log);
+  end
+else % each case uses a (potentially) different set of parameters
+  L = 0;
+  for m=1:ncases
+    if usecell
+      dpvals = cat(1, pev{CPD.dps, m});
+    else
+      dpvals = pev(CPD.dps, m);
+    end
+    i = subv2ind(CPD.sizes(CPD.dps), dpvals(:)');
+    y = self_ev{m};
+    if length(CPD.cps) == 0 
+      L = L + gaussian_prob(y, CPD.mean(:,i), CPD.cov(:,:,i), use_log);
+    else
+      if usecell
+	x = cat(1, pev{CPD.cps, m});
+      else
+	x = pev(CPD.cps, m);
+      end
+      L = L + gaussian_prob(y, CPD.mean(:,i) + CPD.weights(:,:,i)*x, CPD.cov(:,:,i), use_log);
+    end
+  end
+end
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m
new file mode 100644
index 00000000..48447358
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/maximize_params.m
@@ -0,0 +1,147 @@
+function CPD = maximize_params(CPD, temp)
+% MAXIMIZE_PARAMS Set the params of a CPD to their ML values (Gaussian)
+% CPD = maximize_params(CPD, temperature)
+%
+% Temperature is currently only used for entropic prior on Sigma
+
+% For details, see "Fitting a Conditional Gaussian Distribution", Kevin Murphy, tech. report,
+% 1998, available at www.cs.berkeley.edu/~murphyk/papers.html
+% Refering to table 2, we use equations 1/2 to estimate the covariance matrix in the untied/tied case,
+% and equation 9 to estimate the weight matrix and mean.
+% We do not implement spherical Gaussians - the code is already pretty complicated!
+
+if ~adjustable_CPD(CPD), return; end
+
+%assert(approxeq(CPD.nsamples, sum(CPD.Wsum)));
+assert(~any(isnan(CPD.WXXsum)))
+assert(~any(isnan(CPD.WXYsum)))
+assert(~any(isnan(CPD.WYYsum)))
+
+[self_size cpsize dpsize] = size(CPD.weights);
+
+% Append 1s to the parents, and derive the corresponding cross products.
+% This is used when estimate the means and weights simultaneosuly,
+% and when estimatting Sigma.
+% Let x2 = [x 1]'
+XY = zeros(cpsize+1, self_size, dpsize); % XY(:,:,i) = sum_l w(l,i) x2(l) y(l)' 
+XX = zeros(cpsize+1, cpsize+1, dpsize); % XX(:,:,i) = sum_l w(l,i) x2(l) x2(l)' 
+YY = zeros(self_size, self_size, dpsize); % YY(:,:,i) = sum_l w(l,i) y(l) y(l)' 
+for i=1:dpsize
+  XY(:,:,i) = [CPD.WXYsum(:,:,i) % X*Y
+	       CPD.WYsum(:,i)']; % 1*Y
+  % [x  * [x' 1]  = [xx' x
+  %  1]              x'  1]
+  XX(:,:,i) = [CPD.WXXsum(:,:,i) CPD.WXsum(:,i);
+	       CPD.WXsum(:,i)'   CPD.Wsum(i)];
+  YY(:,:,i) = CPD.WYYsum(:,:,i);
+end
+
+w = CPD.Wsum(:);
+% Set any zeros to one before dividing
+% This is valid because w(i)=0 => WYsum(:,i)=0, etc
+w = w + (w==0);
+
+if CPD.clamped_mean
+  % Estimating B2 and then setting the last column (the mean) to the clamped mean is *not* equivalent
+  % to estimating B and then adding the clamped_mean to the last column.
+  if ~CPD.clamped_weights
+    B = zeros(self_size, cpsize, dpsize);
+    for i=1:dpsize
+      if det(CPD.WXXsum(:,:,i))==0
+	B(:,:,i) = 0;
+      else
+	% Eqn 9 in table 2 of TR
+	%B(:,:,i) = CPD.WXYsum(:,:,i)' * inv(CPD.WXXsum(:,:,i));
+	B(:,:,i) = (CPD.WXXsum(:,:,i) \ CPD.WXYsum(:,:,i))';
+      end
+    end
+    %CPD.weights = reshape(B, [self_size cpsize dpsize]);
+    CPD.weights = B;
+  end
+elseif CPD.clamped_weights % KPM 1/25/02
+  if ~CPD.clamped_mean % ML estimate is just sample mean of the residuals
+    for i=1:dpsize
+      CPD.mean(:,i) = (CPD.WYsum(:,i) - CPD.weights(:,:,i) * CPD.WXsum(:,i)) / w(i);
+    end
+  end
+else % nothing is clamped, so estimate mean and weights simultaneously
+  B2 = zeros(self_size, cpsize+1, dpsize);
+  for i=1:dpsize
+    if det(XX(:,:,i))==0  % fix by U. Sondhauss 6/27/99
+      B2(:,:,i)=0;          
+    else                    
+      % Eqn 9 in table 2 of TR
+      %B2(:,:,i) = XY(:,:,i)' * inv(XX(:,:,i));
+      B2(:,:,i) = (XX(:,:,i) \ XY(:,:,i))';
+    end                   
+    CPD.mean(:,i) = B2(:,cpsize+1,i);
+    CPD.weights(:,:,i) = B2(:,1:cpsize,i);
+  end
+end
+
+% Let B2 = [W mu]
+if cpsize>0
+  B2(:,1:cpsize,:) = reshape(CPD.weights, [self_size cpsize dpsize]);
+end
+B2(:,cpsize+1,:) = reshape(CPD.mean, [self_size dpsize]);
+
+% To avoid singular covariance matrices,
+% we use the regularization method suggested in "A Quasi-Bayesian approach to estimating
+% parameters for mixtures of normal distributions", Hamilton 91.
+% If the ML estimate is Sigma = M/N, the MAP estimate is (M+gamma*I) / (N+gamma),
+% where gamma >=0 is a smoothing parameter (equivalent sample size of I prior)
+
+gamma = CPD.cov_prior_weight;
+
+if ~CPD.clamped_cov
+  if CPD.cov_prior_entropic % eqn 12 of Brand AI/Stat 99
+    Z = 1-temp;
+    % When temp > 1, Z is negative, so we are dividing by a smaller
+    % number, ie. increasing the variance.
+  else
+    Z = 0;
+  end
+  if CPD.tied_cov
+    S = zeros(self_size, self_size);
+    % Eqn 2 from table 2 in TR
+    for i=1:dpsize
+      S = S + (YY(:,:,i) - B2(:,:,i)*XY(:,:,i));
+    end
+    %denom = max(1, CPD.nsamples + gamma + Z);
+    denom = CPD.nsamples + gamma + Z;
+    S = (S + gamma*eye(self_size)) / denom;
+    if strcmp(CPD.cov_type, 'diag')
+      S = diag(diag(S));
+    end
+    CPD.cov = repmat(S, [1 1 dpsize]);
+  else 
+    for i=1:dpsize      
+      % Eqn 1 from table 2 in TR
+      S = YY(:,:,i) - B2(:,:,i)*XY(:,:,i);
+      %denom = max(1, w(i) + gamma + Z); % gives wrong answers on mhmm1
+      denom = w(i) + gamma + Z;
+      S = (S + gamma*eye(self_size)) / denom;
+      CPD.cov(:,:,i) = S;
+    end
+    if strcmp(CPD.cov_type, 'diag')
+      for i=1:dpsize      
+	CPD.cov(:,:,i) = diag(diag(CPD.cov(:,:,i)));
+      end
+    end
+  end
+end
+
+
+check_covars = 0;
+min_covar = 1e-5;
+if check_covars % prevent collapsing to a point
+  for i=1:dpsize
+    if min(svd(CPD.cov(:,:,i))) < min_covar
+      disp(['resetting singular covariance for node ' num2str(CPD.self)]);
+      CPD.cov(:,:,i) = CPD.init_cov(:,:,i);
+    end
+  end
+end
+
+
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m
new file mode 100644
index 00000000..988012e2
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_ess.m
@@ -0,0 +1,85 @@
+function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+% UPDATE_ESS Update the Expected Sufficient Statistics of a Gaussian node
+% function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
+
+%if nargin < 6
+%  hidden_bitv = zeros(1, max(fmarginal.domain));
+%  hidden_bitv(find(isempty(evidence)))=1;
+%end
+
+dom = fmarginal.domain;
+self = dom(end);
+ps = dom(1:end-1);
+hidden_self = hidden_bitv(self);
+cps = myintersect(ps, cnodes);
+dps = mysetdiff(ps, cps);
+hidden_cps = all(hidden_bitv(cps));
+hidden_dps = all(hidden_bitv(dps));
+
+CPD.nsamples = CPD.nsamples + 1;            
+[ss cpsz dpsz] = size(CPD.weights); % ss = self size
+
+% Let X be the cts parent (if any), Y be the cts child (self).
+
+if ~hidden_self & (isempty(cps) | ~hidden_cps) & hidden_dps % all cts nodes are observed, all discrete nodes are hidden
+  % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0
+  % Since discrete parents are hidden, we do not need to add evidence to w.
+  w = fmarginal.T(:);
+  CPD.Wsum = CPD.Wsum + w;
+  y = evidence{self};
+  Cyy = y*y';
+  if ~CPD.useC
+     W = repmat(w(:)',ss,1); % W(y,i) = w(i)
+     W2 = repmat(reshape(W, [ss 1 dpsz]), [1 ss 1]); % W2(x,y,i) = w(i)
+     CPD.WYsum = CPD.WYsum +  W .* repmat(y(:), 1, dpsz);
+     CPD.WYYsum = CPD.WYYsum + W2  .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]);
+  else
+     W = w(:)';
+     W2 = reshape(W, [1 1 dpsz]);
+     CPD.WYsum = CPD.WYsum +  rep_mult(W, y(:), size(CPD.WYsum)); 
+     CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum));
+  end
+  if cpsz > 0 % X exists
+    x = cat(1, evidence{cps}); x = x(:);
+    Cxx = x*x';
+    Cxy = x*y';
+    if ~CPD.useC
+       CPD.WXsum = CPD.WXsum + W .* repmat(x(:), 1, dpsz);
+       CPD.WXXsum = CPD.WXXsum + W2 .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]);
+       CPD.WXYsum = CPD.WXYsum + W2 .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]);
+    else
+       CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum));
+       CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum));
+       CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum));
+    end
+  end
+  return;
+end
+
+% general (non-vectorized) case
+fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow!
+
+if dpsz == 1 % no discrete parents
+  w = 1;
+else
+  w = fullm.T(:);
+end
+
+CPD.Wsum = CPD.Wsum + w;
+xi = 1:cpsz;
+yi = (cpsz+1):(cpsz+ss);
+for i=1:dpsz
+  muY = fullm.mu(yi, i);
+  SYY = fullm.Sigma(yi, yi, i);
+  CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
+  CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
+  if cpsz > 0
+    muX = fullm.mu(xi, i);
+    SXX = fullm.Sigma(xi, xi, i);
+    SXY = fullm.Sigma(xi, yi, i);
+    CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
+    CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
+    CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
+  end
+end                
+
diff --git a/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m
new file mode 100644
index 00000000..798c795c
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/CPDs/@gaussian_CPD/Old/update_tied_ess.m
@@ -0,0 +1,118 @@
+function CPD = update_tied_ess(CPD, domain, engine, evidence, ns, cnodes)
+
+if ~adjustable_CPD(CPD), return; end
+nCPDs = size(domain, 2);
+fmarginal = cell(1, nCPDs);
+for l=1:nCPDs
+  fmarginal{l} = marginal_family(engine, nodes(l));
+end
+
+[ss cpsz dpsz] = size(CPD.weights);
+if const_evidence_pattern(engine)
+  dom = domain(:,1);
+  dnodes = mysetdiff(1:length(ns), cnodes);
+  ddom = myintersect(dom, dnodes);
+  cdom = myintersect(dom, cnodes);
+  odom = dom(~isemptycell(evidence(dom)));
+  hdom = dom(isemptycell(evidence(dom)));
+  % If all hidden nodes are discrete and all cts nodes are observed 
+  % (e.g., HMM with Gaussian output)
+  % we can add the observed evidence in parallel
+  if mysubset(ddom, hdom) & mysubset(cdom, odom)
+    [mu, Sigma, T] = add_cts_ev_to_marginals(fmarginal, evidence, ns, cnodes);
+  else
+    mu = zeros(ss, dpsz, nCPDs);
+    Sigma = zeros(ss, ss, dpsz, nCPDs);
+    T = zeros(dpsz, nCPDs);
+    for l=1:nCPDs
+      [mu(:,:,l), Sigma(:,:,:,l), T(:,l)] = add_ev_to_marginals(fmarginal{l}, evidence, ns, cnodes);
+    end
+  end
+end
+CPD.nsamples = CPD.nsamples + nCPDs;            
+
+
+if dpsz == 1 % no discrete parents
+  w = 1;
+else
+  w = fullm.T(:);
+end
+CPD.Wsum = CPD.Wsum + w;
+% Let X be the cts parent (if any), Y be the cts child (self).
+xi = 1:cpsz;
+yi = (cpsz+1):(cpsz+ss);
+for i=1:dpsz
+  muY = fullm.mu(yi, i);
+  SYY = fullm.Sigma(yi, yi, i);
+  CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
+  CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
+  if cpsz > 0
+    muX = fullm.mu(xi, i);
+    SXX = fullm.Sigma(xi, xi, i);
+    SXY = fullm.Sigma(xi, yi, i);
+    CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
+    CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
+    CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
+  end
+end                
+
+
+%%%%%%%%%%%%%
+
+function fullm = add_evidence_to_marginal(fmarginal, evidence, ns, cnodes)
+
+
+dom = fmarginal.domain;
+
+% Find out which values of the discrete parents (if any) are compatible with 
+% the discrete evidence (if any).
+dnodes = mysetdiff(1:length(ns), cnodes);
+ddom = myintersect(dom, dnodes);
+cdom = myintersect(dom, cnodes);
+odom = dom(~isemptycell(evidence(dom)));
+hdom = dom(isemptycell(evidence(dom)));
+
+dobs = myintersect(ddom, odom);
+dvals = cat(1, evidence{dobs});
+ens = ns; % effective node sizes
+ens(dobs) = 1;
+S = prod(ens(ddom));
+subs = ind2subv(ens(ddom), 1:S);
+mask = find_equiv_posns(dobs, ddom);
+subs(mask) = dvals;
+supportedQs = subv2ind(ns(ddom), subs);
+
+if isempty(ddom)
+  Qarity = 1;
+else
+  Qarity = prod(ns(ddom));
+end
+fullm.T = zeros(Qarity, 1);
+fullm.T(supportedQs) = fmarginal.T(:);
+
+% Now put the hidden cts parts into their right blocks,
+% leaving the observed cts parts as 0.
+cobs = myintersect(cdom, odom);
+chid = myintersect(cdom, hdom);
+cvals = cat(1, evidence{cobs});
+n = sum(ns(cdom));
+fullm.mu = zeros(n,Qarity);
+fullm.Sigma = zeros(n,n,Qarity);
+
+if ~isempty(chid)
+  chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom));
+end
+if ~isempty(cobs)
+  cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom));
+end
+
+for i=1:length(supportedQs)
+  Q = supportedQs(i);
+  if ~isempty(chid)
+    fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i);
+    fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i);
+  end
+  if ~isempty(cobs)
+    fullm.mu(cobs_blocks, Q) = cvals(:);
+  end
+end