I have a data distribution in a 2D histogram.
I need this histogram to be fitted by a sum of 2D gaussians.

My 2D gaussian is defined with it’s amplitude A, mu and covariance matrix as follows :

def gaussian_2d(X, A, mu, cov):
        diff = X - mu
        res = A * (1/(2*np.pi*np.sqrt(det(cov))))*np.exp(-0.5*diff @ inv(cov) @ diff.transpose())
    return res

And my model function is a sum of two gaussians with one that has known mu and cov matrix :

def model_func(xy, A1, mu1_x, mu1_y, cov1_00, cov1_01, cov1_10, cov1_11, A2):
            mu1 = (mu1_x, mu1_y)
            cov1 = np.array([[cov1_00, cov1_01], [cov1_10, cov1_11]])
            gaussian1 = gaussian_2d(xy, A1, mu1, cov1)

            mu2 = cut_g[0][iter] #extracted from another function
            cov2 = cut_g[1][iter] #extracted from another function
            gaussian2 = gaussian_2d(xy, A2, mu2, cov2)

            return gaussian1 + gaussian2

I have tried using sklearn.mixture.GaussianMixture. But that doesn’t seem to allow me to use the fact that one of my gaussians is partially known.

I have tried using scipy.optimize.curve_fit.
After a lot of pain to try and get it to work because of vector dimension errors because I do not really know how to feed a 2d problem to curve_fit (I add this detail because it might be where my problem originates from), I can get it to run.
But I always get an “optimal parameters not found” error :

RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 1800.

I have a pretty good guess of the initial parameters fed in there and my distribution is very gaussian so I think the algorithm should be capable of fitting it.