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