This is the variational Bayesian inference method for Gaussian mixture model. Unlike the EM algorithm (maximum likelihood estimation), it can automatically determine the number of the mixture components k. Please try following code for a demo:
close all; clear;
d = 2;
k = 3;
n = 2000;
[X,z] = mixGaussRnd(d,k,n);
plotClass(X,z);
m = floor(n/2);
X1 = X(:,1:m);
X2 = X(:,(m+1):end);
% VB fitting
[y1, model, L] = mixGaussVb(X1,10);
figure;
plotClass(X1,y1);
figure;
plot(L)
% Predict testing data
[y2, R] = mixGaussVbPred(model,X2);
figure;
plotClass(X2,y2);
The data set is of 3 clusters. You only need to set a number (say 10) which is larger than the intrinsic number of clusters. The algorithm will automatically find the proper k.
Detail description of the algorithm can be found in the reference.
Pattern Recognition and Machine Learning by Christopher M. Bishop (P.474)

Upon the request, I provided the prediction function for out-of-sample inference.

This function is now a part of the PRML toolbox (http://www.mathworks.com/matlabcentral/fileexchange/55826-pattern-recognition-and-machine-learning-toolbox).