diff options
Diffstat (limited to 'sourcecodes/bnt-master/BNT/examples/static/Zoubin')
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 |