[SciPy-User] lmfit: confidence intervals issue

classic Classic list List threaded Threaded
9 messages Options
Reply | Threaded
Open this post in threaded view
|

[SciPy-User] lmfit: confidence intervals issue

Antonino Ingargiola
Hi,

I'm not sure this is the right list for this kind of issues.

I'm using lmfit to solve a non-linear least square problem to fit a "noisy" histogram to a model that is a single exponential decay convoluted with a "sharp" peak and a baseline noise floor.

The problem is that unless I set constrains for the 'offset' parameter the confidence interval calculation fails in a strange way in scipy brentq function.

You can see the full example notebook and the error message in this ipython notebook:


What can I do to fix the problem? If this is the normal behavior, a more informative error message would be helpful.

Thanks in advance for any help,
Antonio

PS: Let me congratulate with the lmfit authors as this is a really amazing piece of software that should be included in scipy IMHO!


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Sturla Molden-3
lmfit is not a part of SciPy, but as long as it's scientific use of Python
and no other appropriate list, I don't think anyone will mind. I cannot
answer your question, though.

Sturla

Antonino Ingargiola <[hidden email]> wrote:

> Hi,
>
> I'm not sure this is the right list for this kind of issues.
>
> I'm using lmfit to solve a non-linear least square problem to fit a "noisy"
> histogram to a model that is a single exponential decay convoluted with a
> "sharp" peak and a baseline noise floor.
>
> The problem is that unless I set constrains for the 'offset' parameter the
> confidence interval calculation fails in a strange way in scipy brentq
> function.
>
> You can see the full example notebook and the error message in this ipython
> notebook:
>
> <a
> href="http://nbviewer.ipython.org/urls/gist.githubusercontent.com/tritemio/901c2eb2a43775e81844/raw/755eede8f170cc33b96316601286ae534e901c49/lmfit%20issue">http://nbviewer.ipython.org/urls/gist.githubusercontent.com/tritemio/901c2eb2a43775e81844/raw/755eede8f170cc33b96316601286ae534e901c49/lmfit%20issue</a>
>
> What can I do to fix the problem? If this is the normal behavior, a more
> informative error message would be helpful.
>
> Thanks in advance for any help,
> Antonio
>
> PS: Let me congratulate with the lmfit authors as this is a really amazing
> piece of software that should be included in scipy IMHO!
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> <a
> href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a>

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

josef.pktd
On Wed, Jun 18, 2014 at 7:25 PM, Sturla Molden <[hidden email]> wrote:

> lmfit is not a part of SciPy, but as long as it's scientific use of Python
> and no other appropriate list, I don't think anyone will mind. I cannot
> answer your question, though.
>
> Sturla
>
> Antonino Ingargiola <[hidden email]> wrote:
>> Hi,
>>
>> I'm not sure this is the right list for this kind of issues.
>>
>> I'm using lmfit to solve a non-linear least square problem to fit a "noisy"
>> histogram to a model that is a single exponential decay convoluted with a
>> "sharp" peak and a baseline noise floor.
>>
>> The problem is that unless I set constrains for the 'offset' parameter the
>> confidence interval calculation fails in a strange way in scipy brentq
>> function.
>>
>> You can see the full example notebook and the error message in this ipython
>> notebook:
>>
>> <a
>> href="http://nbviewer.ipython.org/urls/gist.githubusercontent.com/tritemio/901c2eb2a43775e81844/raw/755eede8f170cc33b96316601286ae534e901c49/lmfit%20issue">http://nbviewer.ipython.org/urls/gist.githubusercontent.com/tritemio/901c2eb2a43775e81844/raw/755eede8f170cc33b96316601286ae534e901c49/lmfit%20issue</a>
>>
>> What can I do to fix the problem? If this is the normal behavior, a more
>> informative error message would be helpful.

http://stackoverflow.com/questions/20619156/python-brentq-problems/20625362

I have no idea about the details in lmfit, but it's not an infrequent
problem also outside of lmfit.

Josef

>>
>> Thanks in advance for any help,
>> Antonio
>>
>> PS: Let me congratulate with the lmfit authors as this is a really amazing
>> piece of software that should be included in scipy IMHO!
>>
>> _______________________________________________
>> SciPy-User mailing list
>> [hidden email]
>> <a
>> href="http://mail.scipy.org/mailman/listinfo/scipy-user">http://mail.scipy.org/mailman/listinfo/scipy-user</a>
>
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> http://mail.scipy.org/mailman/listinfo/scipy-user
_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Matt Newville
In reply to this post by Antonino Ingargiola
Hi Antonino,

On Wed, Jun 18, 2014 at 6:12 PM, Antonino Ingargiola <[hidden email]> wrote:
>
> Hi,
>
> I'm not sure this is the right list for this kind of issues.


It's OK with me.  Using github issues would be reasonable too.
 
>
>
> I'm using lmfit to solve a non-linear least square problem to fit a "noisy" histogram to a model that is a single exponential decay convoluted with a "sharp" peak and a baseline noise floor.
>
> The problem is that unless I set constrains for the 'offset' parameter the confidence interval calculation fails in a strange way in scipy brentq function.
>
> You can see the full example notebook and the error message in this ipython notebook:
>
> http://nbviewer.ipython.org/urls/gist.githubusercontent.com/tritemio/901c2eb2a43775e81844/raw/755eede8f170cc33b96316601286ae534e901c49/lmfit%20issue
>
> What can I do to fix the problem? If this is the normal behavior, a more informative error message would be helpful.

I haven't tried your example, but it's possible that the development version on github.com fixes this -- there was a similar problem with  conf_interval() failing when calling brentq() that was fixed relatively recently (slow devel cycle),  but I'm not certain that the fix there will work here.

I would guess that the essential problem is related to the fact that your also getting NaNs from the covariance matrix.  Again, without actually trying your problem, I believe that this is probably related to your use of

     pos_range = (x - offset) >= 0

to select the data range for which y!=0.   This could make the fit insensitive to small changes (ie, well below 1) to offset, making it hard to compute the partial derivative for it.  This could, in turn, make the covariance matrix hard to get (which is what you are seeing).    Without finite, numerical values for stderr from the covariance matrix, it's hard for conf_interval() to make good guesses about what range of values to search around. 

In short, it's a challenge to  have variables that are used essentially as integers.  Still, even if it doesn't work well, it shouldn't be failing this way. 

My advice would be to try the devel version of lmfit.  If it's still failing, try to construct a smaller failing example (or make something complete that can be run as a non-IPython script), and start an Issue.

Cheers,

--Matt Newville


_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Antonino Ingargiola
Hi Matt,

I reply inline...

On Wed, Jun 18, 2014 at 8:19 PM, Matt Newville <[hidden email]> wrote:
[cut]
I haven't tried your example, but it's possible that the development version on github.com fixes this -- there was a similar problem with  conf_interval() failing when calling brentq() that was fixed relatively recently (slow devel cycle),  but I'm not certain that the fix there will work here.

Indeed, upgrading to the latest version solved the error. Now I get the CI in both cases.

I would guess that the essential problem is related to the fact that your also getting NaNs from the covariance matrix.  Again, without actually trying your problem, I believe that this is probably related to your use of

     pos_range = (x - offset) >= 0

to select the data range for which y!=0.   This could make the fit insensitive to small changes (ie, well below 1) to offset, making it hard to compute the partial derivative for it.  This could, in turn, make the covariance matrix hard to get (which is what you are seeing).    Without finite, numerical values for stderr from the covariance matrix, it's hard for conf_interval() to make good guesses about what range of values to search around. 


Thanks for the explanation, now I understand much better the problem. I have a model function with a discontinuity in the origin (i.e. exp(-x - x0) for x > x0 else 0). If I sample it with a step dx, I will always have a problem when x0 changes less than dx. Is there any known trick I can use to avoid this problem?
 
In short, it's a challenge to  have variables that are used essentially as integers.  Still, even if it doesn't work well, it shouldn't be failing this way.  

