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

1 files changed, 203 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/KPMstats/clg_Mstep.m b/sourcecodes/bnt-master/KPMstats/clg_Mstep.m
new file mode 100644
index 00000000..b7da903b
--- /dev/null
+++ b/sourcecodes/bnt-master/KPMstats/clg_Mstep.m
@@ -0,0 +1,203 @@
+function [mu, Sigma, B] = clg_Mstep(w, Y, YY, YTY, X, XX, XY, varargin)
+% MSTEP_CLG Compute ML/MAP estimates for a conditional linear Gaussian
+% [mu, Sigma, B] = Mstep_clg(w, Y, YY, YTY, X, XX, XY, varargin)
+%
+% We fit P(Y|X,Q=i) = N(Y; B_i X + mu_i, Sigma_i) 
+% and w(i,t) = p(M(t)=i|y(t)) = posterior responsibility
+% See www.ai.mit.edu/~murphyk/Papers/learncg.pdf.
+%
+% See process_options for how to specify the input arguments.
+%
+% INPUTS:
+% w(i) = sum_t w(i,t) = responsibilities for each mixture component
+%  If there is only one mixture component (i.e., Q does not exist),
+%  then w(i) = N = nsamples,  and 
+%  all references to i can be replaced by 1.
+% Y(:,i) = sum_t w(i,t) y(:,t) = weighted observations
+% YY(:,:,i) = sum_t w(i,t) y(:,t) y(:,t)' = weighted outer product
+% YTY(i) = sum_t w(i,t) y(:,t)' y(:,t) = weighted inner product
+%   You only need to pass in YTY if Sigma is to be estimated as spherical.
+%
+% In the regression context, we must also pass in the following
+% X(:,i) = sum_t w(i,t) x(:,t) = weighted inputs
+% XX(:,:,i) = sum_t w(i,t) x(:,t) x(:,t)' = weighted outer product
+% XY(i) = sum_t w(i,t) x(:,t) y(:,t)' = weighted outer product
+%
+% Optional inputs (default values in [])
+%
+% 'cov_type' - 'full', 'diag' or 'spherical' ['full']
+% 'tied_cov' - 1 (Sigma) or 0 (Sigma_i) [0]
+% 'clamped_cov' - pass in clamped value, or [] if unclamped [ [] ]
+% 'clamped_mean' - pass in clamped value, or [] if unclamped [ [] ]
+% 'clamped_weights' - pass in clamped value, or [] if unclamped [ [] ]
+% 'cov_prior' - added to Sigma(:,:,i) to ensure psd [0.01*eye(d,d,Q)]
+%
+% If cov is tied, Sigma has size d*d.
+% But diagonal and spherical covariances are represented in full size.
+
+[cov_type, tied_cov, ...
+ clamped_cov, clamped_mean, clamped_weights,  cov_prior, ...
+ xs, ys, post] = ...
+    process_options(varargin, ...
+		    'cov_type', 'full', 'tied_cov', 0,  'clamped_cov', [], 'clamped_mean', [], ...
+		    'clamped_weights', [], 'cov_prior', [], ...
+		    'xs', [], 'ys', [], 'post', []);
+
+[Ysz Q] = size(Y);
+
+if isempty(X) % no regression
+  %B = [];
+  B2 = zeros(Ysz, 1, Q);
+  for i=1:Q
+    B(:,:,i) = B2(:,1:0,i); % make an empty array of size Ysz x 0 x Q
+  end
+  [mu, Sigma] = mixgauss_Mstep(w, Y, YY, YTY, varargin{:});
+  return;
+end
+
+
+N = sum(w);
+if isempty(cov_prior)
+  cov_prior = 0.01*repmat(eye(Ysz,Ysz), [1 1 Q]);
+end
+%YY = YY + cov_prior; % regularize the scatter matrix
+
+% Set any zero weights to one before dividing
+% This is valid because w(i)=0 => Y(:,i)=0, etc
+w = w + (w==0);
+
+Xsz = size(X,1);
+% Append 1 to X to get Z
+ZZ = zeros(Xsz+1, Xsz+1, Q);
+ZY = zeros(Xsz+1, Ysz, Q);
+for i=1:Q
+  ZZ(:,:,i) = [XX(:,:,i)  X(:,i);
+	       X(:,i)'    w(i)];
+  ZY(:,:,i) = [XY(:,:,i);
+	       Y(:,i)'];
+end
+
+
+%%% Estimate mean and regression 
+
+if ~isempty(clamped_weights) & ~isempty(clamped_mean)
+  B = clamped_weights;
+  mu = clamped_mean;
+end
+if ~isempty(clamped_weights) & isempty(clamped_mean)
+  B = clamped_weights;
+  % eqn 5
+  mu = zeros(Ysz, Q);
+  for i=1:Q
+    mu(:,i) = (Y(:,i) - B(:,:,i)*X(:,i)) / w(i);
+  end
+end
+if isempty(clamped_weights) & ~isempty(clamped_mean)
+  mu = clamped_mean;
+  % eqn 3
+  B = zeros(Ysz, Xsz, Q);
+  for i=1:Q
+    tmp = XY(:,:,i)' - mu(:,i)*X(:,i)';
+    %B(:,:,i) = tmp * inv(XX(:,:,i));
+    B(:,:,i) = (XX(:,:,i) \ tmp')';
+  end
+end
+if isempty(clamped_weights) & isempty(clamped_mean)
+  mu = zeros(Ysz, Q);
+  B = zeros(Ysz, Xsz, Q);
+  % Nothing is clamped, so we must estimate B and mu jointly
+  for i=1:Q
+    % eqn 9
+    if rcond(ZZ(:,:,i)) < 1e-10
+      sprintf('clg_Mstep warning: ZZ(:,:,%d) is ill-conditioned', i);
+      % probably because there are too few cases for a high-dimensional input
+      ZZ(:,:,i) = ZZ(:,:,i) + 1e-5*eye(Xsz+1);
+    end
+    %A = ZY(:,:,i)' * inv(ZZ(:,:,i));
+    A = (ZZ(:,:,i) \ ZY(:,:,i))';
+    B(:,:,i) = A(:, 1:Xsz);
+    mu(:,i) = A(:, Xsz+1);
+  end
+end
+
+if ~isempty(clamped_cov)
+  Sigma = clamped_cov;
+  return;
+end
+
+
+%%% Estimate covariance
+
+% Spherical
+if cov_type(1)=='s'
+  if ~tied_cov
+    Sigma = zeros(Ysz, Ysz, Q);
+    for i=1:Q
+      % eqn 16
+      A = [B(:,:,i) mu(:,i)];
+      %s = trace(YTY(i) + A'*A*ZZ(:,:,i) - 2*A*ZY(:,:,i)) / (Ysz*w(i)); % wrong!
+      s = (YTY(i) + trace(A'*A*ZZ(:,:,i)) - trace(2*A*ZY(:,:,i))) / (Ysz*w(i));
+      Sigma(:,:,i) = s*eye(Ysz,Ysz);
+
+      %%%%%%%%%%%%%%%%%%% debug
+      if ~isempty(xs)
+	[nx T] = size(xs);
+	zs = [xs; ones(1,T)];
+	yty = 0;
+	zAAz = 0;
+	yAz = 0;
+	for t=1:T
+	  yty = yty + ys(:,t)'*ys(:,t) * post(i,t);
+	  zAAz = zAAz + zs(:,t)'*A'*A*zs(:,t)*post(i,t);
+	  yAz = yAz + ys(:,t)'*A*zs(:,t)*post(i,t);
+	end
+	assert(approxeq(yty, YTY(i)))
+	assert(approxeq(zAAz, trace(A'*A*ZZ(:,:,i))))
+	assert(approxeq(yAz, trace(A*ZY(:,:,i))))
+	s2 = (yty + zAAz - 2*yAz) / (Ysz*w(i));
+	assert(approxeq(s,s2))
+      end
+      %%%%%%%%%%%%%%% end debug
+      
+    end
+  else
+    S = 0;
+    for i=1:Q
+      % eqn 18
+      A = [B(:,:,i) mu(:,i)];
+      S = S + trace(YTY(i) + A'*A*ZZ(:,:,i) - 2*A*ZY(:,:,i));
+    end
+    Sigma = repmat(S / (N*Ysz), [1 1 Q]);
+  end
+else % Full/diagonal
+  if ~tied_cov
+    Sigma = zeros(Ysz, Ysz, Q);
+    for i=1:Q
+      A = [B(:,:,i) mu(:,i)];
+      % eqn 10
+      SS = (YY(:,:,i) - ZY(:,:,i)'*A' - A*ZY(:,:,i) + A*ZZ(:,:,i)*A') / w(i);
+      if cov_type(1)=='d'
+	Sigma(:,:,i) = diag(diag(SS));
+      else
+	Sigma(:,:,i) = SS;
+      end
+    end
+  else % tied
+    SS = zeros(Ysz, Ysz);
+    for i=1:Q
+      A = [B(:,:,i) mu(:,i)];
+      % eqn 13
+      SS = SS + (YY(:,:,i) - ZY(:,:,i)'*A' - A*ZY(:,:,i) + A*ZZ(:,:,i)*A');
+    end
+    SS = SS / N;
+    if cov_type(1)=='d'
+      Sigma = diag(diag(SS));
+    else
+      Sigma = SS;
+    end
+    Sigma = repmat(Sigma, [1 1 Q]);
+  end
+end
+
+Sigma = Sigma + cov_prior;
+