Hi all.
Can anybody point me a direction how to make a robust fit of nonlinear function using leastsq. Maybe someone have seen ready function doing this. Thanks. Petro. _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Hi,
It seems to me that least squares cannot lead to a ribust fit in the usual sense. You may want to implement your own robust cost function instead of using a L2 norm. Matthieu
2011/5/30 Piter_ <[hidden email]> Hi all. -- Information System Engineer, Ph.D. Blog: http://matt.eifelle.com LinkedIn: http://www.linkedin.com/in/matthieubrucher _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by x.piter
On Mon, May 30, 2011 at 05:58, Piter_ <[hidden email]> wrote:
> Hi all. > Can anybody point me a direction how to make a robust fit of nonlinear > function using leastsq. > Maybe someone have seen ready function doing this. A variety of robust fits are implemented in statsmodels: http://statsmodels.sourceforge.net/trunk/rlm.html -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On Mon, May 30, 2011 at 11:12 AM, Robert Kern <[hidden email]> wrote:
> On Mon, May 30, 2011 at 05:58, Piter_ <[hidden email]> wrote: >> Hi all. >> Can anybody point me a direction how to make a robust fit of nonlinear >> function using leastsq. >> Maybe someone have seen ready function doing this. > > A variety of robust fits are implemented in statsmodels: > > http://statsmodels.sourceforge.net/trunk/rlm.html Unfortunately, rlm includes only linear models, so unless the problem can be converted to a linear in parameters problem RLM will not help directly. I don't know of anything that would be immediately available for this. long answer: >From a quick Google search it seems the same iteratively reweighted regression method can be applied to non-linear models, but I didn't find an internet accessible paper, and I don't know the details about how easy it would be to add non-linear least squares instead of linear weighted least squares to statsmodels.robust. scipy.optimize.curvefit allows for weights, so it might be possible to down-weight observations with large residuals in non-linear least squares. Some ways that shouldn't be too difficult to implement: If there are clear outliers, it might be possible to identify them and remove or trim them, trimmed least squares. (Using optimize fmin with a robust loss function, might be possible, but I have no idea how well it works or how to get estimates for the covariance of the error estimates.) Since I'm not so familiar with these robust methods but know Maximum Likelihood estimation, what I would do is to assume that the errors come from a non-normal distribution, either a mixture model, if some observations might be generated by a different model, or assume that the errors are t-distributed. In the linear examples that I looked at, t-distributed maximum likelihood was very robust to outliers, (error distribution with heavy tails). It should also work quite easily (using GenericLikelihoodModel in statsmodels), if the non-linear model is well behaved and/or good starting values are available. I don't remember whether I have read this specific paper but it's top in my google search http://www.jstor.org/stable/2290063 (Cited by 490 in google) Robust Statistical Modeling Using the t Distribution Kenneth L. Lange, Roderick J. A. Little and Jeremy M. G. Taylor Abstract: "The t distribution provides a useful extension of the normal for statistical modeling of data sets involving errors with longer- than-normal tails. An analytical strategy based on maximum likelihood for a general model with multivariate t errors is suggested and applied to a variety of problems, including linear and nonlinear regression, robust estimation of the mean and covariance matrix with missing data, unbalanced multivariate repeated-measures data, multivariate modeling of pedigree data, and multi- variate nonlinear regression. The degrees of freedom parameter of the t distribution provides a convenient dimension for achieving robust statistical inference, with moderate increases in computational complexity for many models. Estimation of precision from asymptotic theory and the bootstrap is discussed, and graphical methods for checking the appropriateness of the t distribution are presented." Josef > > -- > Robert Kern > > "I have come to believe that the whole world is an enigma, a harmless > enigma that is made terrible by our own mad attempt to interpret it as > though it had an underlying truth." > -- Umberto Eco > _______________________________________________ > 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 |
In reply to this post by x.piter
I've done this for M-estimates using bisquare weights.
If the errors are error(x) = y(x) - fit(x) you want the minimize the sum of squares of weighted errors: werror(x) = weight(error(x)) * error(x) Provide one function that returns werror for each x[:,n],y[n] and another that returns the partial derivatives of werror[:] with respect to x[i,:]. The latter is the Jacobian that you set as Dfun keyword argument. I.e. if you have m parameters and n samples in the array x, the Jacobian is an m x n array. The last thing you need is an initial guess. This will depend on the problem, so its hard to give an advide. Just be beware that a least squares fit will not always work. If the you think deriving the Jacobian is tedious, consider using SymPy or just let leastsq approximate it (i.e. set Dfun to None). Sturla Den 30.05.2011 12:58, skrev Piter_: > Hi all. > Can anybody point me a direction how to make a robust fit of nonlinear > function using leastsq. > Maybe someone have seen ready function doing this. > Thanks. > Petro. > _______________________________________________ > 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 |
In reply to this post by Matthieu Brucher-2
Den 30.05.2011 13:19, skrev Matthieu Brucher:
> > It seems to me that least squares cannot lead to a ribust fit in the > usual sense. Leastsq is actually Levenberg-Marquardt (lmder and lmdif from MINPACK). It can be used to minimize any cost function if you provide the Jacobian. Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
Thanks for all answers.
I will look it through. Best. Petro _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Sturla Molden-2
Den 31.05.2011 16:10, skrev Sturla Molden:
> you want the minimize the sum of squares of weighted > errors: > > werror(x) = weight(error(x)) * error(x) One more thing: weight(error) is actually the square root of the robust weighting function, as we want to minimize the sum of robust_weight(error) * (error**2) Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Sturla Molden-2
On Tue, May 31, 2011 at 10:17 AM, Sturla Molden <[hidden email]> wrote:
> Den 30.05.2011 13:19, skrev Matthieu Brucher: >> >> It seems to me that least squares cannot lead to a ribust fit in the >> usual sense. > > Leastsq is actually Levenberg-Marquardt (lmder and lmdif from MINPACK). > It can > be used to minimize any cost function if you provide the Jacobian. Do you have an example for this? >From the docstring of lmdif and lmder it can only minimize sum of squares, e.g. c the purpose of lmdif is to minimize the sum of the squares of c m nonlinear functions in n variables by a modification of c the levenberg-marquardt algorithm. Josef > > Sturla > _______________________________________________ > 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 |
Den 31.05.2011 17:30, skrev [hidden email]:
> > Do you have an example for this? > > From the docstring of lmdif and lmder it can only minimize sum of squares, e.g. > > c the purpose of lmdif is to minimize the sum of the squares of > c m nonlinear functions in n variables by a modification of > c the levenberg-marquardt algorithm. > Huh? If it can minimize sum(x*x), then it can also minimize sum(w*x*x) by the transformation z = sqrt(w)*x. Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Sturla Molden-2
On Tue, May 31, 2011 at 10:21 AM, Sturla Molden <[hidden email]> wrote:
> Den 31.05.2011 16:10, skrev Sturla Molden: >> you want the minimize the sum of squares of weighted >> errors: >> >> werror(x) = weight(error(x)) * error(x) > > One more thing: weight(error) is actually the square root of the > robust weighting function, as we want to minimize the sum of > > robust_weight(error) * (error**2) Do you have a reference or a full example for this? Your description sounds relatively simple, and I guess now that we (statsmodels) can get the non-linear version with only small(ish) changes and a call to a WNLLS (curve_fit) instead of WLS (linear). Josef > > Sturla > > _______________________________________________ > 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 |
In reply to this post by Sturla Molden-2
On Tue, May 31, 2011 at 11:37 AM, Sturla Molden <[hidden email]> wrote:
> Den 31.05.2011 17:30, skrev [hidden email]: >> >> Do you have an example for this? >> > From the docstring of lmdif and lmder it can only minimize sum of squares, e.g. >> >> c the purpose of lmdif is to minimize the sum of the squares of >> c m nonlinear functions in n variables by a modification of >> c the levenberg-marquardt algorithm. >> > > Huh? > > If it can minimize sum(x*x), then it can also minimize sum(w*x*x) by the > transformation z = sqrt(w)*x. Yes, that's what scipy.optimize.curve_fit is doing, but it is still a sum of squares, your statement was *any* cost function (not just quadratic, sum of squares): > It can > be used to minimize any cost function if you provide the Jacobian. Just getting mislead by the definition of *any*. Josef > > Sturla > > > > _______________________________________________ > 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 |
In reply to this post by josef.pktd
Den 31.05.2011 17:39, skrev [hidden email]:
> Do you have a reference or a full example for this? Yes, I'll look one up for you. > Your description sounds relatively simple, and I guess now that we > (statsmodels) can get the non-linear version with only small(ish) > changes and a call to a WNLLS (curve_fit) instead of WLS (linear). > Yes and no. You have to provide an initial fit and derive the Jacobian (unless you are happy with an approximation). Also the "robust fit" (asked for here) is not just WNLLS, because the weights are a non-linear function of the residuals. You can solve this by iterative WNLLS, or provide the full Jacbobian directly to Levenberg-Marquardt. In the former case, it is just sqrt(w) times the Jacobian for the residuals. But you can also get the M-estimator from a single pass of Levenberg-Marquardt by using the chain rule on sqrt(w(e(x)))*e(x). So there are actually two ways of doing this with leastsq. Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by josef.pktd
On Tue, May 31, 2011 at 9:39 AM, <[hidden email]> wrote:
I've also had good results with Tukey's biweight. As Sturla says, it can be implemented as iterated weighted least squares. There is a whole class of robust methods along that line. The L_1 cost function can also be done that way, and the usual algorithm for the geometric median is one useful application. Chuck _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by josef.pktd
Den 31.05.2011 17:43, skrev [hidden email]:
> > Yes, that's what scipy.optimize.curve_fit is doing, but it is still a > sum of squares, your statement was *any* cost function (not just > quadratic, sum of squares): > I was just answering Matthieu's statement that least sqaures cannot lead to a robust fit in the usual sence. But in this case it can. Sorry for sloppy choise of wording. It cannot optimize any function. But it can optimize any function that can be expressed as a non-linear least-squares problem, and that includes robust M-estimators. Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
In reply to this post by Charles R Harris
Den 31.05.2011 17:57, skrev Charles R Harris:
> > > I've also had good results with Tukey's biweight That is the one I prefer too, particularly with a robust estimate of the standard error (MAD/0.6745). If the residuals are normally distributed, there is hardly any difference from the least squares fit. Sturla _______________________________________________ SciPy-User mailing list [hidden email] http://mail.scipy.org/mailman/listinfo/scipy-user |
On Tue, May 31, 2011 at 12:09 PM, Sturla Molden <[hidden email]> wrote:
> Den 31.05.2011 17:57, skrev Charles R Harris: >> >> >> I've also had good results with Tukey's biweight > > That is the one I prefer too, particularly with a robust estimate of the > standard error (MAD/0.6745). If the residuals are normally distributed, > there is hardly any difference from the least squares fit. just to add some advertising for statsmodels and Skipper's work http://statsmodels.sourceforge.net/devel/rlm.html#norms http://statsmodels.sourceforge.net/devel/rlm_techn1.html http://statsmodels.sourceforge.net/devel/generated/scikits.statsmodels.robust.scale.stand_mad.html#scikits.statsmodels.robust.scale.stand_mad hopefully we get the non-linear extension (if someone finds the time this summer). Josef > > Sturla > _______________________________________________ > 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 |
Free forum by Nabble | Edit this page |