wolffd@0: function a = gmmactiv(mix, x)
wolffd@0: %GMMACTIV Computes the activations of a Gaussian mixture model.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	This function computes the activations A (i.e. the  probability
wolffd@0: %	P(X|J) of the data conditioned on each component density)  for a
wolffd@0: %	Gaussian mixture model.  For the PPCA model, each activation is the
wolffd@0: %	conditional probability of X given that it is generated by the
wolffd@0: %	component subspace. The data structure MIX defines the mixture model,
wolffd@0: %	while the matrix X contains the data vectors.  Each row of X
wolffd@0: %	represents a single vector.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GMM, GMMPOST, GMMPROB
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check that inputs are consistent
wolffd@0: errstring = consist(mix, 'gmm', x);
wolffd@0: if ~isempty(errstring)
wolffd@0:   error(errstring);
wolffd@0: end
wolffd@0: 
wolffd@0: ndata = size(x, 1);
wolffd@0: a = zeros(ndata, mix.ncentres);  % Preallocate matrix
wolffd@0: 
wolffd@0: switch mix.covar_type
wolffd@0:   
wolffd@0: case 'spherical'
wolffd@0:   % Calculate squared norm matrix, of dimension (ndata, ncentres)
wolffd@0:   n2 = dist2(x, mix.centres);
wolffd@0:   
wolffd@0:   % Calculate width factors
wolffd@0:   wi2 = ones(ndata, 1) * (2 .* mix.covars);
wolffd@0:   normal = (pi .* wi2) .^ (mix.nin/2);
wolffd@0:   
wolffd@0:   % Now compute the activations
wolffd@0:   a = exp(-(n2./wi2))./ normal;
wolffd@0:   
wolffd@0: case 'diag'
wolffd@0:   normal = (2*pi)^(mix.nin/2);
wolffd@0:   s = prod(sqrt(mix.covars), 2);
wolffd@0:   for j = 1:mix.ncentres
wolffd@0:     diffs = x - (ones(ndata, 1) * mix.centres(j, :));
wolffd@0:     a(:, j) = exp(-0.5*sum((diffs.*diffs)./(ones(ndata, 1) * ...
wolffd@0:       mix.covars(j, :)), 2)) ./ (normal*s(j));
wolffd@0:   end
wolffd@0:   
wolffd@0: case 'full'
wolffd@0:   normal = (2*pi)^(mix.nin/2);
wolffd@0:   for j = 1:mix.ncentres
wolffd@0:     diffs = x - (ones(ndata, 1) * mix.centres(j, :));
wolffd@0:     % Use Cholesky decomposition of covariance matrix to speed computation
wolffd@0:     c = chol(mix.covars(:, :, j));
wolffd@0:     temp = diffs/c;
wolffd@0:     a(:, j) = exp(-0.5*sum(temp.*temp, 2))./(normal*prod(diag(c)));
wolffd@0:   end
wolffd@0: case 'ppca'
wolffd@0:   log_normal = mix.nin*log(2*pi);
wolffd@0:   d2 = zeros(ndata, mix.ncentres);
wolffd@0:   logZ = zeros(1, mix.ncentres);
wolffd@0:   for i = 1:mix.ncentres
wolffd@0:     k = 1 - mix.covars(i)./mix.lambda(i, :);
wolffd@0:     logZ(i) = log_normal + mix.nin*log(mix.covars(i)) - ...
wolffd@0:       sum(log(1 - k));
wolffd@0:     diffs = x - ones(ndata, 1)*mix.centres(i, :);
wolffd@0:     proj = diffs*mix.U(:, :, i);
wolffd@0:     d2(:,i) = (sum(diffs.*diffs, 2) - ...
wolffd@0:       sum((proj.*(ones(ndata, 1)*k)).*proj, 2)) / ...
wolffd@0:       mix.covars(i);
wolffd@0:   end
wolffd@0:   a = exp(-0.5*(d2 + ones(ndata, 1)*logZ));
wolffd@0: otherwise
wolffd@0:   error(['Unknown covariance type ', mix.covar_type]);
wolffd@0: end
wolffd@0: