Use scipy.optimize.least_squares to fit Gaussian Mixture Model to spectral data

Question

I'm trying to fit the sum of 3 Gaussians to experimental data using scipy.minimize function.

I need your help in properly passing the objective function since that is where error seems to be pointing to. Obviously, I'm new to Python.

The documentation for scipy.minimize does not have clear examples for this case and the questions about GMMs asked here previously were poorly formulated. Please advise...

import numpy
from matplotlib import pyplot
from scipy.optimize import minimize, rosen_der, rosen_hess

x = numpy.arange(0, 55)
y = [-0.00058032, -0.00063992, -0.00057869, -0.00058395, -0.00053528,  -0.0002694, -0.0003716, -0.000284,    0.00104651,  0.00209935,  0.00360213,  0.00502779,  0.00625538,  0.00715873,   0.00753231,  0.00712235,  0.00689089,  0.0061677,  0.00544124,  0.00478251,  0.00487787,  0.00415067,  0.00368579,  0.00370327,  0.00323007,  0.0029862,   0.00250529,  0.00219493,  0.00212242,  0.00209026,  0.0020827,  0.00204044,  0.00218628,  0.00236552,  0.00245056,  0.00282404,  0.0031072,   0.00332862,  0.00351655,  0.00367349,  0.00387923,  0.00395812,  0.00388796,  0.00379902,  0.00369458,  0.00350222,  0.00337815,  0.0032241,  0.00306897,  0.00294152,  0.00276761,  0.00257586,  0.00231613,  0.00211727,  0.00190347]

# experimental data: y
# objective function that is to be minimized: G1 + G2 + G3 - y
def sumGauss(x, y, *args):
    m1, m2, m3, s1, s2, s3, k1, k2, k3 = args
    ret = -y
    ret += k1 * numpy.exp(-(x - m1)**2 / (2 * s1**2))
    ret += k2 * numpy.exp(-(x - m2)**2 / (2 * s2**2))
    ret += k3 * numpy.exp(-(x - m3)**2 / (2 * s3**2))
    return ret

initial_values = [15, 29, 43, 1, 1, 1, 1, 1, 1]

res = minimize(sumGauss(x, y), initial_values, method='trust-exact',
                jac=rosen_der, hess=rosen_hess,
                options={'gtol': 1e-8, 'disp': True})

Here is the error message:

Traceback (most recent call last):
File "fit_gaussian.py", line 64, in <module>    res = minimize(sumGauss(x, y), params, method='trust-exact',
File "fit_gaussian.py", line 41, in sumGauss
m1, m2, m3, s1, s2, s3, k1, k2, k3 = args
ValueError: not enough values to unpack (expected 9, got 0)

Answer 1

There are just a few things to change. If you want to fit data it is better to use least_squares:

from scipy.optimize import least_squares

Make y a numpy array, so you can perform arithmetic operations on it:

y=numpy.array(y)

Do not unpack args, they should remain a single input. Optimisation functions usually work on the 1st input argument, so move it forward:

def sumGauss(args, x, y):
    m1, m2, m3, s1, s2, s3, k1, k2, k3 = args
    ret = -y
    ret += k1 * np.exp(-(x - m1)**2 / (2 * s1**2))
    ret += k2 * np.exp(-(x - m2)**2 / (2 * s2**2))
    ret += k3 * np.exp(-(x - m3)**2 / (2 * s3**2))
    return ret

Then run the optimisation, those jacobians are for different functions so don't use them. Pass the extra params your functions needs as a tuple with keyword args:

res = least_squares(sumGauss, initial_values, method='trf', args=(x,y))

Your optimised args: res.x

Plot, res.fun are the residuals:

import matplotlib.pyplot as plt
plt.plot(y,'o');plt.plot(res.fun+y);plt.plot(res.fun)