wolffd@0: function mix = netlabgmminit(mix, x, options)
wolffd@0: %GMMINIT Initialises Gaussian mixture model from data
wolffd@0: %(Renamed NETLABGMMINIT in MIRtoolbox to avoid conflict with statistics
wolffd@0: %toolbox)
wolffd@0: %	Description
wolffd@0: %	MIX = GMMINIT(MIX, X, OPTIONS) uses a dataset X to initialise the
wolffd@0: %	parameters of a Gaussian mixture model defined by the data structure
wolffd@0: %	MIX.  The k-means algorithm is used to determine the centres. The
wolffd@0: %	priors are computed from the proportion of examples belonging to each
wolffd@0: %	cluster. The covariance matrices are calculated as the sample
wolffd@0: %	covariance of the points associated with (i.e. closest to) the
wolffd@0: %	corresponding centres. For a mixture of PPCA model, the PPCA
wolffd@0: %	decomposition is calculated for the points closest to a given centre.
wolffd@0: %	This initialisation can be used as the starting point for training
wolffd@0: %	the model using the EM algorithm.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GMM
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: [ndata, xdim] = size(x);
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: % Arbitrary width used if variance collapses to zero: make it 'large' so
wolffd@0: % that centre is responsible for a reasonable number of points.
wolffd@0: GMM_WIDTH = 1.0;
wolffd@0: 
wolffd@0: % Use kmeans algorithm to set centres
wolffd@0: options(5) = 1;	
wolffd@0: [mix.centres, options, post] = netlabkmeans(mix.centres, x, options);
wolffd@0: 
wolffd@0: % Set priors depending on number of points in each cluster
wolffd@0: cluster_sizes = max(sum(post, 1), 1);  % Make sure that no prior is zero
wolffd@0: mix.priors = cluster_sizes/sum(cluster_sizes); % Normalise priors
wolffd@0: 
wolffd@0: switch mix.covar_type
wolffd@0: case 'spherical'
wolffd@0:    if mix.ncentres > 1
wolffd@0:       % Determine widths as distance to nearest centre 
wolffd@0:       % (or a constant if this is zero)
wolffd@0:       cdist = dist2(mix.centres, mix.centres);
wolffd@0:       cdist = cdist + diag(ones(mix.ncentres, 1)*realmax);
wolffd@0:       mix.covars = min(cdist);
wolffd@0:       mix.covars = mix.covars + GMM_WIDTH*(mix.covars < eps);
wolffd@0:    else
wolffd@0:       % Just use variance of all data points averaged over all
wolffd@0:       % dimensions
wolffd@0:       mix.covars = mean(diag(cov(x)));
wolffd@0:    end
wolffd@0:   case 'diag'
wolffd@0:     for j = 1:mix.ncentres
wolffd@0:       % Pick out data points belonging to this centre
wolffd@0:       c = x(find(post(:, j)),:);
wolffd@0:       diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
wolffd@0:       mix.covars(j, :) = sum((diffs.*diffs), 1)/size(c, 1);
wolffd@0:       % Replace small entries by GMM_WIDTH value
wolffd@0:       mix.covars(j, :) = mix.covars(j, :) + GMM_WIDTH.*(mix.covars(j, :)<eps);
wolffd@0:     end
wolffd@0:   case 'full'
wolffd@0:     for j = 1:mix.ncentres
wolffd@0:       % Pick out data points belonging to this centre
wolffd@0:       c = x(find(post(:, j)),:);
wolffd@0:       diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
wolffd@0:       mix.covars(:,:,j) = (diffs'*diffs)/(size(c, 1));
wolffd@0:       % Add GMM_WIDTH*Identity to rank-deficient covariance matrices
wolffd@0:       if rank(mix.covars(:,:,j)) < mix.nin
wolffd@0: 	mix.covars(:,:,j) = mix.covars(:,:,j) + GMM_WIDTH.*eye(mix.nin);
wolffd@0:       end
wolffd@0:     end
wolffd@0:   case 'ppca'
wolffd@0:     for j = 1:mix.ncentres
wolffd@0:       % Pick out data points belonging to this centre
wolffd@0:       c = x(find(post(:,j)),:);
wolffd@0:       diffs = c - (ones(size(c, 1), 1) * mix.centres(j, :));
wolffd@0:       [tempcovars, tempU, templambda] = ...
wolffd@0: 	ppca((diffs'*diffs)/size(c, 1), mix.ppca_dim);
wolffd@0:       if length(templambda) ~= mix.ppca_dim
wolffd@0: 	error('Unable to extract enough components');
wolffd@0:       else 
wolffd@0:         mix.covars(j) = tempcovars;
wolffd@0:         mix.U(:, :, j) = tempU;
wolffd@0:         mix.lambda(j, :) = templambda;
wolffd@0:       end
wolffd@0:     end
wolffd@0:   otherwise
wolffd@0:     error(['Unknown covariance type ', mix.covar_type]);
wolffd@0: end
wolffd@0: