wolffd@0: function [data, indices] = mixgauss_sample(mu, Sigma, mixweights, Nsamples)
wolffd@0: %  mixgauss_sample Sample data from a mixture of Gaussians
wolffd@0: % function [data, indices] = mixgauss_sample(mu, Sigma, mixweights, Nsamples)
wolffd@0: %
wolffd@0: % Model is P(X) = sum_k mixweights(k) N(X; mu(:,k), Sigma(:,:,k)) or Sigma(k) for scalar
wolffd@0: % data(:,i) is the i'th sample from P(X)
wolffd@0: % indices(i) is the component from which sample i was drawn
wolffd@0: 
wolffd@0: [D K] = size(mu);
wolffd@0: data = zeros(D, Nsamples);
wolffd@0: indices = sample_discrete(mixweights, 1, Nsamples);
wolffd@0: for k=1:K
wolffd@0:   if ndims(Sigma) < 3
wolffd@0:     sig = Sigma(k);
wolffd@0:   else
wolffd@0:     sig = Sigma(:,:,k);
wolffd@0:   end
wolffd@0:   ndx = find(indices==k);
wolffd@0:   if length(ndx) > 0
wolffd@0:     data(:,ndx) = sample_gaussian(mu(:,k), sig, length(ndx))';
wolffd@0:   end
wolffd@0: end