Incapable of getting lmfit to fit an inverse Lorentz (Cauchy) Distribution

[Alexandros Gloor]

Dear all,

I am attempting to get information about the center of a trough of a gaussian that is subtracting from a constant value. I have, as you see below, attempted to create a model where the Lorentzian is subtracting from the constant and tried to run .fit to get more accurate results, but all lmfit does is return the initialized values and a graph that does not make much sense (I believe it's plotting something tiny at 0,0). 

The error I am getting when I run this code is not a Python error, but an lmfit error. Along the printout it includes a statement which I believe to be relevant;

##  Warning: uncertainties could not be estimated:

The code I am running goes as follows:

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker 

matData = pd.read_csv('materialNoCav.txt',sep=" ")
matData.rename(columns={'kpar (rad/m)': 'kPar','freq (Hz)':'freq','omega (rad/s)':'omega','Power outflow, time average (W/m), Boundary Probe 1':'power'}, inplace=True)
matPowerVals = matData.pivot(index='freq',columns='kPar',values='power')
matPowerArray = ((np.asarray(matData['power'])).reshape(21,51))
matKList = uniqueList(matData['kPar']) # extracts the kParallel values for our 21 value parameter sweep
matKParDict = {}
for kval in matKList:
    df = matData[matData.kPar==kval]
    matKParDict[kval] = df

 

from lmfit import Parameters
from lmfit.models import ConstantModel, LorentzianModel 

matResF = []
for i in matKList:
    kParTemp = pd.DataFrame(matKParDict[i]).reset_index(drop=True) # temporary DF of our dictionary section
    x=kParTemp['freq']
    y=kParTemp['power']
    
    const = ConstantModel(prefix='c_')
    pars = const.make_params()
    pars['c_c'].set(value=1,min=0.99,max=1)
    lor = LorentzianModel(prefix='l_')
    pars.update(lor.make_params())
    pars['l_amplitude'].set(value=1,min=0,max=1)
    pars['l_center'].set(value=5.6e14,min=5.2e14,max=6.6e14)
    
    model = const - lor
    out = model.fit(y, pars, x=x)
    print(out.fit_report(min_correl=0.5))
    comps = out.eval_components(x=x)
    plt.plot(x, y, comps['c_'], label='Constant')
    plt.plot(x, y, comps['l_'], label='Lorentzian')
    plt.legend(loc='best')
    plt.show()

I have attached the data file I am using. While it is not necessary, I am doing a sweep across all 21 datasets. If this issue can be resolved for i=0, it would be equally useful. 

The green attached plot is my current result, the orange-blue one is a pure lorentzian with a completely useless and wrong answer (since it's only running Lorentz and nothing else) that I am posting to show the raw datapoints. 

Please let me know if there's a better way to fit an inverse function. 

[mpm...@gmail.com]

Dear Alexandros,

Posting the full output of your script would be helpful for debugging. 

I suspect that your parameters might not have been passed to model.fit. The call signature is

model.fit(data, params=None).

Another possible problem is that you've initialized parameter(s) at its limit, such that lmfit cannot evaluate the function +- w.r.t. its starting value.

pars['c_c'].set(value=1,min=0.99,max=1) and similar to l_amplitude

Best,

Mark

[Matt Newville]

Alexandros, 

On Thu, Aug 13, 2020 at 8:04 PM Alexandros Gloor <alexand...@gmail.com> wrote:

Dear all,

I am attempting to get information about the center of a trough of a gaussian that is subtracting from a constant value. I have, as you see below, attempted to create a model where the Lorentzian is subtracting from the constant and tried to run .fit to get more accurate results, but all lmfit does is return the initialized values and a graph that does not make much sense (I believe it's plotting something tiny at 0,0). 

The error I am getting when I run this code is not a Python error, but an lmfit error. Along the printout it includes a statement which I believe to be relevant;

##  Warning: uncertainties could not be estimated:

Nope, that is not the error your code gives.  The error your code gives is that `pd` is not defined.  Even assuming your code should start with

    import pandas as pd

    import numpy as np

the error is that `uniqueList` is not defined.  

Really: post a complete, minimal example that you actually rand and the error you actually get.  A minimal example will certainly not loop over twenty datasets.  It is alarming how often we say this.  

Like Mark suggests, certainly doing 

  pars['c_c'].set(value=1,min=0.99,max=1)

is a really bad idea.  Perhaps you meant 

   pars['c_c'].set(value=1, vary=False)

?

 

Do you want to add a constant and a gaussian (or lorentzian -- you mention both) with a negative amplitude?

If so, you're definitely going through a lot of confusing steps to do that.  

The code I am running goes as follows:

Please let me know if there's a better way to fit an inverse function. 

--
You received this message because you are subscribed to the Google Groups "lmfit-py" group.
To unsubscribe from this group and stop receiving emails from it, send an email to lmfit-py+u...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/a69fe708-22b7-4e87-bce6-3dc67adedb22o%40googlegroups.com.

--Matt 

[Alexandros Gloor]

Matt, the first two issues you listed were probably an error due to my copy-pasting parts of my code into the mailing list. I wanted to avoid sending a whole jupyter file, but must have missed some imports/definitions.

That being said, it didn't cross my mind to make the amplitude negative, I've been subtracting the lorentzian with a positive amplitude which in the end has the same effect. That seems a simple enough solution that it might just resolve my problem. I'll be checking back on the full mailing list and reporting with any further information requested then.

Best,

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CA%2B7ESbohiZs9%3D0zjT4YSpiZ%3DeoFZbUgSNoJ%2BTs%2BFt7%2BkQJAO0w%40mail.gmail.com.

[Alexandros Gloor]

Dear Mark, 

This code is running fine for me (with the same dataframe as posted before). I missed one function in my first email. As per Matt's recommendation I have fixed the values I needed constant. When you plot this graph you will see the green shape I posted previously. The data, and the fit, should only be spanning from ~5.2e14 to 6.6e14 hence my limiting min/max values for the frequency parameter. I am not sure why model.fit would not be reading the parameters I have provided as both the model output and pars itself is giving me the initial values I entered, without fitting them to anything. Hopefully the attached file is more helpful for you! 

Best,

Alexandros

[mpm...@gmail.com]

Dear Alexandros,

I can see the solution now, but I have to come back to the fact that the initial and follow up posts do not follow the instructions, even for issues we mentioned. You still didn't post the fit report for example. This makes it harder for people to help and probably impossible for other starting users to read and understand the solution. 

I hope this does not come across as grumpy, but please can you re-post a new thread that follows the instructions. I'll answer it there, which should leave a clear solution for others. 

Best,

Mark

If you are posting a question about lmfit not working for you as you would expect, please include in your post: 

• a minimal and complete example.  Please see https://stackoverflow.com/help/mcve for more details.
  
• the fit report from lmfit, if you have one.
• If an error occurs, include the full traceback.
• It's OK to post an image of data and fit, but if you are including code or output, this must be in plain text.  Do not post an image of text.