robust fit

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

robust fit

x.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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Matthieu Brucher-2
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.
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



--
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Robert Kern-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

josef.pktd
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

x.piter
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

josef.pktd
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

josef.pktd
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

josef.pktd
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Charles R Harris
In reply to this post by josef.pktd


On Tue, May 31, 2011 at 9:39 AM, <[hidden email]> wrote:
On Tue, May 31, 2011 at 10:21 AM, Sturla Molden <[hidden email]> wrote:
> Den <a href="tel:31.05.2011%2016" value="+13105201116">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).


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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

Sturla Molden-2
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
Reply | Threaded
Open this post in threaded view
|

Re: robust fit

josef.pktd
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