comparison toolboxes/FullBNT-1.0.7/netlab3.3/gtmem.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [net, options, errlog] = gtmem(net, t, options)
+%GTMEM	EM algorithm for Generative Topographic Mapping.
+%
+%	Description
+%	[NET, OPTIONS, ERRLOG] = GTMEM(NET, T, OPTIONS) uses the Expectation
+%	Maximization algorithm to estimate the parameters of a GTM defined by
+%	a data structure NET. The matrix T represents the data whose
+%	expectation is maximized, with each row corresponding to a vector.
+%	It is assumed that the latent data NET.X has been set following a
+%	call to GTMINIT, for example.    The optional parameters have the
+%	following interpretations.
+%
+%	OPTIONS(1) is set to 1 to display error values; also logs error
+%	values in the return argument ERRLOG. If OPTIONS(1) is set to 0, then
+%	only warning messages are displayed.  If OPTIONS(1) is -1, then
+%	nothing is displayed.
+%
+%	OPTIONS(3) is a measure of the absolute precision required of the
+%	error function at the solution. If the change in log likelihood
+%	between two steps of the EM algorithm is less than this value, then
+%	the function terminates.
+%
+%	OPTIONS(14) is the maximum number of iterations; default 100.
+%
+%	The optional return value OPTIONS contains the final error value
+%	(i.e. data log likelihood) in OPTIONS(8).
+%
+%	See also
+%	GTM, GTMINIT
+%
+
+%	Copyright (c) Ian T Nabney (1996-2001)
+
+% Check that inputs are consistent
+errstring = consist(net, 'gtm', t);
+if ~isempty(errstring)
+  error(errstring);
+end
+
+% Sort out the options
+if (options(14))
+  niters = options(14);
+else
+  niters = 100;
+end
+
+display = options(1);
+store = 0;
+if (nargout > 2)
+  store = 1;	% Store the error values to return them
+  errlog = zeros(1, niters);
+end
+test = 0;
+if options(3) > 0.0
+  test = 1;	% Test log likelihood for termination
+end
+
+% Calculate various quantities that remain constant during training
+[ndata, tdim] = size(t);
+ND = ndata*tdim;
+[net.gmmnet.centres, Phi] = rbffwd(net.rbfnet, net.X);
+Phi = [Phi ones(size(net.X, 1), 1)];
+PhiT = Phi';
+[K, Mplus1] = size(Phi);
+
+A = zeros(Mplus1, Mplus1);
+cholDcmp = zeros(Mplus1, Mplus1);
+% Use a sparse representation for the weight regularizing matrix.
+if (net.rbfnet.alpha > 0)
+  Alpha = net.rbfnet.alpha*speye(Mplus1);
+  Alpha(Mplus1, Mplus1) = 0;
+end 
+
+for n = 1:niters
+   % Calculate responsibilities
+   [R, act] = gtmpost(net, t);
+     % Calculate error value if needed
+   if (display | store | test)
+      prob = act*(net.gmmnet.priors)';
+      % Error value is negative log likelihood of data
+      e = - sum(log(max(prob,eps)));
+      if store
+         errlog(n) = e;
+      end
+      if display > 0
+         fprintf(1, 'Cycle %4d  Error %11.6f\n', n, e);
+      end
+      if test
+         if (n > 1 & abs(e - eold) < options(3))
+            options(8) = e;
+            return;
+         else
+            eold = e;
+         end
+      end
+   end
+
+   % Calculate matrix be inverted (Phi'*G*Phi + alpha*I in the papers).
+   % Sparse representation of G normally executes faster and saves
+   % memory
+   if (net.rbfnet.alpha > 0)
+      A = full(PhiT*spdiags(sum(R)', 0, K, K)*Phi + ...
+         (Alpha.*net.gmmnet.covars(1)));
+   else
+      A = full(PhiT*spdiags(sum(R)', 0, K, K)*Phi);
+   end
+   % A is a symmetric matrix likely to be positive definite, so try
+   % fast Cholesky decomposition to calculate W, otherwise use SVD.
+   % (PhiT*(R*t)) is computed right-to-left, as R
+   % and t are normally (much) larger than PhiT.
+   [cholDcmp singular] = chol(A);
+   if (singular)
+      if (display)
+         fprintf(1, ...
+            'gtmem: Warning -- M-Step matrix singular, using pinv.\n');
+      end
+      W = pinv(A)*(PhiT*(R'*t));
+   else
+      W = cholDcmp \ (cholDcmp' \ (PhiT*(R'*t)));
+   end
+   % Put new weights into network to calculate responsibilities
+   % net.rbfnet = netunpak(net.rbfnet, W);
+   net.rbfnet.w2 = W(1:net.rbfnet.nhidden, :);
+   net.rbfnet.b2 = W(net.rbfnet.nhidden+1, :);
+   % Calculate new distances
+   d = dist2(t, Phi*W);
+   
+   % Calculate new value for beta
+   net.gmmnet.covars = ones(1, net.gmmnet.ncentres)*(sum(sum(d.*R))/ND);
+end
+
+options(8) = -sum(log(gtmprob(net, t)));
+if (display >= 0)
+  disp(maxitmess);
+end