SineModel fit does not work in certain cases

[lneuhaus]

Description

Certain perfect sines are imperfectly fitted by the SineModel. This is a result of constraining the phase shift of the model to the range of 0 to 2pi. If this constraint is lifted, the fit works properly.

The ideal solution to this problem would be a post-fit-hook, that allows to run the optimization without constraint on the range of the phase shift, but nevertheless returns a phase shift in a well-defined range. Picking a smarter optimization algorithm would also be a nice solution, though probably a less powerful one.

A Minimal, Complete, and Verifiable example

import numpy as np
from lmfit.models import SineModel, ConstantModel

for unconstrained in [False, True]:
    print(f"CONSTAINT LIFTED: {unconstrained}")
    for alpha in np.linspace(-1, 1, 10):
        x = np.linspace(0, 25, 20)
        data = np.sin(alpha + 0.3 * x)
        sine_model = SineModel()
        if unconstrained:
            sine_model.set_param_hint('shift', min=-5, max=10)
        s_params = sine_model.guess(data, x=x)
        background = ConstantModel()
        b_params = background.guess(data)
        model = sine_model + background
        pars = s_params + b_params
        result = model.fit(data, pars, x=x, method="least_squares")
        params = result.best_values
        shift = params["shift"]
        frequency = params["frequency"]
        diff = (alpha - shift) % (2.0 * np.pi)
        if np.isclose(diff, 0) or np.isclose(diff, 2*np.pi):
            print(f"alpha {alpha}: ok")
        else:
            print (f"alpha {alpha}, alpha%2pi {alpha%(2*np.pi):.3f}, shift%2pi {shift%(2*np.pi):.3f}")

Error message:

This is the output I get:

CONSTAINT LIFTED: False
alpha -1.0: ok
alpha -0.7777777777777778, alpha%2pi 5.505, shift%2pi 0.000
alpha -0.5555555555555556, alpha%2pi 5.728, shift%2pi 0.000
alpha -0.33333333333333337, alpha%2pi 5.950, shift%2pi 0.000
alpha -0.11111111111111116, alpha%2pi 6.172, shift%2pi 0.000
alpha 0.11111111111111116: ok
alpha 0.33333333333333326: ok
alpha 0.5555555555555554: ok
alpha 0.7777777777777777: ok
alpha 1.0: ok
CONSTAINT LIFTED: True
alpha -1.0: ok
alpha -0.7777777777777778: ok
alpha -0.5555555555555556: ok
alpha -0.33333333333333337: ok
alpha -0.11111111111111116: ok
alpha 0.11111111111111116: ok
alpha 0.33333333333333326: ok
alpha 0.5555555555555554: ok
alpha 0.7777777777777777: ok
alpha 1.0: ok

Version information

Python: 3.6.8 |Anaconda, Inc.| (default, Dec 30 2018, 01:22:34) 
[GCC 7.3.0]

lmfit: 1.0.1+88.gc3d0772, scipy: 1.4.1, numpy: 1.18.1, asteval: 0.9.17, uncertainties: 3.1

[newville]

@lneuhaus I'm afraid I never understood the actual issue here. Can you clarify?

[lneuhaus]

Sure, sorry if it wasn't clear.

Above example code generates data points that lie on a perfect (=no noise) sine (np.sin(alpha + 0.3 * x)) and sweeps the phase shift alpha from -1 to +1 radian. It then uses the SineModel somewhat naively to fit a sine to this data. Certain phase shifts slightly below 2*pi don't give a good fit. I thought I had such cases handled by unit tests but it seems I did not.

A fix of this issue should likely add above example code as a unit test, which would fail with the current code base. My question is how to improve the SineModel to work in these cases.

1. One solution is to remove the constraint the model defines for the phase shift of the sine (range from 0 to 2 pi), as the SineModel starts to give reasonable fits in this case. This is the solution I implemented in my application code. If you agree I could make a PR with this solution as it would already be an improvement.

2. A better solution would add a post-fit hook that would map the fitted phase shift back to a well-defined interval (0 to 2 pi), saving the user code one or two lines. We had discussed this in #676 and you were not against it, but felt it wasn't required as we could get away by using constraints. I think this example demonstrates that we do not get away so easily, and that there is a difference between fitting with a truly periodic phase shift and one that is constrained to the interval 0 to 2 pi. If you point me to where such a post-fit hook could be added I could give it a try. But I share your concern that this could become a lot of work due to handling of errors etc.