Yes, but now I'm curious. How did you get around the problem? Are the confidence interval still accurate in this case? I assume that if the CI for the offset (x0) is < dx they don't make to much sense, is it right?
 
Thanks,
Antonio

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Matt Newville
Hi Antonio,


On Fri, Jun 20, 2014 at 3:37 PM, Antonino Ingargiola <[hidden email]> wrote:
Hi Matt,

I reply inline...

On Wed, Jun 18, 2014 at 8:19 PM, Matt Newville <[hidden email]> wrote:
[cut]

I haven't tried your example, but it's possible that the development version on github.com fixes this -- there was a similar problem with  conf_interval() failing when calling brentq() that was fixed relatively recently (slow devel cycle),  but I'm not certain that the fix there will work here.

Indeed, upgrading to the latest version solved the error. Now I get the CI in both cases.

OK, good.  Now, whether those CI are reliable...
 

I would guess that the essential problem is related to the fact that your also getting NaNs from the covariance matrix.  Again, without actually trying your problem, I believe that this is probably related to your use of

     pos_range = (x - offset) >= 0

to select the data range for which y!=0.   This could make the fit insensitive to small changes (ie, well below 1) to offset, making it hard to compute the partial derivative for it.  This could, in turn, make the covariance matrix hard to get (which is what you are seeing).    Without finite, numerical values for stderr from the covariance matrix, it's hard for conf_interval() to make good guesses about what range of values to search around. 


Thanks for the explanation, now I understand much better the problem. I have a model function with a discontinuity in the origin (i.e. exp(-x - x0) for x > x0 else 0). If I sample it with a step dx, I will always have a problem when x0 changes less than dx. Is there any known trick I can use to avoid this problem?

I'm not sure that there is a robust way to have the range of data considered in the fitting to be a parameter.  I might suggest (but haven't looked at your situation in great detail) to have you consider using an "offset" that shifts the origin for the model, then interpolate that onto the grid of data.  That way, you might be able to set a fit range in the data coordinates before the fit, and not have it change.  The model can be shifted in "x" (assuming there is such a thing -- your data appears to have an obvious "x" axis), but is fitting a fixed data range.  Again, I'm not sure that would fix all problems, but it might help.
 
In short, it's a challenge to  have variables that are used essentially as integers.  Still, even if it doesn't work well, it shouldn't be failing this way.  

Yes, but now I'm curious. How did you get around the problem?

Least-squares fitting with values used as integers are going to have poorly defined derivatives, and are bound to cause problems, at least sometimes.  I don't know of a general solution, but perhaps someone does.
 
Are the confidence interval still accurate in this case?

Practically speaking, if the fit with leastsq() is not sensitive to one of the variables and results in NaNs in the covariance matrix, the confidence_intervals aren't going to be easy to find.  The error you were seeing being raised by brentq() was from conf_interval() trying to find a suitable range of values to explore -- it has to find values that make the fit worse by 1-sigma, and so forth.  If there are NaNs in the covariance matrix, it won't be able determine whether one fit is worse than another, and won't really work to find better confidence intervals than those from the covariance matrix. 

 
I assume that if the CI for the offset (x0) is < dx they don't make to much sense, is it right?


I think that can be OK for a continuous/analog value, but will fail for a discrete value when dx is less than the step in discrete values.

Hope that helps,

--Matt Newville

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Antonino Ingargiola
Hi Matt,


On Fri, Jun 20, 2014 at 2:12 PM, Matt Newville <[hidden email]> wrote:
[cut]
Thanks for the explanation, now I understand much better the problem. I have a model function with a discontinuity in the origin (i.e. exp(-x - x0) for x > x0 else 0). If I sample it with a step dx, I will always have a problem when x0 changes less than dx. Is there any known trick I can use to avoid this problem?

I'm not sure that there is a robust way to have the range of data considered in the fitting to be a parameter.  I might suggest (but haven't looked at your situation in great detail) to have you consider using an "offset" that shifts the origin for the model, then interpolate that onto the grid of data.  That way, you might be able to set a fit range in the data coordinates before the fit, and not have it change.  The model can be shifted in "x" (assuming there is such a thing -- your data appears to have an obvious "x" axis), but is fitting a fixed data range.  Again, I'm not sure that would fix all problems, but it might help.

