Covariance matrix is NoneType in lmfit Python3-6

Question

I want to fit a f(x,y) function using lmfit. Dataset is small and there are many fitting parameters (6 points on x-axis, 11 points on y-axis and 16 unconstrained fitting parameters). Using all defaults from Model.fit I cannot obtain covariance matrix and during the fitting process the values of free parameters are not being changed at all.

I tried to change initial values for the parameters. However, when I set the same kind of problem in OriginPro Surface Fitting functionality, the Levenberg-Marquardt algorithm manages to fit data and estimate the errors (although quite large-valued for certain parameters). This means that there has to be some problem with my code. I can't find where the problem lies. I'm not Python master.

The MWE is as below.

import numpy as np
from lmfit import Model, Parameters
import numdifftools # not calling this doesn't change anything

x, y = np.array([226.5, 361.05, 404.41, 589, 632.8, 1013.98]), np.linspace(0,100,11)
X, Y = np.meshgrid(x, y)

Z = np.array([[1.3945, 1.34896, 1.34415, 1.33432, 1.33306, 1.32612],\
[1.39422, 1.3487, 1.34389, 1.33408, 1.33282, 1.32591],\
[1.39336, 1.34795, 1.34315, 1.33336, 1.33211, 1.32524],\
[1.39208, 1.34682, 1.34205, 1.3323, 1.33105, 1.32424],\
[1.39046, 1.3454, 1.34065, 1.33095, 1.32972, 1.32296],\
[1.38854, 1.34373, 1.33901, 1.32937, 1.32814, 1.32145],\
[1.38636, 1.34184, 1.33714, 1.32757, 1.32636, 1.31974],\
[1.38395, 1.33974, 1.33508, 1.32559, 1.32438, 1.31784],\
[1.38132, 1.33746, 1.33284, 1.32342, 1.32223, 1.31576],\
[1.37849, 1.33501, 1.33042, 1.32109, 1.31991, 1.31353],\
[1.37547, 1.33239, 1.32784, 1.31861, 1.31744, 1.31114]])

#This has to be defined beforehand (otherwise parameters names are not defined error)
a1,a2,a3,a4 = 1.3208, -1.2325E-5, -1.8674E-6, 5.0233E-9
b1,b2,b3,b4 = 5208.2413, -0.5179, -2.284E-2, 6.9608E-5
c1,c2,c3,c4 = -2.5551E8, -18341.336, -920, 2.7729
d1,d2,d3,d4 = 9.3495, 2E-3, 3.6733E-5, -1.2932E-7

# Function to fit
def model(x, y, *args):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)

# This is the callable that is passed to Model.fit. M is a (2,N) array
# where N is the total number of data points in Z, which will be ravelled
# to one dimension.
def _model(M, **args):
    x, y = M
    arr = model(x, y, params)
    return arr

# We need to ravel the meshgrids of X, Y points to a pair of 1-D arrays.
xdata = np.vstack((X.ravel(), Y.ravel()))

# Fitting parameters.
fmodel = Model(_model)
params = Parameters()
params.add_many(('a1',1.3208,True,1,np.inf,None,None),\
                ('a2',-1.2325E-5,True,-np.inf,np.inf,None,None),\
                ('a3',-1.8674E-6,True,-np.inf,np.inf,None,None),\
                ('a4',5.0233E-9,True,-np.inf,np.inf,None,None),\
                ('b1',5208.2413,True,-np.inf,np.inf,None,None),\
                ('b2',-0.5179,True,-np.inf,np.inf,None,None),\
                ('b3',-2.284E-2,True,-np.inf,np.inf,None,None),\
                ('b4',6.9608E-5,True,-np.inf,np.inf,None,None),\
                ('c1',-2.5551E8,True,-np.inf,np.inf,None,None),\
                ('c2',-18341.336,True,-np.inf,np.inf,None,None),\
                ('c3',-920,True,-np.inf,np.inf,None,None),\
                ('c4',2.7729,True,-np.inf,np.inf,None,None),\
                ('d1',9.3495,True,-np.inf,np.inf,None,None),\
                ('d2',2E-3,True,-np.inf,np.inf,None,None),\
                ('d3',3.6733E-5,True,-np.inf,np.inf,None,None),\
                ('d4',-1.2932E-7,True,-np.inf,np.inf,None,None))

