lmfit not curve fitting

[Brennan Rodgers]

Hello,

I am trying to use lmfit to do some basic curve fitting to a gaussian but none of it works. I follow through the examples on the website and adapt the code straight from there but instead of being a nice bell curve, it's a jagged line that isn't much of a fit at all. I'm not sure what I'm doing wrong. I've tried writing my own function, I've tried using the built in Gaussian, the Voigt, etc, and none of it works. I've got some astronomy data that I want to curve fit that should be fairly straight forward and some exoplanet data that needs to be fit just like the Gaussians on the webpage. Any idea what I'm doing wrong?

[Brennan Rodgers]

Here is a sample of the form of exoplanet data I'm hoping to curve fit too. It will be a straight line and then a dip so kind of like an upside down gaussian. Curve fit and lmfit havn't been able to handle this at all so Ive been using a polyfit which is obviously not good enough. Any ideas on how to get a nice fit to this kind of data please let me know.

[Matt Newville]

Hi Brennan,

On Sun, Jul 15, 2018 at 3:39 PM Brennan Rodgers <brennanr...@gmail.com> wrote:

Hello,

I am trying to use lmfit to do some basic curve fitting to a gaussian but none of it works. I follow through the examples on the website and adapt the code straight from there but instead of being a nice bell curve, it's a jagged line that isn't much of a fit at all. I'm not sure what I'm doing wrong. I've tried writing my own function, I've tried using the built in Gaussian, the Voigt, etc, and none of it works. I've got some astronomy data that I want to curve fit that should be fairly straight forward and some exoplanet data that needs to be fit just like the Gaussians on the webpage. Any idea what I'm doing wrong?

​Uh, well, what have you tried?    That's kind of a lot of "nothing I tried work" without any indication of a) what it is you actually tried, and b) what you mean by "works".   With what you have posted, I can't imagine how anyone would ever be able to help you.  

Please attach ​a complete, minimal example of what you have tried.  Attach actual code, *NOT* an image of code.  Real Text.  

The plot you attach here doesn't look very Gaussian or Peak-like to me.   Does it look that way to you?  Maybe I'm not understanding what you mean. 

--Matt 

[Matt Newville]

Brennan, 

On Sun, Jul 15, 2018 at 3:44 PM Brennan Rodgers <brennanr...@gmail.com> wrote:

Here is a sample of the form of exoplanet data I'm hoping to curve fit too. It will be a straight line and then a dip so kind of like an upside down gaussian. Curve fit and lmfit havn't been able to handle this at all so Ive been using a polyfit which is obviously not good enough. Any ideas on how to get a nice fit to this kind of data please let me know.

​Again: what have you tried?  Post code and the results you get.   The plot you post here looks nothing like your data from 10 minuutes ago.  This does look slightly more peak-like.  You could certainly build a model that was a line or low order polynomial plus a Gaussian with a negative amplitude....  it sort of sounds like that is what you think you want.  

--Matt

​

[Brennan Rodgers]

Hello,

I was hoping you could gather what I was trying to do from the code I posted and the console output. I gave pictures of what I tried, that is an indication. Here is the code copied and pasted instead of a picture.Please forgive my ignorance, I am very new to this. 

import numpy as np

from numpy import exp,pi,sqrt

from scipy.optimize import curve_fit

import matplotlib.pyplot as plt

import random

import csv

from lmfit import Model

from lmfit.models import VoigtModel

#inputting the raw data as a CSV file

with open(input("Please enter the name of your CSV file "), 'r') as f:

    reader = csv.reader(f, delimiter=',')

    headers = next(reader)

    data = [row for row in reader]

    x_given = np.array([float(row[0]) for row in data])

    y_given = np.array([float(row[1]) for row in data])

#accounting for the data being in Julian Days

x_p = np.linspace(2458238,2458300,1054)

#this is the gaussian from the lmfit website I was trying but was unsuccessful with