3. You probably know the details of the default optimizer used in lmfit much better than I, so maybe you have an intuition as to why it is getting stuck in a local optimum at one end of the constraint interval 0 to 2 pi (by going slightly negative it would reduce the error but it just stays at a phase shift of 0 due to the constraint - it seems to be too dumb to try jumping to 2 pi and continuing from there). If there is another optimization algorithm to be used for periodic constraints that should behave better, maybe that would be the better solution.

[newville]

On Fri, Nov 20, 2020 at 12:32 AM lneuhaus notifications@github.com wrote:

Sure, sorry if it wasn't clear.

Above example code generates data points that lie on a perfect (=no noise) sine (np.sin(alpha + 0.3 x)) and sweeps the phase shift alpha from -1 to +1 radian. It then uses the SineModel somewhat naively to fit a sine to this data. Certain phase shifts slightly below 2pi don't give a good fit. I thought I had such cases handled by unit tests but it seems I did not.

Sorry for being dense here, but I am actually not seeing this. Normally an Issue should show minimal code that shows the problem and give complete output. And, well explain what the problem is. Where is the "I expected A and got B" here? Maybe I'm just not parsing your script and output.

If the issue is that an initial guess close to a parameter's bound is giving poor results, show that.

A fix of this issue should likely add above example code as a unit test,

which would fail with the current code base. My question is how to improve the SineModel to work in these cases.

1.

One solution is to remove the constraint the model defines for the phase shift of the sine (range from 0 to 2 pi), as the SineModel starts to give reasonable fits in this case. This is the solution I implemented in my application code. If you agree I could make a PR with this solution as it would already be an improvement.

That would seem OK. Or at least it might be reasonable to better inform a user how to lift the bounds set by parameter hints that are meant to be "model specific".

1.

A better solution would add a post-fit hook that would map the fitted phase shift back to a well-defined interval (0 to 2 pi), saving the user code one or two lines. We had discussed this in #676 https://github.com/lmfit/lmfit-py/pull/676 and you were not against it, but felt it wasn't required as we could get away by using constraints. I think this example demonstrates that we do not get away so easily, and that there is a difference between fitting with a truly periodic phase shift and one that is constrained to the interval 0 to 2 pi. If you point me to where such a post-fit hook could be added I could give it a try. But I share your concern that this could become a lot of work due to handling of errors etc.

Well, a post-fit hook would be significant work. We don't have that for any other code. It is worth considering, but would not be a quick fix for any model. I'm not sure it is "better", but it might be worth thinking about. Unfortunately, that won't come from me anytime soon.

1. 2.

You probably know the details of the default optimizer used in lmfit much better than I, so maybe you have an intuition as to why it is getting stuck in a local optimum at one end of the constraint interval 0 to 2 pi (by going slightly negative it would reduce the error but it just stays at a phase shift of 0 due to the constraint - it seems to be too dumb to try jumping to 2 pi and continuing from there). If there is another optimization algorithm to be used for periodic constraints that should behave better, maybe that would be the better solution.

None of the methods used for least-squares will try a large "jump" in a parameter value. They are iterative refinement methods. Some of the slow, pessimistic methods like differential evolution, basin hopping, and ampgo may make large jumps, as they (mostly) ignore the derivative of cost vs parameter_value.

If a parameter has bounds of [0, 1], the best value is really 0.25, and you guess 0.999, the solution may not be found. The short version of advice is: don't guess values very close to the parameter bounds.

[lneuhaus]

Ok, I will work out a PR for 1 and propose to add a new method best_values_postprocessed to access post-processed fit results. That wouldn't require much work to implement and allow to test the usefulness of the concept.

[newville]

@lneuhaus sorry for continuing to be dense.... a PR to do what?

[lneuhaus]

I was referring to:

application code. If you agree I could make a PR with this solution as it would already be an improvement.

That would seem OK.

[newville]

@lneuhaus We have gone back and forth a few times here and I still do not understand the exact problem you are having or the proposed fix. It seems like this really ought to be not that complicated. Please describe and give a simple, minimal example that clearly shows the problem you are having. Like, honestly you posted a bunch of code above and everything says "ok". So: what the heck is the problem?

[lneuhaus]

I will start my PR with a unit test. You will understand from the test failure.

[newville]

@lneuhaus Since you are unable to provide a case that shows the problem, I tried to do so. I failed to find a problem. All cases work fine for me.

import numpy as np
from lmfit.models import SineModel

