lmfit minimize weighting initial data

Question

I am fairly new to python and I am trying to do some curve fitting using lmfit. The code works pretty well but I would like to force the fit through the origin (0,0). I've seen here in stakoverlow that using "curve_fit" you can add and attribute "sigma" that does the trick but this does not work in "minimize". Do you have any workaround to this???

So far, this is my code:

##############################################################################
###################### IMPORT MODULES ########################################
##############################################################################
import numpy as np
import matplotlib.pyplot as plt
from lmfit import minimize, Parameters
from lmfit import report_fit

##############################################################################
##################### DATA ##################################
##############################################################################
x=np.array([0,1, 4, 6, 8, 9, 11, 13, 15, 18, 20, 27])
data=np.array([0, 67, 208, 339, 460, 539, 599, 635, 659, 685, 701, 731])

##############################################################################
######################## GOMPERTZ FUNCTION DEFINITION#########################
##############################################################################

def BMP_Gompertz(params, x, data):
    BMPmax=params['BMPmax'].value
    Rmax=params['Rmax'].value
    Lambda=params['Lambda'].value

    model=BMPmax*np.exp(-np.exp(Rmax*np.exp(1)/BMPmax*(Lambda-x)+1))

    global model_data
    model_data=[BMPmax,Rmax,Lambda]

    return data-model

params=Parameters()
params.add('BMPmax', value=300., min=0.)
params.add('Rmax', value=25., min=0.) 
params.add('Lambda',value=0.5, min=0.)

model_data=[]


##############################################################################
###################### GOMPERTZ CURVE FITTING & REPORT #######################
##############################################################################
result=minimize(BMP_Gompertz, params, args=(x,data))
report_fit(result.params)

##############################################################################
########################## GOMPERTZ PLOTTING #################################
##############################################################################
plt.plot(x, data, 'bo')
t=np.linspace(0.,100,10000)
plt.plot(t,model_data[0]*np.exp(-np.exp(model_data[1]*np.exp(1)/model_data[0]*(model_data[2]-t)+1)))
plt.xlabel('Days')
plt.ylabel('BMP (NLbiogas/kg SV)')
plt.show()

as you can seee, the fit does not start at (0,0) but at around (0,10) and I would like to always force it to start at (0,0)... It looks like I am still not able to upload images (don't have persmissions yet)... anyway, I think you can get the point.

Also another question, is there another way to store the parameters and plot the results?. Right now to store the the parameters returned by the model I am doing it by storing them ina global variable called "model_data". Then, in the plot section I've created a new "x" array called "t" using linspace and then plotting the t vs BMP_Gompertz (myfunction) using the parameters stored in "model_data". Works perfectly, but it looks like it should be other nicer ways to do that.

Thank you very much for your help

Answer 1

I'm not sure all your different questions can be fully answered, but here are a few comments:

1. You can add weights to the fit. With your example using minimize you could pass in an array sigma holding uncertainties in the data, and change return data-model to return (data-model)/sigma. Below, I will recommend using lmfit.Model which has a slightly different way to specify weights.

2. Getting your current model to go through (0, 0) will be challenging, even using weighting. That is, exponential functions decay, but do not ever reach 0 so you may need to decide what is "zero enough". If this is a physical requirement, you may need to modify the model.

3. The fit results are stored, and you do not need to use the model_data. The returned result.params contains the best fit parameters, and you can get result.params['Rmax'].value etc.

Since you are doing curve-fitting and using lmfit, I suggest using lmfit.Model which will simplify your code and make it easier to extract the predicted model for plotting. With lmfit.Model much of the work you do in your objective function is done for you, and your script would become:

import numpy as np
import matplotlib.pyplot as plt
from lmfit import Model

x=np.array([0,1, 4, 6, 8, 9, 11, 13, 15, 18, 20, 27])
data=np.array([0, 67, 208, 339, 460, 539, 599, 635, 659, 685, 701, 731])

# This function now returns the model. The function arguments are 
# named so that Model() will know what to name the Parameters.
def BMP_Gompertz(x, bmpmax, rmax, xlambda):
    return bmpmax *np.exp(-np.exp((rmax*np.exp(1)/bmpmax)*(xlambda-x)+1))

# create a Model from the model function
bmodel = Model(BMP_Gompertz)

# create Parameters, giving initial values
params = bmodel.make_params(bmpmax=300, rmax=25, xlambda=0.5)
params['bmpmax'].min = 0
params['rmax'].min = 0
params['xlambda'].min = 0

# do fit, print result
result = bmodel.fit(data, params, x=x)
print(result.fit_report())

# plot results -- note that `best_fit` is already available
plt.plot(x, data, 'bo')
plt.plot(x, result.best_fit, 'r--')

t=np.linspace(0.,100,10000)

# use `result.eval()` to evaluate model given params and x
plt.plot(t, bmodel.eval(result.params, x=t), 'k-')
plt.xlabel('Days')
plt.ylabel('BMP (NLbiogas/kg SV)')
plt.show()

To add weights to a Model fit, you can define an array of weights to be used to multiply data-model, and pass that to Model.fit. To emphasize the data with the lowest values (and so push the fit towards (0, 0), you might use something like this:

weights = 1.0 / np.sqrt(data + 1)
result = bmodel.fit(data, params, x=x, weights=weights)

Again, that will emphasize the smaller y values, and tend to push down the model value at x=0, but will not really force the result to exactly (0, 0).