#def gaussian(x,amplitude,center,width):

      #return (amplitude / (sqrt(2*pi) * width)) * exp(-(x-center)**2 / (width**2))

gmodel = VoigtModel()

pars = gmodel.guess(y_given,x=x_given)

result = gmodel.fit(y_given,pars,x=x_given)

print(result.fit_report())

plt.title('Nova Persei 2018')

plt.xlabel(headers[0])

plt.ylabel(headers[1])

plt.plot(x_given,y_given, 'D', markersize=1.8, color='green')

plt.plot(x_p,result.best_fit,'r-', markersize=100)

plt.legend(['Data','fit'], loc='best')

plt.ylim(5,13)

plt.xlim(2458238,2458307)

plt.gca().invert_yaxis()

So you can see here I'm trying to use the Voigt model. I've tried most of the built in models that are listed on the lmfit website/ "Working" in the context that I'm looking for would be like the example pictures I've attached here that are outputs from the lmfit website. The first example would be what I'm looking to do with the nova persei data and then the 2nd is what I'm looking for with the exoplanet stuff. Sorry for the confusion and being vague. Please forgive my ignorance as this is a new world for me. 

 Many thanks,

     -Brennan

[Brennan Rodgers]

When I change x_p = np.linspace(2458238,2458300,1054) to x_p = np.linspace(0,10,11) then my mock exoplanet test file looks like this (attached)

[Matt Newville]

On Sun, Jul 15, 2018 at 4:34 PM Brennan Rodgers <brennanr...@gmail.com> wrote:

Hello,

I was hoping you could gather what I was trying to do from the code I posted and the console output. I gave pictures of what I tried, that is an indication. Here is the code copied and pasted instead of a picture.Please forgive my ignorance, I am very new to this. 

​For everyone here: We seem to be having a lot of people being very confused about how to ask questions on this forum.  I've included instructions in the Welcome page to the Google group.

 

So you can see here I'm trying to use the Voigt model.

​Then why does your code include a copy-and-pasted definition for gaussian?

​

I've tried most of the built in models that are listed on the lmfit website/​​

"Working" in the context that I'm looking for would be like the example pictures I've attached here that are outputs from the lmfit website.

The first example would be what I'm looking to do with the nova persei data and then the 2nd is what I'm looking for with the exoplanet stuff. Sorry for the confusion and being vague. Please forgive my ignorance as this is a new world for me. 

​

You don't say what the starting values for your parameters are.   You don't include the fit report.  You keep posting plots of different data.

> When I change x_p = np.linspace(2458238,2458300,1054) to x_p = np.linspace(0,10,11) then my mock exoplanet test file looks like this (attached)

That is changing the data, of course and will certainly change the results...

I would suggest that you begin by plotting your data the way you intend to actually view it.  Then make a few test Gaussian or Voigt functions with parameters for center, amplitue, and sigma around whay you except,  and plot them with your data to see if these look close to right.   Once you get a peak that is remotely close, the fit should then be able to refine your starting values.  

--Matt 

[Brennan Rodgers]

I'm not sure why you're so confused. There are 2 different graphs that I am trying to fit to, I've said that quite clearly and provided picture examples. I'm not sure how to make it any clearer for you at this point. The definition for the gaussian was commented out, if you actually read the comment you will understand why I included it. I feel like you're not very interested in being helpful.

[Brennan Rodgers]

I have been testing with the Voigt function and have included the console output from those. If you actually look at those pictures you can see it is not remotely close and so I am unsure of where to go from there.

[Brennan Rodgers]

The change in linespace is to account for a new data set, that's why I said it was for the exoplanet data. It is a DIFFERENT set of data, which is why I had to change the values. The point is to show you a set of values and then the voigt fit which is not fitting very well at all. Please read the message carefully.

[Brennan Rodgers]

I am plotting it the way I intend to view it. You are not understanding me.

[Brennan Rodgers]

Let's try this again.

import numpy as np

from numpy import exp,pi,sqrt

from scipy.optimize import curve_fit

import matplotlib.pyplot as plt

