Let's say that we have two different curves with the same variables. How can I optimize the variables to have the best fit for both curves at the same time? I can optimize each curve separately but the obtained values will not provide the best fit for the other curve necessarily.
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from lmfit import Model, Parameters
def f(wdata, a, b, c, d):
return ( a+b*wdata + c*d*wdata )
def g(wdata, a, b, c, d):
return ( a**2*(wdata)**2/a*b*wdata*(c)**2 )
wdata = (1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100)
wdata= np.asarray(wdata)
ydata1 = f(wdata, -19, -60, 9, 100)
ydata2 = g(wdata, -19, -60, 9, 100)
fmodel = Model(f)
gmodel = Model(g)
params = Parameters()
params.add('a', value=-19, max = 0, vary=False)
params.add('b', value=-60, vary=True)
params.add('c', value=9, min = 0, vary=True)
params.add('xangle', value=0.05, vary=True, min=-np.pi/2, max=np.pi/2)
params.add('d', expr='(c*a/b)*sin(xangle)')
resultf = fmodel.fit(ydata1, params, wdata=wdata, method='leastsq')
print(resultf.fit_report())
plt.plot(wdata, ydata1, 'bo')
plt.plot(wdata, resultf.init_fit, 'k--')
plt.plot(wdata, resultf.best_fit, 'r-')
resultg = gmodel.fit(ydata2, params, wdata=wdata, method='leastsq')
print(resultg.fit_report())
plt.plot(wdata, ydata2, 'bo')
plt.plot(wdata, resultg.init_fit, 'k--')
plt.plot(wdata, resultg.best_fit, 'r-')
plt.show()
I think what you want to do is create a model function that models the concatenatation the arrays of f and g data. Perhaps something like this:
def f_and_g(wdata, a, b, c, d):
fmodel = f(wdata, a, b, c, d)
gmodel = g(wdata, a, b, c, d)
return np.concatenate((fmodel, gmodel))
# turn that into a model
fg_model = Model(f_and_g)
# same parameters as before:
params = Parameters()
params.add('a', value=-19, max = 0, vary=False)
params.add('b', value=-60, vary=True)
params.add('c', value=9, min = 0, vary=True)
params.add('xangle', value=0.05, vary=True, min=-np.pi/2, max=np.pi/2)
params.add('d', expr='(c*a/b)*sin(xangle)')
# concatenate data to build a 1D array to be modeled:
fdata = f(wdata, -19, -60, 9, 100) + np.normal(scale=0.1, size=len(wdata))
gdata = g(wdata, -19, -60, 9, 100) + np.normal(scale=0.1, size=len(wdata))
fg_data = np.concatenate((fdata, gdata))
# now fit the concatenated data to the concatenated model
result = fg_model.fit(fg_data, params, wdata=wdata, method='leastsq')
print(result.fit_report())
# to plot individual results, you'll have to untangle the concatenated data:
plt.plot(wdata, fdata, 'bo', label='f data')
plt.plot(wdata, result.best_fit[:len(wdata)], 'b--', label='f fit')
plt.plot(wdata, rgdata, 'ro', label='g data')
plt.plot(wdata, result.best_fit[len(wdata):], 'r--', label='g fit')
plt.legend()
plt.show()
One thing that is sort of assumed here is that the f data and g data are equally scaled and important. That is not always the case. When it is not the case, you may want to add a weighting to the Model (see lmfit docs and examples). Typically, such weighting would reflect variations in the uncertainties at each data point (ie, weight=1./stderr, where stderr is the standard error in the data). But for multi-data-set fitting, you may want to alter that so that one data set (or even one part of one data set) is emphasized more than other parts.