result = fmodel.fit(Z.ravel(), params, M=xdata)
fit = model(X, Y, result.params)

print(result.covar)

This code results in covariance being NoneType. I expect that it will after all be calculated, because Origin can somehow manage. If it is needed I can provide all parameters from Origin Surface Fitting Parameters.

When plotting Z-fit difference, there is quite large discrepancy for low x-values (not happening in Origin).

that cannot be a minimal working example because `def _model(M, *args)` is not supported by `lmfit.Model`. Your `model()` and `_model()` functions are both in need of re-writing, and you probably do not need `_model()` at all. — M Newville, Oct 10 '19 at 14:23
I obtained this error using lmfit version 0.9.11. The only change I know made is another asterisk and now the definition reads: def _model(M, **args). I edit the code above. The code now reproduces the error i 0.9.14. So what shall I do then? — Radek D, Oct 11 '19 at 13:44

score 1 · Accepted Answer · answered Oct 11 '19 at 23:11

You are not defining your model function in a way that can be used sensibly by lmfit. You have:

def _model(M, **args):
    x, y = M
    arr = model(x, y, params)
    return arr


def model(x, y, *args):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)


model = Model(_model)

Which has a few problems:

args is not used in _model, and params is not defined in the function so will be module-level.
Similarly in model, args is not used and a1, a2, etc will be taken from the module-level (programming) variables and (importantly!!) these will not be updated in the fit.

In short, your model function never sees varying values for the parameters.

lmfit.Model takes the named function arguments and turns those into parameter names. It does not turn **kws or *position_args into parameter names. So I think that what you want to do is write a model function like this:

def model(x, y, a1, a2, a2, a4, b1, b2, b3 ,b3, c1, c2, c3, c4,
          d1, d2, d3, d4):
    return a1+a2*y+a3*np.power(y,2)+a4*np.power(y,3)+\
    (b1+b2*y+b3*np.power(y,2)+b4*np.power(y,3))/np.power(x,2)+\
    (c1+c2*y+c3*np.power(y,2)+c4*np.power(y,3))/np.power(x,4)+\
    (d1+d2*y+d3*np.power(y,2)+d4*np.power(y,3))/np.power(x,6)

Then create a model from that with:

# Note: don't give a function and Model instance the same name!!
my_model = Model(model, independent_vars=('x', 'y'))

With that model defined you can run the fit, and without having to unravel your data (the independent data in lmfit can be of almost any data type, and data arrays can be multi-dimensional):

result = my_model.fit(Z, params, x=X, y=Y)

For what it is worth, making such changes works for me in the sense that the fit runs to completion. The fit still gets stuck with some of the parameters not updating from their initial values, but that is sort of a separate question from the mechanics of setting up and running the fit, and is probably due to polynomials being pretty unstable or poor initial estimates.

As an aside: np.power(y,n) can be spelled y**n and readability counts. Also, numerical stability is sometimes improved with replaced

a + b*x + c*x**2 + d*x**3

with

a + x*(b + x*(c + x*d))

Though I do not know if that would help in your case.

Thank you for your help. This worked in terms of the fitting process. The only parameters that do not change are those for x^6 dependence. The covariance is again being calculated once I removed the x^6 dependence from the equations. If the equation is kept without (a-d)4, this still causes the covariance to be None. Now, the fact that `result.redchi` is very small (order of 10^-7), after covariance is being calculated, leads me to the conclusion that the model itself is overfitting the data, which is obviously something not wanted in general. But this is a separate question. — Radek D, Oct 14 '19 at 08:50

Covariance matrix is NoneType in lmfit Python3-6

1 Answers1