import random

import csv

from lmfit import Model

from lmfit.models import VoigtModel

with open(input("Please enter the name of your CSV file "), 'r') as f:

    reader = csv.reader(f, delimiter=',')

    headers = next(reader)

    data = [row for row in reader]

    x_given = np.array([float(row[0]) for row in data])

    y_given = np.array([float(row[1]) for row in data])

x_p = np.linspace(0,10,11)

gmodel = VoigtModel()

pars = gmodel.guess(y_given,x=x_given)

result = gmodel.fit(y_given,pars,x=x_given)

print(result.fit_report())

plt.title('Exoplanet')

plt.xlabel(headers[0])

plt.ylabel(headers[1])

plt.plot(x_given,y_given, 'D', markersize=1.8, color='green')

plt.plot(x_p,result.best_fit,'r-', markersize=100)

plt.legend(['Data','fit'], loc='best')

plt.ylim(5,13)

plt.gca().invert_yaxis()

I'm trying to fit to a simple test file. Here is the fit report: 

Please enter the name of your CSV file testcsv.csv

[[Model]]

    Model(voigt)

[[Fit Statistics]]

    # fitting method   = leastsq

    # function evals   = 21

    # data points      = 11

    # variables        = 3

    chi-square         = 5.97326314

    reduced chi-square = 0.74665789

    Akaike info crit   = -0.71662096

    Bayesian info crit = 0.47706486

[[Variables]]

    sigma:      3.11837373 +/- 0.34804360 (11.16%) (init = 0.65)

    center:     6.00000000 +/- 0.30294095 (5.05%) (init = 6)

    amplitude:  133.991093 +/- 11.4024767 (8.51%) (init = 18)

    gamma:      3.11837373 +/- 0.34804360 (11.16%) == 'sigma'

    fwhm:       11.2302305 +/- 1.25341290 (11.16%) == '3.6013100*sigma'

    height:     8.96787463 +/- 0.38487555 (4.29%) == 'amplitude*wofz((1j*gamma)/(sigma*sqrt(2))).real/(sigma*sqrt(2*pi))'

[[Correlations]] (unreported correlations are < 0.100)

    C(sigma, amplitude) =  0.910

Attached is a picture of the console output. Are there any suggestions that pertain to curve fitting and how to modify the code to fit better please?

[Matt Newville]

On Sun, Jul 15, 2018 at 5:40 PM Brennan Rodgers <brennanr...@gmail.com> wrote:

Let's try this again.

 

Thank you for including the fit report.  

​  It looks to me like the fit "worked" -- it ran to completion, altered the starting values and even estimated uncertainties in the parameters.   The automagically guessed starting values were even close to the right order of magnitude.   The fit might not be "good" compared to what you were expecting, but I consider that this "worked".

But, why are you plotting the data again `x_given` and the fit agaist `x_p`?  It sure looks like `x_given` starts at 1, not 0, while `x_p` starts at 0.

​Also, ​

There is certainly a large offset

​ to the data -- it does not go to 0 far from the peak center.  You would almost certainly want to include a constant or linear component to acount for that.

--Matt

[Brennan Rodgers]

Interesting! So how do I go from what I have now to something like in the picture that I have attached? Again apologies for my ignorance on the subject. You make a good point about the mismatch at starting x values, I believe that is just a mistake on my part. 

import numpy as np

from numpy import exp,pi,sqrt

from scipy.optimize import curve_fit

import matplotlib.pyplot as plt

import csv

from lmfit import Model

from lmfit.models import VoigtModel

with open(input("Please enter the name of your CSV file "), 'r') as f:

    reader = csv.reader(f, delimiter=',')

    headers = next(reader)

    data = [row for row in reader]

    x_given = np.array([float(row[0]) for row in data])

    y_given = np.array([float(row[1]) for row in data])

gmodel = VoigtModel()

pars = gmodel.guess(y_given,x=x_given)

result = gmodel.fit(y_given,pars,x=x_given)

print(result.fit_report())

