1

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]

Hiddenguy
  • 547
  • 13
  • 33

1 Answers1

2

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...

M Newville
  • 7,486
  • 2
  • 16
  • 29
  • `amp` and `e` are completely correlated in your definition of Lorentzian. Please (and, always) report the full and exact traceback. I don't see how it is possible for Python2.7 to give a RuntimeWarning there. – M Newville Jun 02 '18 at 13:17
  • Well sorry to disappoint, but that's the whole error. Look at my question, there You can find PS of my desktop. – Hiddenguy Jun 02 '18 at 21:04
  • I don't understand how that could be the entire message or how it could be interlaced with the output text and colorized. But then, I don't understand why you post a screenshot as the output text... why not just post the text? It looks like this may be from some Python-ish IDE or environment -- I don't know what that is, but whatever is colorizing the text is probably "interpreting" the error message. Of course, a complete minimal example and complete output are always a good idea. Then again, the actual problem is already identified that `amp` and `e` are completely correlated. – M Newville Jun 02 '18 at 22:50
  • Well this is what I have done, I copied the whole text, but You assumed that I didn't post the whole error. You say that `amp` and `e` are correlated. What should I change then? Note that in my former code, I had the same values as here for `amp` and `e`, but everything worked well, and there were no errors at all. – Hiddenguy Jun 02 '18 at 23:05
  • Well I've found out what is wrong. The problem was in the bounds. When I thrown them away, everything worked just fine, showing me values I was expected. I marked Your answer as the one which solves my problem, because without Your post I wouldn't finish it. So thank You for cooperation :) – Hiddenguy Jun 03 '18 at 09:51