x = np.linspace(0, 25, 101)
smodel = SineModel()

for xshift in (0.5, 1, 1.5, 2.0, 2.5):
    params = smodel.make_params(amplitude=1.25, shift=xshift, frequency=0.25)
    print("## Guess shift = {:.4f}".format(xshift))
    print("ActualShift  FittedShift    Reduced Chi-square")
    for dshift in np.arange(0.1, 3.1, 0.2):
        data = np.sin(dshift + 0.3 * x) + np.random.normal(scale=0.05, size=len(x))
        result = smodel.fit(data, params, x=x)
        print("{:.4f}, {:.4f}, {:.4f}".format(dshift,
                                              result.params['shift'].value,
                                              result.redchi))

which gives

## Guess shift = 0.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1273, 0.0027
0.3000, 0.2777, 0.0023
0.5000, 0.5056, 0.0020
0.7000, 0.7187, 0.0028
0.9000, 0.8842, 0.0025
1.1000, 1.0848, 0.0027
1.3000, 1.3163, 0.0029
1.5000, 1.5101, 0.0024
1.7000, 1.7077, 0.0030
1.9000, 1.8714, 0.0027
2.1000, 2.0881, 0.0023
2.3000, 2.3116, 0.0026
2.5000, 2.4718, 0.0025
2.7000, 2.7099, 0.0018
2.9000, 2.9052, 0.0022
## Guess shift = 1.0000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.0954, 0.0022
0.3000, 0.2663, 0.0020
0.5000, 0.4741, 0.0033
0.7000, 0.6917, 0.0026
0.9000, 0.9056, 0.0030
1.1000, 1.0808, 0.0023
1.3000, 1.3114, 0.0024
1.5000, 1.4948, 0.0024
1.7000, 1.6835, 0.0024
1.9000, 1.9056, 0.0023
2.1000, 2.1148, 0.0024
2.3000, 2.3005, 0.0023
2.5000, 2.5188, 0.0028
2.7000, 2.7180, 0.0024
2.9000, 2.0802, 0.4682
## Guess shift = 1.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1093, 0.0024
0.3000, 0.2969, 0.0022
0.5000, 0.4967, 0.0024
0.7000, 0.7072, 0.0022
0.9000, 0.8918, 0.0024
1.1000, 1.0754, 0.0029
1.3000, 1.3000, 0.0030
1.5000, 1.4814, 0.0027
1.7000, 1.6947, 0.0031
1.9000, 1.8852, 0.0025
2.1000, 2.1085, 0.0025
2.3000, 2.2940, 0.0030
2.5000, 2.5043, 0.0024
2.7000, 2.6920, 0.0028
2.9000, 2.9221, 0.0018
## Guess shift = 2.0000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1130, 0.0025
0.3000, 0.3042, 0.0029
0.5000, 0.5262, 0.0022
0.7000, 0.6868, 0.0025
0.9000, 0.9357, 0.0023
1.1000, 1.1024, 0.0027
1.3000, 1.2868, 0.0026
1.5000, 1.5179, 0.0023
1.7000, 1.6964, 0.0024
1.9000, 1.8968, 0.0028
2.1000, 2.1155, 0.0022
2.3000, 2.2976, 0.0015
2.5000, 2.5066, 0.0024
2.7000, 2.6892, 0.0023
2.9000, 2.8727, 0.0024
## Guess shift = 2.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1181, 0.0026
0.3000, 0.2979, 0.0032
0.5000, 0.4984, 0.0027
0.7000, 0.7044, 0.0029
0.9000, 0.9040, 0.0027
1.1000, 1.1043, 0.0024
1.3000, 1.3075, 0.0027
1.5000, 1.5060, 0.0024
1.7000, 1.6954, 0.0025
1.9000, 1.8969, 0.0022
2.1000, 2.0941, 0.0021
2.3000, 2.3169, 0.0030
2.5000, 2.4894, 0.0022
2.7000, 2.6971, 0.0025
2.9000, 2.8990, 0.0022

All fits give good results. If you see a problem, then explain what that is.

[lneuhaus]

Very easy: replace 3.1 by 6.2 in your code and you get cases where things go off:

## Guess shift = 0.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1244, 0.0022
0.3000, 0.2894, 0.0018
0.5000, 0.4924, 0.0029
0.7000, 0.6919, 0.0018
0.9000, 0.9034, 0.0027
1.1000, 1.1076, 0.0031
1.3000, 1.2978, 0.0027
1.5000, 1.5069, 0.0022
1.7000, 1.6898, 0.0027
1.9000, 1.9041, 0.0027
2.1000, 2.1065, 0.0025
2.3000, 2.3262, 0.0022
2.5000, 2.5023, 0.0030
2.7000, 2.7010, 0.0023
2.9000, 2.8951, 0.0022
3.1000, 1.0056, 0.4782
3.3000, 0.6887, 0.5101
3.5000, 0.3208, 0.5501
3.7000, 0.0000, 0.5060
3.9000, 0.1924, 0.5711
4.1000, 0.0000, 0.4417
4.3000, 0.0000, 0.4173
4.5000, 0.0000, 0.3577
4.7000, 0.0000, 0.2812
4.9000, 0.0000, 0.2343
5.1000, 0.0000, 0.1830
5.3000, 0.0000, 0.1402
5.5000, 0.0000, 0.0897
5.7000, 0.0000, 0.0526
5.9000, 0.0000, 0.0249
6.1000, 0.0000, 0.0067
## Guess shift = 1.0000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.0882, 0.0018
0.3000, 0.2954, 0.0029
0.5000, 0.5050, 0.0032
0.7000, 0.7180, 0.0027
0.9000, 0.9018, 0.0024
1.1000, 1.1071, 0.0021
1.3000, 1.3151, 0.0026
1.5000, 1.5036, 0.0020
1.7000, 1.6937, 0.0018
1.9000, 1.8993, 0.0024
2.1000, 2.0684, 0.0030
2.3000, 2.2902, 0.0018
2.5000, 2.5180, 0.0022
2.7000, 2.1128, 0.4567
2.9000, 2.1044, 0.4785
3.1000, 1.9371, 0.4819
3.3000, 1.7395, 0.5121
3.5000, 1.5263, 0.5347
3.7000, 1.2908, 0.5697
3.9000, 1.0686, 0.5641
4.1000, 0.7723, 0.5954
4.3000, 0.6915, 0.5710
4.5000, 0.5774, 0.5556
4.7000, 0.0000, 0.2916
4.9000, 0.0000, 0.2294
5.1000, 0.0000, 0.1864
5.3000, 0.0000, 0.1323
5.5000, 0.0000, 0.0943
5.7000, 0.0000, 0.0567
5.9000, 0.0000, 0.0259
6.1000, 0.0000, 0.0063
## Guess shift = 1.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.0997, 0.0023
0.3000, 0.2762, 0.0029
0.5000, 0.5084, 0.0022
0.7000, 0.7030, 0.0024
0.9000, 0.8637, 0.0026
1.1000, 1.0733, 0.0020
1.3000, 1.3102, 0.0022
1.5000, 1.5101, 0.0022
1.7000, 1.7008, 0.0025
1.9000, 1.8947, 0.0020
2.1000, 2.0937, 0.0027
2.3000, 2.3173, 0.0024
2.5000, 2.4986, 0.0024
2.7000, 2.7140, 0.0025
2.9000, 2.9192, 0.0023
3.1000, 2.3894, 0.4900
3.3000, 2.4014, 0.5311
3.5000, 2.4096, 0.5384
3.7000, 2.2640, 0.5599
3.9000, 2.0272, 0.5756
4.1000, 1.8777, 0.5937
4.3000, 1.7223, 0.5925
4.5000, 1.5353, 0.5519
4.7000, 1.3713, 0.5310
4.9000, 1.2442, 0.5181
5.1000, 0.0000, 0.1869
5.3000, 0.0000, 0.1426
5.5000, 0.0000, 0.0903
5.7000, 0.0000, 0.0573
5.9000, 0.0000, 0.0239
6.1000, 0.0000, 0.0069
## Guess shift = 2.0000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1139, 0.0024
0.3000, 0.2803, 0.0022
0.5000, 0.4878, 0.0028
0.7000, 0.6751, 0.0022
0.9000, 0.9148, 0.0026
1.1000, 1.1181, 0.0028
1.3000, 1.2933, 0.0034
1.5000, 1.5114, 0.0019
1.7000, 1.6575, 0.0026
1.9000, 1.8934, 0.0029
2.1000, 2.0977, 0.0029
2.3000, 2.2947, 0.0024
2.5000, 2.5090, 0.0025
2.7000, 2.6650, 0.0029
2.9000, 2.9068, 0.0017
3.1000, 3.1206, 0.0022
3.3000, 3.2920, 0.0022
3.5000, 3.4708, 0.0023
3.7000, 2.8389, 0.5606
3.9000, 2.7851, 0.5756
4.1000, 2.7327, 0.5715
4.3000, 2.6460, 0.5674
4.5000, 2.4820, 0.5503
4.7000, 2.3461, 0.5537
4.9000, 2.1542, 0.5185
5.1000, 2.0273, 0.4957
5.3000, 1.8914, 0.4722
5.5000, 0.2878, 0.4552
5.7000, 0.0000, 0.0562
5.9000, 0.0000, 0.0267
6.1000, 0.0000, 0.0074
## Guess shift = 2.5000
ActualShift  FittedShift    Reduced Chi-square
0.1000, 0.1010, 0.0028
0.3000, 0.3033, 0.0028
0.5000, 0.4782, 0.0028
0.7000, 0.7152, 0.0022
0.9000, 0.9318, 0.0025
1.1000, 1.1127, 0.0022
1.3000, 1.3015, 0.0021
1.5000, 1.4927, 0.0031
1.7000, 1.6822, 0.0029
1.9000, 1.8878, 0.0024
2.1000, 2.1077, 0.0022
2.3000, 2.3217, 0.0029
2.5000, 2.5041, 0.0027
2.7000, 2.6883, 0.0029
2.9000, 2.8834, 0.0023
3.1000, 3.1335, 0.0030
3.3000, 3.3080, 0.0024
3.5000, 3.4968, 0.0029
3.7000, 3.6813, 0.0031
3.9000, 3.9022, 0.0023
4.1000, 4.0919, 0.0022
4.3000, 4.3026, 0.0025
4.5000, 3.4006, 0.5667
4.7000, 4.7069, 0.0022
4.9000, 3.0933, 0.5192
5.1000, 2.9417, 0.4934
5.3000, 2.7379, 0.4712
5.5000, 2.5793, 0.4634
5.7000, 2.4445, 0.4482
5.9000, 2.2939, 0.4632
6.1000, 0.0000, 0.0069

