wolffd@0: function [data, label] = gmmsamp(mix, n)
wolffd@0: %GMMSAMP Sample from a Gaussian mixture distribution.
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %
wolffd@0: %	DATA = GSAMP(MIX, N) generates a sample of size N from a Gaussian
wolffd@0: %	mixture distribution defined by the MIX data structure. The matrix X
wolffd@0: %	has N rows in which each row represents a MIX.NIN-dimensional sample
wolffd@0: %	vector.
wolffd@0: %
wolffd@0: %	[DATA, LABEL] = GMMSAMP(MIX, N) also returns a column vector of
wolffd@0: %	classes (as an index 1..N) LABEL.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	GSAMP, GMM
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: % Check input arguments
wolffd@0: errstring = consist(mix, 'gmm');
wolffd@0: if ~isempty(errstring)
wolffd@0:   error(errstring);
wolffd@0: end
wolffd@0: if n < 1
wolffd@0:   error('Number of data points must be positive')
wolffd@0: end
wolffd@0: 
wolffd@0: % Determine number to sample from each component
wolffd@0: priors = rand(1, n);
wolffd@0: 
wolffd@0: % Pre-allocate data array
wolffd@0: data = zeros(n, mix.nin);
wolffd@0: if nargout > 1
wolffd@0:   label = zeros(n, 1);
wolffd@0: end
wolffd@0: cum_prior = 0;		% Cumulative sum of priors
wolffd@0: total_samples = 0;	% Cumulative sum of number of sampled points
wolffd@0: for j = 1:mix.ncentres
wolffd@0:   num_samples = sum(priors >= cum_prior & ...
wolffd@0:     priors < cum_prior + mix.priors(j));
wolffd@0:   % Form a full covariance matrix
wolffd@0:   switch mix.covar_type
wolffd@0:     case 'spherical'
wolffd@0:       covar = mix.covars(j) * eye(mix.nin);
wolffd@0:     case 'diag'
wolffd@0:       covar = diag(mix.covars(j, :));
wolffd@0:     case 'full'
wolffd@0:       covar = mix.covars(:, :, j);
wolffd@0:     case 'ppca'
wolffd@0:       covar = mix.covars(j) * eye(mix.nin) + ...
wolffd@0:         mix.U(:, :, j)* ...
wolffd@0:         (diag(mix.lambda(j, :))-(mix.covars(j)*eye(mix.ppca_dim)))* ...
wolffd@0:         (mix.U(:, :, j)');
wolffd@0:     otherwise
wolffd@0:       error(['Unknown covariance type ', mix.covar_type]);
wolffd@0:   end
wolffd@0:   data(total_samples+1:total_samples+num_samples, :) = ...
wolffd@0:     gsamp(mix.centres(j, :), covar, num_samples);
wolffd@0:   if nargout > 1
wolffd@0:     label(total_samples+1:total_samples+num_samples) = j;
wolffd@0:   end
wolffd@0:   cum_prior = cum_prior + mix.priors(j);
wolffd@0:   total_samples = total_samples + num_samples;
wolffd@0: end
wolffd@0: