Gaussian Fit using lmfit

Question

I am trying to fit a gaussian. I tried to fit using OriginPro and Python. The fit in OriginPro is better than that obtained through Python and I would like to do it using Python.

OriginPro:

Python:

The code I used is:

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

data=data
mod=lmfit.models.GaussianModel()
x=np.arange(140,510,1)
y= np.array(np.log(data[140:510]), dtype=np.float32)
idx= (np.isfinite(y))
pars = mod.guess(y[idx], x=x[idx])
out = mod.fit(y[idx], pars, x=x[idx])
plt.plot(x[idx], out.best_fit, 'r', label='Gaussian Model')
plt.plot(x,np.log(data[140:510]))
plt.show()

data:

[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
        0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., 34., 33., 30.,
       25., 47., 36., 37., 33., 41., 42., 40., 49., 47., 44., 35.,  0.,
        0.,  0., 30., 32., 37., 36., 50., 40., 30., 41., 29., 33., 34.,
       42., 41., 39., 43., 34., 30., 27., 48., 38., 20., 53., 38., 32.,
       41., 40., 45., 39., 34., 33., 41., 32., 41., 41., 30., 42., 41.,
       33., 62., 36., 47., 44., 35., 44., 36., 40., 35., 43., 38., 36.,
       48., 50., 31., 46., 50., 43., 50., 40., 35., 39., 50., 47., 46.,
       35., 39., 33., 49., 43., 56., 45., 56.,  0.,  0.,  0.,  0., 42.,
       41., 43., 46., 40., 40., 38., 40., 39., 49., 33., 39., 57., 51.,
       51., 43., 47., 55., 67., 68., 44., 48.,  0.,  0.,  0.,  0.,  0.,
       48., 66., 55., 60., 63., 58., 72., 53., 62., 55., 63., 64., 60.,
       61., 56., 58., 71., 51., 74., 58., 69., 61., 50., 59., 66., 57.,
       56., 68., 63., 53., 70., 73., 49., 65., 56., 62., 45., 55.,  0.,
        0.,  0.,  0., 66., 53., 52., 52., 61., 42., 53., 49., 81., 67.,
       63., 69., 58., 80., 74., 51., 55., 65., 68., 56., 71., 58., 59.,
       84., 54., 49., 46., 69., 57., 54., 74., 68., 60., 56., 51., 58.,
       47., 57.,  0.,  0.,  0.,  0., 65., 47., 59., 39., 34., 57., 45.,
       43., 57., 61., 39., 52., 48., 42., 53., 41.,  0.,  0.,  0.,  0.,
       48., 35., 53., 44., 46., 52., 50., 51., 44., 55., 42., 58., 53.,
       42.,  0.,  0.,  0.,  0., 41., 39., 47., 46., 46., 46., 41., 44.,
       49., 48., 47., 43., 46., 38., 47., 39., 29., 53., 40., 45., 43.,
       41., 23., 41., 39., 46., 37., 42., 24., 40., 46., 39., 37., 49.,
       47., 37., 28., 36., 34., 34., 43., 41., 48., 42., 44., 45., 43.,
       45., 36., 33., 48., 39., 54., 35., 30., 42., 34., 51., 39., 50.,
       36., 34., 50., 37., 48., 40., 39., 21., 41., 40., 42., 44., 48.,
       42., 42., 39., 43., 57., 37., 36., 38., 34., 50., 40., 30., 29.,
       38., 44., 43., 47., 35., 37., 44., 44., 32., 36., 41., 44., 28.,
       43., 36., 37., 34., 37.,  0.,  0.,  0.,  0., 54., 34., 45.,  0.,
        0.,  0.,  0.,  0., 44., 32., 28., 34., 42., 35., 40., 45., 42.,
       47., 39., 35., 52., 34., 55., 38.,  1.]

Output for fit_report:

In [37]:out.fit_report()
Out[37]: "[[Model]]\n    Model(gaussian)\n[[Fit Statistics]]\n    # fitting method   = leastsq\n    # function evals   = 17\n    # data points      = 333\n    # variables        = 3\n    chi-square         = 12.5923742\n    reduced chi-square = 0.03815871\n    Akaike info crit   = -1084.59201\n    Bayesian info crit = -1073.16758\n[[Variables]]\n    amplitude:  3590.97370 +/- 125.195985 (3.49%) (init = 781.4034)\n    center:     317.738344 +/- 3.57625327 (1.13%) (init = 322.4771)\n    sigma:      361.440601 +/- 13.6760525 (3.78%) (init = 181.5)\n    fwhm:       851.127556 +/- 32.2046420 (3.78%) == '2.3548200*sigma'\n    height:     3.96355944 +/- 0.01657428 (0.42%) == '0.3989423*amplitude/max(1e-15, sigma)'\n[[Correlations]] (unreported correlations are < 0.100)\n    C(amplitude, sigma) =  0.997"

Has it got something to do with the guess function I used?

See https://math.stackexchange.com/q/1292889/295281. I've implemented this method in a package called scikit-guess, which you can find on pypi. It's still experimental — Mad Physicist, May 14 '21 at 18:07
your code is not complete enough to comment on what might have gone wrong. provide a complete, minimal example and show the results (from fit_report()) as well a plot. From a quick look at the images, it appears your OriginPro code has included a constant offset as well as a Gaussian, so you might want to include that in the Python code as well. — M Newville, May 14 '21 at 20:50
@MNewville Now I have added all the code that is relevant, do you want me to show the data in 2theta as well? I have also added the fit_report() and the plot was already there. — Pranav45, May 15 '21 at 04:40
@Pranav45 that cannot be all the code. You left out imports, though those are probably obvious. More importantly, `theta_up` is undefined. As a result, no one can diagnose what your code does. See https://stackoverflow.com/help/minimal-reproducible-example. But also: your model does not include a constant (perhaps varying?) offset. Doesn't your data need that? — M Newville, May 15 '21 at 14:40
@MNewville Thank you for letting me know. Is it alright now? — Pranav45, May 15 '21 at 18:45

score 0 · Accepted Answer · answered May 16 '21 at 02:06

Hints give in the comments seem to be insufficient. You need to include an offset in your model. A Gaussian function goes to 0 far from the peak intensity - I don't know what OriginPro is doing, but clearly, it is modeling more than a Gaussian.

Try making a model that is a Gaussian + a Constant, as with:

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

