Do not force sigma to be positive in the step model

[Mark Dean]

Dear all,

This line in the fit function

```    out = (x - center)/max(tiny, sigma)```

https://github.com/mpmdean/lmfit-py/blob/776e14d73028548f783288c37468b520e459050e/lmfit/lineshapes.py#L429

forces `sigma` to be positive.

In my opinion, it would be better to force `sigma` to be nonzero, as negative sigma would allow fits high-to-low fits like this

x = np.linspace(0, 5, 100)
y = np.hstack([[0]*50, [1]*50])

fig, ax = plt.subplots()

model = lmfit.models.StepModel(form='erf')

result = model.fit(y, x=x)

which currently fails and fits sigma \approx zero. The fix is quite simple. It would involve changing

```    out = (x - center)/max(tiny, sigma)```

to

```

    if sigma >= 0:
        sigma = max(tiny, sigma)
    elif sigma < 0:
        sigma = min(tiny, sigma)
   
    out = (x - center)/sigma

```

Best,

Mark

[Matt Newville]

Hi Mark, 

(sorry for the delay in getting this posted).

On Sun, Oct 6, 2024 at 10:11 AM Mark Dean <mpm...@gmail.com> wrote:

Dear all,

This line in the fit function

```    out = (x - center)/max(tiny, sigma)```

https://github.com/mpmdean/lmfit-py/blob/776e14d73028548f783288c37468b520e459050e/lmfit/lineshapes.py#L429

forces `sigma` to be positive.

In my opinion, it would be better to force `sigma` to be nonzero, as negative sigma would allow fits high-to-low fits like this

Hm,  what does a negative value for sigma mean?

x = np.linspace(0, 5, 100)
y = np.hstack([[0]*50, [1]*50])

fig, ax = plt.subplots()

model = lmfit.models.StepModel(form='erf')

result = model.fit(y, x=x)

which currently fails and fits sigma \approx zero.

Well, you did not provide starting values for any of the parameters.  As it turns out, the default value for center is 0 and that for sigma is 1.  And your step is 1 'x' point wide.  So, the initial values are terrible, making this example of not finding a good solution not super-convincing of "a problem" to me.  

The fix is quite simple. It would involve changing

```    out = (x - center)/max(tiny, sigma)```

to

```

    if sigma >= 0:
        sigma = max(tiny, sigma)
    elif sigma < 0:
        sigma = min(tiny, sigma)
   
    out = (x - center)/sigma

Should sigma be allowed to change signs during a fit?     

It seems to me that this use of sigma<0 would mean "step down" instead of "step up".   That can be handled with an offset and a negative amplitude -- I'm not saying that is the only way to handle "step down", but "sigma<0" seems sort of odd for that too.

But, I'm not sure that is what you are asking, because your step is going up.

It's entirely possible that I'm missing the point you're trying to make...

--Matt

[Mark Dean]

Dear Matt,

Sorry for the error in the initial code snippet. Let me restate the comment with some further additions. The text below should contain answers to your questions.

import lmfit
import numpy as np
from scipy.special import erf
import matplotlib.pyplot as plt

fig, axs = plt.subplots(1, 2, figsize=(10, 5))

# Let's say I have a curve of this general form
x = np.linspace(-5, 5, 100)
y = (1 + erf(-x))/2
# Note that this a very common pump-probe trace shape in ultrafast science

# Reading the docs of lmfit, I thought I could fit this as follows

model = lmfit.models.StepModel(form='erf')

params = model.make_params()
params['sigma'].set(value=-1)

result = model.fit(y, x=x, params=params)
result.plot_fit(ax=axs[0], show_init=True)

# It took me quite a while to realize that lmfit was forcing +ve sigma 

# without this being documented.

# As you say, it is possible to fit like this

model = lmfit.models.StepModel(form='erf') + lmfit.models.ConstantModel()
params = model.make_params()
params['amplitude'].set(value=-1)
params['c'].set(value=1)

result = model.fit(y, x=x, params=params)
result.plot_fit(ax=axs[1], show_init=True)

# I think this is much more cumbersome and less natural

# as I have to both flip the step upside-down as 

# well as offset it with a constant, rather than

# just recognizing that the step has been flipped in x.  

# Regarding whether sigma should be allowed to change during a fit
# I think this should be up to the user.

# I would say that lmfit should allow all curve shapes that 

# make mathematical sense, as it cannot reasonably 

# hope to guess which shapes make physical sense

# . The default amplitude behavior allows sign switching

# even though most physical situations probably have 

# either positive or negative amplitudes. 

--
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/CA%2B7ESbpA0t78wwSJchifJvW%3Dp53tJg%3DaTkkz3hQ11E4CLWZUsQ%40mail.gmail.com.

--

--------------------------------------------

Mark P. M. Dean

Brookhaven National Laboratory

Condensed Matter Physics Bldg 734

Upton NY 11973

USA

Phone: 001 631 344 7847

Email: md...@bnl.gov

[Matt Newville]

Hi Mark, 

Thanks. I think I understand better that the idea of having a negative sigma does imply a step-down.   That seems OK,t though we'll have to document that, and I worry a bit about what happens if sigma is very small (like, on the order of tiny), but I'll assume that is pretty unlikely in practice. I would guess that most of the time, the user would know ahead of time if a step was required to be a step up or down and could impose that with a min/max boundary.

I'll plan on changing both StepModel and RectangleModel.

To view this discussion on the web visit https://groups.google.com/d/msgid/lmfit-py/CABMFFfpxut7j__BZCwAZK4J%2Bme%3DN___NntNWUWcvK16s1AUZJA%40mail.gmail.com.

--Matt