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

12 files changed, 507 insertions, 0 deletions

diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Entries b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Entries
new file mode 100644
index 00000000..31efa304
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Entries
@@ -0,0 +1,10 @@
+/README/1.1.1.1/Wed May 29 15:59:54 2002//
+/csum.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/ffa.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mfa.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mfa_cl.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/mfademo.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/rdiv.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/rprod.m/1.1.1.1/Wed May 29 15:59:54 2002//
+/rsum.m/1.1.1.1/Wed May 29 15:59:54 2002//
+D
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Repository b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Repository
new file mode 100644
index 00000000..15fbcd8e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Repository
@@ -0,0 +1 @@
+FullBNT/BNT/examples/static/Zoubin
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Root b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Root
new file mode 100644
index 00000000..f3bd14a6
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/CVS/Root
@@ -0,0 +1 @@
+:ext:nsaunier@bnt.cvs.sourceforge.net:/cvsroot/bnt
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/README b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/README
new file mode 100644
index 00000000..0fe8214b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/README
@@ -0,0 +1,61 @@
+This software was downloaded from
+   http://www.gatsby.ucl.ac.uk/~zoubin/software.html
+with permission of the author.
+
+
+This software was written by 
+
+Zoubin Ghahramani
+Dept of Computer Science
+University of Toronto
+zoubin@cs.toronto.edu
+
+This software is written in Matlab 4.2c and should run on all platforms
+supporting this version of Matlab. Matlab is a commercial software
+package available from The MathWorks (http://www.mathworks.com/). 
+
+This software is meant for free non-commercial use and distribution. See the
+copyright notice at the bottom of this page.
+
+If you use it, please refer to the accompanying technical report: 
+
+Ghahramani, Z. and Hinton, G.E. (1996) The EM Algorithm for Mixtures
+of Factor Analyzers. University of Toronto Technical Report CRG-TR-96-1. 
+Available at ftp://ftp.cs.toronto.edu/pub/zoubin/tr-96-1.ps.gz  
+
+If you find bugs, or would like to see if I've implemented any
+extensions, please send me email at zoubin@cs.toronto.edu. The
+software is provided "as is", and I cannot guarantee I will be able
+to fix all problems or answer all inquiries. 
+
+See mfademo.m for a demo.
+
+Hope you find it useful. Please send me email if you find it useful
+and I will put you on a mailing list announcing releases of other
+statistical machine learning software in Matlab.
+
+
+----------------------------------------------------------------------
+	Copyright (c) 1996 by Zoubin Ghahramani
+                Toronto, Ontario, Canada. 
+                   All Rights Reserved 
+
+Permission to use, copy, modify, and distribute this software and its
+documentation for non-commercial purposes only is hereby granted
+without fee, provided that the above copyright notice appears in all
+copies and that both the copyright notice and this permission notice
+appear in supporting documentation, and that my name not be used in
+advertising or publicity pertaining to distribution of the software
+without specific, written prior permission. I make no representations
+about the suitability of this software for any purpose. It is provided
+"as is" without express or implied warranty.
+
+I disclaim all warranties with regard to this software, including all
+implied warranties of merchantability and fitness. In no event shall I
+be liable for any special, indirect or consequential damages or any
+damages whatsoever resulting from loss of use, data or profits,
+whether in an action of contract, negligence or other tortious action,
+arising out of or in connection with the use or performance of this
+software.
+
+Zoubin Ghahramani					 Dec 17, 1996
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/csum.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/csum.m
new file mode 100644
index 00000000..2fba6ca5
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/csum.m
@@ -0,0 +1,11 @@
+% column sum
+% function Z=csum(X)
+
+function Z=csum(X)
+
+N=length(X(:,1));
+if (N>1)
+  Z=sum(X);
+else
+  Z=X;
+end;
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/ffa.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/ffa.m
new file mode 100644
index 00000000..e4caa3e0
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/ffa.m
@@ -0,0 +1,75 @@
+% function [L,Ph,LL]=ffa(X,K,cyc,tol);
+% 
+% Fast Maximum Likelihood Factor Analysis using EM
+%
+% X - data matrix
+% K - number of factors
+% cyc - maximum number of cycles of EM (default 100)
+% tol - termination tolerance (prop change in likelihood) (default 0.0001)
+%
+% L - factor loadings 
+% Ph - diagonal uniquenesses matrix
+% LL - log likelihood curve
+%
+% Iterates until a proportional change < tol in the log likelihood 
+% or cyc steps of EM 
+%
+
+function [L,Ph,LL]=ffa(X,K,cyc,tol);
+
+if nargin<4  tol=0.0001; end;
+if nargin<3  cyc=100; end;
+
+N=length(X(:,1));
+D=length(X(1,:));
+tiny=exp(-700);
+
+X=X-ones(N,1)*mean(X);
+XX=X'*X/N;
+diagXX=diag(XX);
+
+randn('seed', 0);
+cX=cov(X);
+scale=det(cX)^(1/D);
+L=randn(D,K)*sqrt(scale/K);
+Ph=diag(cX);
+
+I=eye(K);
+
+lik=0; LL=[];
+
+const=-D/2*log(2*pi);
+
+
+for i=1:cyc;
+
+  %%%% E Step %%%%
+  Phd=diag(1./Ph);
+  LP=Phd*L;
+  MM=Phd-LP*inv(I+L'*LP)*LP';
+  dM=sqrt(det(MM));
+  beta=L'*MM;
+  XXbeta=XX*beta';
+  EZZ=I-beta*L +beta*XXbeta;
+
+  %%%% Compute log likelihood %%%%
+  
+  oldlik=lik;
+  lik=N*const+N*log(dM)-0.5*N*sum(diag(MM*XX));
+  fprintf('cycle %i lik %g \n',i,lik);
+  LL=[LL lik];
+  
+  %%%% M Step %%%%
+
+  L=XXbeta*inv(EZZ);
+  Ph=diagXX-diag(L*XXbeta');
+
+  if (i<=2)    
+    likbase=lik;
+  elseif (lik<oldlik)     
+    disp('VIOLATION');
+  elseif ((lik-likbase)<(1+tol)*(oldlik-likbase)||~isfinite(lik))  
+    break;
+  end;
+
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa.m
new file mode 100644
index 00000000..2060e331
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa.m
@@ -0,0 +1,153 @@
+% function [Lh,Ph,Mu,Pi,LL]=mfa(X,M,K,cyc,tol);
+% 
+% Maximum Likelihood Mixture of Factor Analysis using EM
+%
+% X - data matrix
+% M - number of mixtures (default 1)
+% K - number of factors in each mixture (default 2)
+% cyc - maximum number of cycles of EM (default 100)
+% tol - termination tolerance (prop change in likelihood) (default 0.0001)
+%
+% Lh - factor loadings 
+% Ph - diagonal uniquenesses matrix
+% Mu - mean vectors
+% Pi - priors
+% LL - log likelihood curve
+%
+% Iterates until a proportional change < tol in the log likelihood 
+% or cyc steps of EM 
+
+function [Lh, Ph,  Mu, Pi, LL] = mfa(X,M,K,cyc,tol)
+
+if nargin<5   tol=0.0001; end;
+if nargin<4   cyc=100; end;
+if nargin<3   K=2; end;
+if nargin<2   M=1; end;
+
+N=length(X(:,1));
+D=length(X(1,:));
+tiny=exp(-700);
+
+%rand('state',0);
+
+fprintf('\n');
+
+if (M==1)
+  [Lh,Ph,LL]=ffa(X,K,cyc,tol);
+  Mu=mean(X);
+  Pi=1;
+else
+  if N==1
+    mX = X;
+  else
+    mX=mean(X);
+  end
+  cX=cov(X);
+  scale=det(cX)^(1/D);
+  randn('state',0); 
+  Lh=randn(D*M,K)*sqrt(scale/K);
+  Ph=diag(cX)+tiny;
+  Pi=ones(M,1)/M;
+  %randn('state',0); 
+  Mu=randn(M,D)*sqrtm(cX)+ones(M,1)*mX;
+  oldMu=Mu;
+  I=eye(K);
+
+  lik=0;
+  LL=[];
+
+  H=zeros(N,M); 	% E(w|x) 
+  EZ=zeros(N*M,K);
+  EZZ=zeros(K*M,K);
+  XX=zeros(D*M,D);
+  s=zeros(M,1);
+  const=(2*pi)^(-D/2);
+  %%%%%%%%%%%%%%%%%%%%
+  for i=1:cyc;
+
+    %%%% E Step %%%%
+
+    Phi=1./Ph;
+    Phid=diag(Phi);
+    for k=1:M
+      Lht=Lh((k-1)*D+1:k*D,:);
+      LP=Phid*Lht;
+      MM=Phid-LP*inv(I+Lht'*LP)*LP';
+      dM=sqrt(det(MM));      	
+      Xk=(X-ones(N,1)*Mu(k,:)); 
+      XM=Xk*MM;
+      H(:,k)=const*Pi(k)*dM*exp(-0.5*rsum(XM.*Xk)); 	
+      EZ((k-1)*N+1:k*N,:)=XM*Lht;
+    end;
+    
+    Hsum=rsum(H);
+    oldlik=lik;
+    lik=sum(log(Hsum+(Hsum==0)*exp(-744)));
+
+    Hzero=(Hsum==0); Nz=sum(Hzero); 
+    H(Hzero,:)=tiny*ones(Nz,M)/M; 
+    Hsum(Hzero)=tiny*ones(Nz,1);
+    
+    H=rdiv(H,Hsum); 				
+    s=csum(H);
+    s=s+(s==0)*tiny;
+    s2=sum(s)+tiny;
+    
+    for k=1:M  
+      kD=(k-1)*D+1:k*D;
+      Lht=Lh(kD,:);
+      LP=Phid*Lht;
+      MM=Phid-LP*inv(I+Lht'*LP)*LP';
+      Xk=(X-ones(N,1)*Mu(k,:)); 
+      XX(kD,:)=rprod(Xk,H(:,k))'*Xk/s(k); 
+      beta=Lht'*MM;
+      EZZ((k-1)*K+1:k*K,:)=I-beta*Lht +beta*XX(kD,:)*beta'; 
+    end;
+
+    %%%% log likelihood %%%%
+
+    LL=[LL lik];
+    fprintf('cycle %g   \tlog likelihood %g ',i,lik);
+    
+    if (i<=2)
+      likbase=lik;
+    elseif (lik<oldlik) 
+      fprintf(' violation');
+    elseif ((lik-likbase)<(1 + tol)*(oldlik-likbase)||~isfinite(lik)) 
+      break;
+    end;
+
+    fprintf('\n');
+    
+    %%%% M Step %%%%
+    
+    % means and covariance structure
+    
+    Ph=zeros(D,1);
+    for k=1:M
+      kD=(k-1)*D+1:k*D;
+      kK=(k-1)*K+1:k*K;
+      kN=(k-1)*N+1:k*N;
+
+      T0=rprod(X,H(:,k));
+      T1=T0'*[EZ(kN,:) ones(N,1)];
+      XH=EZ(kN,:)'*H(:,k);
+      T2=inv([s(k)*EZZ(kK,:) XH; XH' s(k)]);
+      T3=T1*T2;
+      Lh(kD,:)=T3(:,1:K);
+      Mu(k,:)=T3(:,K+1)';
+      T4=diag(T0'*X-T3*T1')/s2;
+      Ph=Ph+T4.*(T4>0); 
+    end;
+
+    Phmin=exp(-700);
+    Ph=Ph.*(Ph>Phmin)+(Ph<=Phmin)*Phmin; % to avoid zero variances
+
+    % priors
+    Pi=s'/s2;
+    
+  end;
+  fprintf('\n');
+end;
+
+
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa_cl.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa_cl.m
new file mode 100644
index 00000000..b90bab18
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfa_cl.m
@@ -0,0 +1,54 @@
+% function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);
+% 
+% Calculates log likelihoods of a data set under a mixture of factor
+% analysis model.
+%
+% X - data matrix
+% Lh - factor loadings 
+% Ph - diagonal uniquenesses matrix
+% Mu - mean vectors
+% Pi - priors
+%
+% lik - log likelihood of X 
+% likv - vector of log likelihoods
+% 
+% If 0 or 1 output arguments requested, lik is returned. If 2 output
+% arguments requested, [lik likv] is returned.
+
+function [lik, likv]=mfa_cl(X,Lh,Ph,Mu,Pi);
+
+N=length(X(:,1));
+D=length(X(1,:));
+K=length(Lh(1,:));
+M=length(Pi);
+
+if (abs(sum(Pi)-1) > 1e-6) 
+  disp('ERROR: Pi should sum to 1');
+  return;
+elseif ((size(Lh) ~= [D*M K]) | (size(Ph) ~= [D 1]) | (size(Mu) ~= [M D]) ...
+  | (size(Pi) ~= [M 1] & size(Pi) ~= [1 M]))   
+  disp('ERROR in input matrix sizes');
+  return;
+end;  
+
+tiny=exp(-744);
+const=(2*pi)^(-D/2);
+
+I=eye(K);
+Phi=1./Ph;
+Phid=diag(Phi);
+for k=1:M  
+  Lht=Lh((k-1)*D+1:k*D,:);
+  LP=Phid*Lht;
+  MM=Phid-LP*inv(I+Lht'*LP)*LP';
+  dM=sqrt(det(MM));      	
+  Xk=(X-ones(N,1)*Mu(k,:)); 
+  XM=Xk*MM; 
+  H(:,k)=const*Pi(k)*dM*exp(-0.5*sum((XM.*Xk)'))'; 	
+end;
+
+Hsum=rsum(H); 				
+
+likv=log(Hsum+(Hsum==0)*tiny);
+lik=sum(likv);
+
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfademo.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfademo.m
new file mode 100644
index 00000000..508d840b
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/mfademo.m
@@ -0,0 +1,81 @@
+echo on;
+
+clc;
+
+% This is a very basic demo of the mixture of factor analyzer software
+% written in Matlab by	Zoubin Ghahramani
+%			Dept of Computer Science
+%			University of Toronto
+
+pause;		% Hit any key to continue 
+
+% To demonstrate the software we generate a sample data set
+% from a mixture of two Gaussians
+
+pause;		% Hit any key to continue 
+
+X1=randn(300,5);	% zero mean 5 dim Gaussian data 
+X2=randn(200,5)+2;	% 5 dim Gaussian data with mean [1 1 1 1 1]
+X=[X1;X2];		% total 500 data points from mixture
+
+% Fitting the model is very easy. For example to fit a mixture of 2
+% factor analyzers with three factors each...
+
+pause;		% Hit any key to continue 
+
+
+[Lh,Ph,Mu,Pi,LL]=mfa(X,2,3);
+
+% Lh, Ph, Mu, and Pi are the factor loadings, observervation
+% variances, observation means for each mixture, and mixing
+% proportions. LL is the vector of log likelihoods (the learning
+% curve). For more information type: help mfa
+
+% to plot the learning curve (log likelihood at each step of EM)...
+
+pause;		% Hit any key to continue 
+
+plot(LL);
+
+% you get a more informative picture of convergence by looking at the
+% log of the first difference of the log likelihoods...
+
+pause;		% Hit any key to continue 
+
+semilogy(diff(LL)); 
+
+% you can look at some of the parameters of the fitted model... 
+
+pause;		% Hit any key to continue 
+
+Mu
+
+Pi
+
+% ...to see whether they make any sense given that me know how the
+% data was generated. 
+
+% you can also evaluate the log likelihood of another data set under
+% the model we have just fitted using the mfa_cl (for Calculate
+% Likelihood) function. For example, here we generate a test from the
+% same distribution. 
+
+
+X1=randn(300,5);
+X2=randn(200,5)+2;
+Xtest=[X1; X2];
+
+pause;		% Hit any key to continue 
+
+mfa_cl(Xtest,Lh,Ph,Mu,Pi)
+
+% we should expect the log likelihood of the test set to be lower than
+% that of the training set.
+
+% finally, we can also fit a regular factor analyzer using the ffa
+% function (Fast Factor Analysis)...
+
+pause;		% Hit any key to continue 
+
+[L,Ph,LL]=ffa(X,3);
+  
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rdiv.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rdiv.m
new file mode 100644
index 00000000..3128061e
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rdiv.m
@@ -0,0 +1,25 @@
+% function Z=rdiv(X,Y)
+%
+% row division: Z = X / Y row-wise
+% Y must have one column 
+
+function Z=rdiv(X,Y)
+
+[N M]=size(X);
+[K L]=size(Y);
+if(N ~= K | L ~=1)
+  disp('Error in RDIV');
+  return;
+end
+
+Z=zeros(N,M);
+
+if M<N,
+  for m=1:M
+    Z(:,m)=X(:,m)./Y;
+  end
+else
+  for n=1:N
+    Z(n,:)=X(n,:)/Y(n);
+  end;
+end;
\ No newline at end of file
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rprod.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rprod.m
new file mode 100644
index 00000000..95d3565d
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rprod.m
@@ -0,0 +1,15 @@
+% row product
+% function Z=rprod(X,Y)
+
+function Z=rprod(X,Y)
+
+if(length(X(:,1)) ~= length(Y(:,1)) | length(Y(1,:)) ~=1)
+  disp('Error in RPROD');
+  return;
+end
+
+Z=zeros(size(X));
+
+for i=1:length(X(1,:))
+  Z(:,i)=X(:,i).*Y;
+end
diff --git a/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rsum.m b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rsum.m
new file mode 100644
index 00000000..0af53fde
--- /dev/null
+++ b/sourcecodes/bnt-master/BNT/examples/static/Zoubin/rsum.m
@@ -0,0 +1,20 @@
+% row sum
+% function Z=rsum(X)
+
+function Z=rsum(X)
+
+[N M]=size(X);
+
+Z=zeros(N,1);
+
+if M==1,
+  Z=X;
+elseif M<2*N,
+  for m=1:M,
+    Z=Z+X(:,m);
+  end;
+else
+  for n=1:N
+    Z(n)=sum(X(n,:));
+  end;
+end