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