CODEWITHSUNDEEP
pythonscipyscikit-learn
user144464
user144464

Reputation: 113

Python - Using lmfit to fit a function

I am trying to use lmfit, but getting sporadic results for the parameters:

import numpy as np
import scipy.stats as sp
from scipy.optimize import curve_fit
from lmfit import minimize, Parameters, Parameter, report_fit

x = [0.01,
0.02,
0.03,
0.04,
0.05,
0.06,
0.07,
0.08,
0.09,
0.1,
0.2,
0.3,
0.4,
0.5,
0.6,
0.7,
0.8,
0.9,
0.99999999,]

#some paramters to try to estimate
sigma1 = 6
scale1 = np.exp(5)

#generate values to use to fit
y = sp.lognorm(sigma1,scale=scale1).ppf(x)

#function set-up to use in lmfit
def invlognorm(params,x,data):
    sigma = params['sigma']
    mu = params['mu']

    model = sp.lognorm(sigma,scale=mu).ppf(x)
    return model - data

params = Parameters()
params.add('sigma', value= 1,)
params.add('mu', value= 1, )
# do fit, here with leastsq model
result = minimize(invlognorm,params, method = 'least-squares',args=(x, y))

and finally checking the results 

result.params.pretty_print()
Name      Value      Min      Max   Stderr     Vary     Expr Brute_Step
mu        2.161     -inf      inf     None     True     None     None
sigma     6.754     -inf      inf     None     True     None     None

but these are nowhere near the original values.

Any help on what's going on here and how I can fix this so it gives reasonable results?

Upvotes: 0

Views: 1692

Answers (1)

M Newville
M Newville

Reputation: 7862

You will almost certainly need better starting values for the sigma and mu parameters.

The lognorm().ppf() function diverges at x=1, giving huge values which will completely dominate any measure of misfit such as chi-square. In addition, small changes in parameter values will have essentially no effect on the total misfit, and all fitting algorithms will quickly give up. The huge value at x=1 will also make any fit insensitive to the other data. Perhaps you actually meant some other method of lognorm such pdf or cdf?

If not, you may want to fit "in log space" -- fit the log of your data to the log of the model. That will reduce the importance of the x=1 datum.

Also, though it is not the cause of the problem, your fit did not actually use the leastsq method as your comment says. To use the leastsq (Levenberg-Marquardt method) use:

# do fit, here with leastsq model
result = minimize(invlognorm, params, args=(x, y))

To use scipy.optimize.least_squares (which actually use trust region reflective use

# do fit, here with least_squares model
result = minimize(invlognorm,params, method='least_squares', args=(x, y))

(note spelling. Your example used least-squares which is not recognized, and so causes the Nelder-Mead method to be used).

Upvotes: 2

Related Questions