plt.title('Exoplanet')

plt.xlabel(headers[0])

plt.ylabel(headers[1])

plt.plot(x_given,y_given, 'D', markersize=1.8, color='green')

plt.plot(x_given,result.best_fit,'r-', markersize=100)

plt.legend(['Data','fit'], loc='best')

plt.ylim(5,13)

plt.gca().invert_yaxis()

Fit report for this is now: 

Please enter the name of your CSV file testcsv.csv

[[Model]]

    Model(voigt)

[[Fit Statistics]]

    # fitting method   = leastsq

    # function evals   = 21

    # data points      = 11

    # variables        = 3

    chi-square         = 5.97326314

    reduced chi-square = 0.74665789

    Akaike info crit   = -0.71662096

    Bayesian info crit = 0.47706486

[[Variables]]

    sigma:      3.11837373 +/- 0.34804360 (11.16%) (init = 0.65)

    center:     6.00000000 +/- 0.30294095 (5.05%) (init = 6)

    amplitude:  133.991093 +/- 11.4024767 (8.51%) (init = 18)

    gamma:      3.11837373 +/- 0.34804360 (11.16%) == 'sigma'

    fwhm:       11.2302305 +/- 1.25341290 (11.16%) == '3.6013100*sigma'

    height:     8.96787463 +/- 0.38487555 (4.29%) == 'amplitude*wofz((1j*gamma)/(sigma*sqrt(2))).real/(sigma*sqrt(2*pi))'

[[Correlations]] (unreported correlations are < 0.100)

    C(sigma, amplitude) =  0.910

And the graph is attached. Looks better! How do I go from here to having a nice fit like what's shown on the lmfit website? 

Many thanks for your help. Again, apologies for my ignorance on the subject.

[Matt Newville]

On Sun, Jul 15, 2018 at 6:16 PM Brennan Rodgers <brennanr...@gmail.com> wrote:

Interesting! So how do I go from what I have now to something like in the picture that I have attached?

Um, I think you might really mean "how do I make the fit better",  not "how do I make the fit look better".  There's a difference ;).

I would guess that you want to add a linear or even a constant background, say following the examples at 

    https://lmfit.github.io/lmfit-py/model.html#composite-models-adding-or-multiplying-models

​

[Brennan Rodgers]

Interesting! Thanks for the link! I'll check it out. 

[Brennan Rodgers]

import numpy as np

from numpy import exp,pi,sqrt

from scipy.optimize import curve_fit

import matplotlib.pyplot as plt

import csv

from lmfit import Model

from lmfit.models import VoigtModel

with open(input("Please enter the name of your CSV file "), 'r') as f:

    reader = csv.reader(f, delimiter=',')

    headers = next(reader)

    data = [row for row in reader]

    x_given = np.array([float(row[0]) for row in data])

    y_given = np.array([float(row[1]) for row in data])

#A line

def line(x,slope,intercept):

    return slope*x + intercept

gmodel = VoigtModel() + Model(line)

pars = gmodel.make_params(amplitude=5, center=5, width=1, slope=0, intercept=1)

result = gmodel.fit(y_given,pars,x=x_given)

print(result.fit_report())

plt.title('Exoplanet')

plt.xlabel(headers[0])

plt.ylabel(headers[1])

plt.plot(x_given,y_given, 'D', markersize=1.8, color='green')

plt.plot(x_given,result.best_fit,'r-', markersize=1)

plt.legend(['Data','fit'], loc='best')

plt.ylim(5,13)

plt.gca().invert_yaxis()

I added a linear function to the model, it looks much better now thank you! It's still not quite perfect though, do you have any suggestions on how to get an even better fit?

[Matt Newville]

​The standard approach is to try several different models (say, Gaussian, Voigt, etc) and different backgrounds (say, linear, constant, quadratic), and maybe try different starting values and compare the goodness-of-fit statistics like reduced chi-square, AIC, and BIC to find which model and set of parameter values best explains the data.

--Matt