Fitting 2 2D gaussians to data

I'm trying to fit 2 2D Gaussians to some data using the following fit function (I was playing around with the options and lowerbound to see if I could get a better fit hence them being commented out):

%This is the file that fits 2D Gaussians to data (z is the data)
function fitfunction = TwoDimensionGaussianFitter(X, Y, z)
multiple2DGaussians = @(c,XY) c(1)*exp(-(((XY(:,:,1)-c(2)).^2)/(2*(c(3).^2))) - (((XY(:,:,1)-c(4)).^2)/(2*(c(5).^2)))) + c(6)*exp(-(((XY(:,:,1)-c(7)).^2)/(2*(c(8).^2))) - (((XY(:,:,1)-c(9)).^2)/(2*(c(10).^2))));
c0 = [1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1];
XY(:,:,1) = X;
XY(:,:,2) = Y;

%lowerBound = [0,-Inf,-Inf,-Inf,-Inf,0,-Inf,-Inf,-Inf,-Inf, -Inf];
%options = optimoptions('lsqcurvefit', 'FunctionTolerance', 1e-16, 'OptimalityTolerance', 1e-16, 'StepTolerance', 1e-20, 'MaxFunctionEvaluations', 10000);
fitfunction = lsqcurvefit(multiple2DGaussians, c0, XY, z, [], []);
fitfunction

end

To check it was working (which it isn't) I was applying it to generated Gaussians first:

% This script will add two 2D Gaussians

clear
clc
%Mean and sd of first normal distribtion
mu = [0 0];
sigma = [1 0.3; 0.3 1];

%Mean and sd of first normal distribtion
mu2 = [5 0];
sigma2 = [1 0.3; 0.3 1];

%Set up the grid
x1 = -10:0.1:10;
x2 = -10:0.1:10;
[X1,X2] = meshgrid(x1,x2);
X = [X1(:) X2(:)];

y = mvnpdf(X,mu,sigma);
y = reshape(y,length(x2),length(x1));

y2 = mvnpdf(X,mu2,sigma2);
y2 = reshape(y2,length(x2),length(x1));

z = y + y2;

surf(x1,x2,z)
caxis([min(z(:))-0.5*range(z(:)),max(z(:))])
axis([-10 10 -10 10 0 0.4])
xlabel('x1')
ylabel('x2')
zlabel('Probability Density')

fitFunction = TwoDimensionGaussianFitter(X1, X2, z);

figure
fittedGaussian1 = fitFunction(1).*exp(-((X1-fitFunction(2)).^2/(2*(fitFunction(3).^2)) + (X2-fitFunction(4)).^2/(2*(fitFunction(5)^2))));
surf(X1,X2,fittedGaussian1)
figure
fittedGaussian2 = fitFunction(6).*exp(-((X1-fitFunction(7)).^2/(2*(fitFunction(8).^2)) + (X2-fitFunction(9)).^2/(2*(fitFunction(10)^2))));
surf(X1,X2,fittedGaussian2)
title('fittedGaussian2')
figure
fittedGaussian = fittedGaussian1 + fittedGaussian2;
surf(X1,X2,fittedGaussian)

My results are completely out and are not stable at all.

[Akira Agata]

I believe it should be a "fit Gaussian mixture model" problem. How about trying fitgmdist function?

[Andrew McMahon]

fitgmdist (as far as i can see) only fits 2D gaussians and doesn't work with 3D data?