[lneuhaus]

Though to be fair, I think this example is too minimalistic: it does not use the SineModel's guess function.

[newville]

@lneuhaus I don't think we can really consider a fit failing because of poor initial guesses to be an error. One can almost always make a fit fail, if that is the goal. So basically, the issue title "does not work in certain cases" strikes me as sort of normal.

I cannot tell if you are saying that the issue is that SineModel.guess() is not always good enough? I definitely did not pick that up from the conversation so far. I'm pretty sure we can't really guarantee that automated guessing will always work well in all cases either. It does work in very many cases.

[lneuhaus]

I cannot tell if you are saying that the issue is that SineModel.guess() is not always good enough?

Yes, exactly my point. And that the issue is not the fact that the guesses are way off, but that the best guess sometimes lies close to a constraint.

I'm pretty sure we can't really guarantee that automated guessing will always work well in all cases either. It does work in very many cases.

While this is certainly true for many complicated or noisy cases, I would expect from a fitting library to properly fit noiseless perfect sines without issues out of the box, including appropriate guessing.

I claim that the current implementation in lmfit does not guarantee this, but that there is a way to guarantee this (by not constraining the phase shift variable).

I will demonstrate this in a unit test, no need to loose time in constructing a parallel example here.

[lneuhaus]

See #682

[newville]

@lneuhaus

While this is certainly true for many complicated or noisy cases, I would expect from a fitting library to properly fit noiseless perfect sines without issues out of the box, including appropriate guessing.

I think your expectation needs recalibration. We simply cannot guarantee that fits will succeed independently from the initial parameter values. Really, I do mean "can not".

I claim that the current implementation in lmfit does not guarantee this, but that there is a way to guarantee this (by not constraining the phase shift variable).

You are right that the current implementation does not guarantee this.

It may be that automated guesses could be better. Automated guesses will not ever be guaranteed to work.

For the specific case of the phase-shift for a sine wave, if the constraints are 0 and 2*pi, guessing pi is probably going to work for a great many cases.

I will demonstrate this in a unit test, no need to loose time in constructing a parallel example here.

We have already spent too much time on this not-exactly-a-problem situation.

[Tillsten]

Completely agreeing with @newville. Still there is an easy enchanment to the guess function. Since it already assumes uniform spacing and uses an FFT to estimate the frequency, it should also use the resulting phase instead of the used iterative approach.