[R] nlrob problem
Martin Maechler
maechler at stat.math.ethz.ch
Fri Dec 23 17:21:22 CET 2011
>>>>> Pascal A Niklaus <pascal.niklaus at ieu.uzh.ch>
>>>>> on Mon, 19 Dec 2011 20:46:57 +0100 writes:
> Dear all, I am not sure if this mail is for R-help or
> should be sent to R-devel instead, and therefore post to
> both.
> While using nlrob from package 'robustbase', I ran into
> the following problem:
> For psi-functions that can become zero
> (e.g. psi.bisquare), weights in the internal call to nls
> can become zero. Example:
> d <- data.frame(x=1:5,y=c(2,3,5,10,9))
> d.nlrob <- nlrob(y ~ k*x,start=list(k=1),data=d,psi=psi.bisquare)
> I think the problem is as follows: After the call to nls, a weighted
> residual vector is calculated by dividing by sqrt(w). This generates
> NaN's for weights of zero, which leads to problems in the subsequent
> call to nls during the next iteration:
> for (iiter in 1:maxit) {
> ...
> w <- psi(resid/Scale, ...)
> ...
> data$w <- sqrt(w)
> ...
> out <- nls( ....., data=data ....... )
> ...
> resid <- -residuals(out)/sqrt(w) # NaN for w=0
> ...
> }
> I wonder whether this problem shouldn't be dealt with by setting 'w' to
> 0 for the NaN cases.
> I can get a fit by calling nlrob with na.action=na.exclude, but I'd have
> intuitively assumed that na.action applies to the NAs in the data set
> passed to nlrob, not to the one internally generated by nlrob and passed
> to nls.
You are right.
The next version of robustbase (0.8-1) will have this fixed.
Note however, that for me, in your example and in other more
interesting ones, as soon as I use a redescending psi() function,
the IRLS iterations do *not* converge...
but rather strangely ``cycle'' around the true solution.
As a redescending psi() has a non-convex rho(), and non-linear
problems depend quite a bit on "feasible"/"good" starting
values, I currently think I would discourage from using such psi().
As others, some possibly more expert in robust non-linear fitting,
may be interested, and have interesting feedback,
I'm crossposting this reply to the R-SIG-robust mailing
list.
(Further replies should typically only go there!)
> The second 'issue' is that the weights are passed as 'w', whereas the
> documentation of 'nls' says weights should be given as 'weights' :
> data: an optional data frame in which to evaluate the variables in
> ‘formula’ and ‘weights’. Can also be a list or an
> environment, but not a matrix.
> I think it would be good to mention in the documentation of 'nls' that
> weights can be passed as 'w' as well.
You have overlooked that nlrob uses a different *formula* for nls()
--- nowadays unnecessarily ---
and that formula has a 'w'.
So this part of your report is not correct.
Martin Maechler, ETH Zurich
> Pascal Niklaus
More information about the R-help
mailing list