Unless I misunderstood, I do already what you suggest. My function is exp(- x + x0) (note that I wrote "-x0" before by mistake) and x0 is "continuous", regardless of the x discretization. The problem is that the function is 0 for x < x0 and therefore there is a discontinuity at x=x0. When the function is evaluated on the save discrete x arrays, changing smoothly x0 does not result in a smooth translation of the function.
 
  
In short, it's a challenge to  have variables that are used essentially as integers.  Still, even if it doesn't work well, it shouldn't be failing this way.  

Yes, but now I'm curious. How did you get around the problem?

Least-squares fitting with values used as integers are going to have poorly defined derivatives, and are bound to cause problems, at least sometimes.  I don't know of a general solution, but perhaps someone does.

I think that fitting a function with a discontinuity using an offset (x translation) as a parameter should be a quite common problem. 

Now that I think about this, maybe, constraining the offset variable within a single step of x will result in a smooth behavior and would allow to find the offset with accuracy of a fraction of the discretization step. I'll try...

Thinking loud, would be nice to have a 2-step fit. In step 1, you find the offset varying it with the same discretization of the x axis. In the second step, you vary the offset only within one x bin to find the fractional part.

Does this makes sense? Any idea if would it be feasible with current lmfit/scipy?

Thanks,
Antonio

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Matt Newville
Hi Antonio,


On Fri, Jun 20, 2014 at 6:36 PM, Antonino Ingargiola <[hidden email]> wrote:
Hi Matt,


On Fri, Jun 20, 2014 at 2:12 PM, Matt Newville <[hidden email]> wrote:
[cut]
Thanks for the explanation, now I understand much better the problem. I have a model function with a discontinuity in the origin (i.e. exp(-x - x0) for x > x0 else 0). If I sample it with a step dx, I will always have a problem when x0 changes less than dx. Is there any known trick I can use to avoid this problem?

I'm not sure that there is a robust way to have the range of data considered in the fitting to be a parameter.  I might suggest (but haven't looked at your situation in great detail) to have you consider using an "offset" that shifts the origin for the model, then interpolate that onto the grid of data.  That way, you might be able to set a fit range in the data coordinates before the fit, and not have it change.  The model can be shifted in "x" (assuming there is such a thing -- your data appears to have an obvious "x" axis), but is fitting a fixed data range.  Again, I'm not sure that would fix all problems, but it might help.

Unless I misunderstood, I do already what you suggest. My function is exp(- x + x0) (note that I wrote "-x0" before by mistake) and x0 is "continuous", regardless of the x discretization. The problem is that the function is 0 for x < x0 and therefore there is a discontinuity at x=x0. When the function is evaluated on the save discrete x arrays, changing smoothly x0 does not result in a smooth translation of the function.

Ah, sorry, I see that now.   I think I must have missed the use of "offset" in   "np.exp(-(x[pos_range] - offset)/tau)" earlier.      Yes, I think that should work, unless I'm missing something else....

Still, the main issue is getting NaNs in the covariance matrix for the "ampl" and "offset" parameters, which is especially strange since the resulting fit with best-fit values looks pretty reasonable.    You might try temporarily simplifying the model (turn weights on/off, turn convolution step on/off) and/or printing values in the residual or model function to see if you can figure out what conditions cause that to happen.   

A different (untested) guess: exponential decays are often surprisingly difficult for leastsq().   You might try a different algorithm (say, Nelder-Mead) as a first pass, then use those results as starting values for leastsq() (as that will estimate uncertainties).  I hear people who have good success with this approach when fitting noisy exponentially decaying data.

  
In short, it's a challenge to  have variables that are used essentially as integers.  Still, even if it doesn't work well, it shouldn't be failing this way.  

Yes, but now I'm curious. How did you get around the problem?

Least-squares fitting with values used as integers are going to have poorly defined derivatives, and are bound to cause problems, at least sometimes.  I don't know of a general solution, but perhaps someone does.

