0

I want to make a global fit of two data sets and plotting the results. I found the answer to "Python and lmfit: How to fit multiple datasets with shared parameters?" answered by @M Newville. The modified code is:

params =  Parameters()
params.add('center_1',    5.0, vary=True)
params.add('amplitude_1', 10.0, vary=True)
params.add('sigma_1',    1.0, vary=True)

params.add('center_2',    8.0, vary=True)
params.add('amplitude_2', 3.0, vary=True)
params.add('sigma_2',    2.0, vary=True)

x = linspace(-10, 10, 101)
data1 = gaussian(x, 5.33, 3.21, 1.51) + random.normal(0, 0.2, x.size)
data2 = gaussian(x, 3.1, 1.21, 1.51) + random.normal(0, 0.2, x.size)
datasets = [data1, data2]
def residual(params, x, datasets):
    model1 = gaussian(x, params['amplitude_1'], params['center_1'], params['sigma_1'])
    model2 = gaussian(x, params['amplitude_2'], params['center_2'], params['sigma_2'])

    resid1 = datasets[0] - model1
    resid2 = datasets[1] - model2
    return np.concatenate((resid1, resid2))

fit = minimize(residual, params, args=(x,), kws={'datasets': datasets})
print(fit_report(fit))

Só, my problem is that I want to present the data and the fit results in a graph, but I don't know how to do it? Please, does anybody can give me a clue?

Edit:

My new code will be:

params =  Parameters()
params.add('cen',    5.0, vary=True)
params.add('amp', 10.0, vary=True)
params.add('sig',    1.0, vary=True)

params.add('pen',    8.0, vary=True)
params.add('inter', 3.0, vary=True)

def reta(x, a, c):
        return a * x + c

x = linspace(-10, 10, 101)
data1 = gaussian(x, 5.33, 3.21, 1.51) + random.normal(0, 0.2, x.size)
data2 = reta(x, 3.1, 1) + random.normal(0, 0.2, x.size)
datasets = [data1, data2]
def residual(params, x, datasets):
    model1 = gaussian(x, params['amp'], params['cen'], params['sig'])
    model2 = reta(x, params['pen'], params['inter'])

    resid1 = datasets[0] - model1
    resid2 = datasets[1] - model2
    return np.concatenate((resid1, resid2))

fit = minimize(residual, params, args=(x,), kws={'datasets': datasets})
print(fit_report(fit))

plt.figure()
plt.plot(x, data1)
plt.plot(x, data2)

If I had only one function, I can plot making this:

plt.plot(x, residual(fit.params) + data1, 'r', label='best fit')

But for two models, I am having problems.

Charlie Vargas Sarmiento
  • 15
  • 1
  • 8
  • 1
    The example at https://stackoverflow.com/questions/20339234/ includes plotting of the results. If that is not sufficient, please clarify the question and show what you have tried that is not doing what you are looking for. – M Newville Apr 02 '20 at 12:42
  • My problem is that I don't know how to replace the fitting results in the functions to make the graph. In your answer, you create the gauss_dataset function to replace the results. I could use it, but I want to have two different models... só I want to know how to replace my fitting results. I will edit my models in the code. – Charlie Vargas Sarmiento Apr 02 '20 at 14:36
  • Well, yeah having a method to calculate the model for each dataset, would make it much simpler, no? You would then just call that method with your `fit.params`. Otherwise, you'd have to do that by hand... – M Newville Apr 02 '20 at 15:57
  • Yes, but the problem will be if the functions have shared parameters. Please, do you know how I can replace the fitting results in my case? Do I need to change the models to show the problem? – Charlie Vargas Sarmiento Apr 02 '20 at 16:37
  • write a function that takes the "current parameters" and calculates the *models* for the two data sets. Then alter your residual function to call that function and use the results to generate the residual. – M Newville Apr 02 '20 at 21:48

0 Answers0