Fitting Gaussian Mixture Model

Question

I have six bivariate normal distributions and I want to combine them as a Gaussian mixture model. I calculated the mean and covariance matrices below. When I sample random data (mvnrnd) for given distribution parameters, gmdistribution.fit gives different results for different sample sizes. In other words, random sampling sizes n=50 and n=1000 converge different gaussian distributions. My underlying data contains 30 samples for each cluster. So what is the best way to fit gaussian mixture model to my data? Any ideas?

mu1=[log(0.29090) log(0.0038)]
mu2=[log(0.4017) log(0.0053)]
mu3=[log(0.4477) log(0.0051)]
mu4=[log(0.5396) log(0.0072)]
mu5=[log(0.6881) log(0.0090)]
mu6=[log(0.8091) log(0.0099)]

cov1=[0.052 0.0011;0.0011 0.044]
cov2=[0.054 0.0010;0.0010 0.078]
cov3=[0.126 0.011;0.011 0.23]
cov4=[0.092 0.0061;0.0061 0.12]
cov5=[0.113 0.0092;0.0092 0.14]
cov6=[0.1047 0.0217;0.0217 0.35]


X = [mvnrnd(mu1,cov1,50);mvnrnd(mu2,cov2,50);mvnrnd(mu3,cov3,50);mvnrnd(mu4,cov4,50);mvnrnd(mu5,cov5,50);mvnrnd(mu6,cov6,50)];
scatter(X(:,1),X(:,2),'g')

options = statset('MaxIter',200,'Display','final','TolFun',1e-6)


obj = gmdistribution.fit(X,6,'Options',options)

hold on
ezcontour(@(x,y)pdf(obj,[x y]),[-2.5 1],[-7 -2.5],300);
hold off
ezsurfc(@(x,y) pdf(obj,[x y]))

x = -2.5:0.1:1.5; y = -7.0:0.1:-3; n=length(x); a=zeros(n,n);
for i = 1:n,
for j = 1:n,
    gaussPDF(i,j) = pdf(obj,[x(i) y(j)]);

end;
end;