I think that fitting a function with a discontinuity using an offset (x translation) as a parameter should be a quite common problem. 

Yes, I do this fairly often myself (including, like you, forcing the model for x < x0 to zero).   I think I must have missed the use of "offset" as both the selection of the valid range and the shift itself. 
 
Now that I think about this, maybe, constraining the offset variable within a single step of x will result in a smooth behavior and would allow to find the offset with accuracy of a fraction of the discretization step. I'll try...

Thinking loud, would be nice to have a 2-step fit. In step 1, you find the offset varying it with the same discretization of the x axis. In the second step, you vary the offset only within one x bin to find the fractional part.

Does this makes sense? Any idea if would it be feasible with current lmfit/scipy?


I don't know how  feasible that would be, but I hope it is not necessary.
 
--Matt Newville

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|

Re: lmfit: confidence intervals issue

Antonino Ingargiola
Hi Matt,


On Sat, Jun 21, 2014 at 5:39 AM, Matt Newville <[hidden email]> wrote:
Hi Antonio,


On Fri, Jun 20, 2014 at 6:36 PM, Antonino Ingargiola <[hidden email]> wrote:
Hi Matt,


On Fri, Jun 20, 2014 at 2:12 PM, Matt Newville <[hidden email]> wrote:
[cut]
Thanks for the explanation, now I understand much better the problem. I have a model function with a discontinuity in the origin (i.e. exp(-x - x0) for x > x0 else 0). If I sample it with a step dx, I will always have a problem when x0 changes less than dx. Is there any known trick I can use to avoid this problem?

I'm not sure that there is a robust way to have the range of data considered in the fitting to be a parameter.  I might suggest (but haven't looked at your situation in great detail) to have you consider using an "offset" that shifts the origin for the model, then interpolate that onto the grid of data.  That way, you might be able to set a fit range in the data coordinates before the fit, and not have it change.  The model can be shifted in "x" (assuming there is such a thing -- your data appears to have an obvious "x" axis), but is fitting a fixed data range.  Again, I'm not sure that would fix all problems, but it might help.

Unless I misunderstood, I do already what you suggest. My function is exp(- x + x0) (note that I wrote "-x0" before by mistake) and x0 is "continuous", regardless of the x discretization. The problem is that the function is 0 for x < x0 and therefore there is a discontinuity at x=x0. When the function is evaluated on the save discrete x arrays, changing smoothly x0 does not result in a smooth translation of the function.

Ah, sorry, I see that now.   I think I must have missed the use of "offset" in   "np.exp(-(x[pos_range] - offset)/tau)" earlier.      Yes, I think that should work, unless I'm missing something else....

Still, the main issue is getting NaNs in the covariance matrix for the "ampl" and "offset" parameters, which is especially strange since the resulting fit with best-fit values looks pretty reasonable.    You might try temporarily simplifying the model (turn weights on/off, turn convolution step on/off) and/or printing values in the residual or model function to see if you can figure out what conditions cause that to happen.    

I suspect that the pseudo-periodic behaviour or the residuals as a function of the offset caused by the time axis discretization is causing problems here. BTW how can I see the covariance matrix?
 
A different (untested) guess: exponential decays are often surprisingly difficult for leastsq().   You might try a different algorithm (say, Nelder-Mead) as a first pass, then use those results as starting values for leastsq() (as that will estimate uncertainties).  I hear people who have good success with this approach when fitting noisy exponentially decaying data.

Oh, well, that's a very good suggestion. In the few trials I did Nelder-Mead works, then leastsq does not move the solution anymore but I can get the CI. Neat.

As I told you the original error is gone when updating to current master. However I find some new combinations of data and initial parameters that give errors when computing the CI. The errors are:

    AttributeError: 'int' object has no attribute 'copy'

and

    ValueError: f(a) and f(b) must have different signs

I opened an issue to track them:


Best,
Antonio

_______________________________________________
SciPy-User mailing list
[hidden email]
http://mail.scipy.org/mailman/listinfo/scipy-user