data = np.array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
              34., 33., 30., 25., 47., 36., 37., 33., 41., 42., 40.,
              49., 47., 44., 35., 0., 0., 0., 30., 32., 37., 36., 50.,
              40., 30., 41., 29., 33., 34., 42., 41., 39., 43., 34.,
              30., 27., 48., 38., 20., 53., 38., 32., 41., 40., 45.,
              39., 34., 33., 41., 32., 41., 41., 30., 42., 41., 33.,
              62., 36., 47., 44., 35., 44., 36., 40., 35., 43., 38.,
              36., 48., 50., 31., 46., 50., 43., 50., 40., 35., 39.,
              50., 47., 46., 35., 39., 33., 49., 43., 56., 45., 56.,
              0., 0., 0., 0., 42., 41., 43., 46., 40., 40., 38., 40.,
              39., 49., 33., 39., 57., 51., 51., 43., 47., 55., 67.,
              68., 44., 48., 0., 0., 0., 0., 0., 48., 66., 55., 60.,
              63., 58., 72., 53., 62., 55., 63., 64., 60., 61., 56.,
              58., 71., 51., 74., 58., 69., 61., 50., 59., 66., 57.,
              56., 68., 63., 53., 70., 73., 49., 65., 56., 62., 45.,
              55., 0., 0., 0., 0., 66., 53., 52., 52., 61., 42., 53.,
              49., 81., 67., 63., 69., 58., 80., 74., 51., 55., 65.,
              68., 56., 71., 58., 59., 84., 54., 49., 46., 69., 57.,
              54., 74., 68., 60., 56., 51., 58., 47., 57., 0., 0., 0.,
              0., 65., 47., 59., 39., 34., 57., 45., 43., 57., 61.,
              39., 52., 48., 42., 53., 41., 0., 0., 0., 0., 48., 35.,
              53., 44., 46., 52., 50., 51., 44., 55., 42., 58., 53.,
              42., 0., 0., 0., 0., 41., 39., 47., 46., 46., 46., 41.,
              44., 49., 48., 47., 43., 46., 38., 47., 39., 29., 53.,
              40., 45., 43., 41., 23., 41., 39., 46., 37., 42., 24.,
              40., 46., 39., 37., 49., 47., 37., 28., 36., 34., 34.,
              43., 41., 48., 42., 44., 45., 43., 45., 36., 33., 48.,
              39., 54., 35., 30., 42., 34., 51., 39., 50., 36., 34.,
              50., 37., 48., 40., 39., 21., 41., 40., 42., 44., 48.,
              42., 42., 39., 43., 57., 37., 36., 38., 34., 50., 40.,
              30., 29., 38., 44., 43., 47., 35., 37., 44., 44., 32.,
              36., 41., 44., 28., 43., 36., 37., 34., 37., 0., 0., 0.,
              0., 54., 34., 45., 0., 0., 0., 0., 0., 44., 32., 28.,
              34., 42., 35., 40., 45., 42., 47., 39., 35., 52., 34.,
              55., 38., 1.])

data = data[140:510]
ipos = np.where(data>0)[0]
data = np.log(data[ipos])
x    = np.arange(140,510,1.0)[ipos]

mod  = lmfit.models.GaussianModel()  + lmfit.models.ConstantModel()

pars = mod.make_params(c=data.mean(),
                       center=x.mean(),
                       sigma=x.std(),
                       amplitude=x.std()*data.ptp())

out = mod.fit(data, pars, x=x)

print(out.fit_report())

plt.plot(x, out.best_fit, 'r', label='Gaussian Model')
plt.plot(x, data)
plt.show()

This will print a report of

[Model]]
    (Model(gaussian) + Model(constant))
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 48
    # data points      = 333
    # variables        = 4
    chi-square         = 9.50683177
    reduced chi-square = 0.02889615
    Akaike info crit   = -1176.19189
    Bayesian info crit = -1160.95932
[[Variables]]
    amplitude:  58.7086539 +/- 5.27950257 (8.99%) (init = 153.376)
    center:     302.462920 +/- 2.79518589 (0.92%) (init = 322.6336)
    sigma:      48.1997878 +/- 3.66967747 (7.61%) (init = 106.8759)
    c:          3.64131210 +/- 0.01679060 (0.46%) (init = 3.796494)
    fwhm:       113.501824 +/- 8.64142989 (7.61%) == '2.3548200*sigma'
    height:     0.48592258 +/- 0.02641186 (5.44%) == '0.3989423*amplitude/max(1e-15, sigma)'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, c)     = -0.832
    C(amplitude, sigma) =  0.798
    C(sigma, c)         = -0.664

and a plot of

Is there any way I can access the uncertainty of variables?. For example, by uncertainty of amplitude I mean 5.27950257 — Pranav45, May 19 '21 at 14:00
best-fit parameters and their uncertainties will be in `out.params`, a dictionary with parameter names as keys, and values being Parameter objects with attributes including `value`, `stderr`, `vary`, `min`, `max`, `correl`, etc. So use `out.params['amplitude'].value`, `out.params['amplitude'].stderr`, etc — M Newville, May 19 '21 at 14:18

Gaussian Fit using lmfit

1 Answers1