Search code examples
python-2.7curve-fittinglmfit

Joining of curve fitting models


I have this 7 quasi-lorentzian curves which are fitted to my data. enter image description here

and I would like to join them, to make one connected curved line. Do You have any ideas how to do this? I've read about ComposingModel at lmfit documentation, but it's not clear how to do this.

Here is a sample of my code of two fitted curves.

for dataset in [Bxfft]:
    dataset = np.asarray(dataset)
    freqs, psd = signal.welch(dataset, fs=266336/300, window='hamming', nperseg=16192, scaling='spectrum')
    plt.semilogy(freqs[0:-7000], psd[0:-7000]/dataset.size**0, color='r', label='Bx')
    x = freqs[100:-7900]
    y = psd[100:-7900]

    # 8 Hz
    model = Model(lorentzian)
    params = model.make_params(amp=6, cen=5, sig=1, e=0)
    result = model.fit(y, params, x=x)
    final_fit = result.best_fit
    print "8 Hz mode"
    print(result.fit_report(min_correl=0.25))
    plt.plot(x, final_fit, 'k-', linewidth=2)

    # 14 Hz
    x2 = freqs[220:-7780]
    y2 = psd[220:-7780]

    model2 = Model(lorentzian)
    pars2 = model2.make_params(amp=6, cen=10, sig=3, e=0)
    pars2['amp'].value = 6
    result2 = model2.fit(y2, pars2, x=x2)
    final_fit2 = result2.best_fit
    print "14 Hz mode"
    print(result2.fit_report(min_correl=0.25))
    plt.plot(x2, final_fit2, 'k-', linewidth=2)

UPDATE!!!

I've used some hints from user @MNewville, who posted an answer and using his code I got this: enter image description here

So my code is similar to his, but extended with each peak. What I'm struggling now is replacing ready LorentzModel with my own.

The problem is when I do this, the code gives me an error like this.

C:\Python27\lib\site-packages\lmfit\printfuncs.py:153: RuntimeWarning: invalid value encountered in double_scalars [[Model]] spercent = '({0:.2%})'.format(abs(par.stderr/par.value))

About my own model:

    def lorentzian(x, amp, cen, sig, e):
         return (amp*(1-e)) / ((pow((1.0 * x - cen), 2)) + (pow(sig, 2)))

    peak1 = Model(lorentzian, prefix='p1_')
    peak2 = Model(lorentzian, prefix='p2_')
    peak3 = Model(lorentzian, prefix='p3_')

    # make composite by adding (or multiplying, etc) components
    model = peak1 + peak2 + peak3

    # make parameters for the full model, setting initial values
    # using the prefixes
    params = model.make_params(p1_amp=6, p1_cen=8, p1_sig=1, p1_e=0,
                               p2_ampe=16, p2_cen=14, p2_sig=3, p2_e=0,
                               p3_amp=16, p3_cen=21, p3_sig=3, p3_e=0,)

rest of the code is similar like at @MNewville

[![enter image description here][3]][3]


Solution

  • A composite model for 3 Lorentzians would look like this:

    from lmfit import Model, LorentzianModel
    peak1 = LorentzianModel(prefix='p1_')
    peak2 = LorentzianModel(prefix='p2_')
    peak3 = LorentzianModel(prefix='p3_')
    
    # make composite by adding (or multiplying, etc) components
    model = peak1 + peaks2 + peak3
    
    # make parameters for the full model, setting initial values 
    # using the prefixes
    params = model.make_params(p1_amplitude=10, p1_center=8, p1_sigma=3,
                               p2_amplitude=10, p2_center=15, p2_sigma=3,
                               p3_amplitude=10, p3_center=20, p3_sigma=3)
    
    # perhaps set bounds to prevent peaks from swapping or crazy values
    params['p1_amplitude'].min = 0
    params['p2_amplitude'].min = 0
    params['p3_amplitude'].min = 0
    params['p1_sigma'].min = 0
    params['p2_sigma'].min = 0
    params['p3_sigma'].min = 0
    params['p1_center'].min = 2
    params['p1_center'].max = 11
    params['p2_center'].min = 10
    params['p2_center'].max = 18
    params['p3_center'].min = 17
    params['p3_center'].max = 25
    
    # then do a fit over the full data range
    result = model.fit(y, params, x=x)
    

    I think the key parts you were missing were: a) just add models together, and b) use prefix to avoid name collisions of parameters.

    I hope that is enough to get you started...