[Rd] Using the nls package

[Telford Tendys]

On Mon, Jul 31, 2000 at 07:20:32AM +0100, Prof Brian D Ripley wrote:
> 
> How can you be sure?  What are `the same numbers'?  The parameters
> unchanged to machine precision, yes, but that's not what you print.

The only safe way I can think of would be to go for machine precision
as the default and then have some user-settable threshold for:

	( current residual / previous residual )

> > Yes, this is understandable but if most people are like me (questionable
> > assumption I admit) they try a few artificial examples first to get
> > the general idea of how these things work. I guess there is a certain
> > illogical assumption that if I try what seems to me to be an obviously
> > easy problem and it cannot handle it then something must be wrong and
> > trying harder problems is not going to help.
> 
> Perhaps you need to read up on the algorithmic theory so you are not so
> easily confused.  NLS problems with an exact-fit solution are a different
> sub-class with different effective algorithms.

Undoubtably there is a lot more reading that I can do to improve
my understanding. I'm just making a suggestion based on the observation
that the algorithm actually found the right answer very quickly.

> > SQUARES of the modulus of the residuals. Using the modulus itself
> > would be similar to using the absolute value of the residuals in the
> > non-complex case (not what we want I believe).
> 
> Sum of squares is common, and ises for complex linear regression in those
> cases where it makes sense (e.g. in time series regression one
> FFT-ed series on others).

As it turns out, I am working with FFT data, I am looking for a transfer
function that will fit some admittance measurements.

> > In either case, the documentation should mention that complex
> > numbers are not supported.
> 
> Why should documentation mention everything the function does not do? There
> are *lots* of places in R where numeric quantities must be real.  Linear
> regression (lm) for a start.  What nls says is
> 
>  formula: a nonlinear model formula including variables and parameters
> 
> which is none too explicit, but the reference (Bates & Watts) is,
> so I assume you did not read it.

The UTS library does have a copy but it isn't due back till September.

I searched a few of the local bookshops (over the internet) and
Diana Bates (with her children's stories) came up infinitely more often
than Douglas Bates and his regression. I will have to search the
University Co-op bookshop with my hands because the recent introduction of
GST has blown up their webpage.

Douglas, if you are reading this, what say you help me understand the
code in nls.R and nls.c and I'll promise to buy your book (though it could
take me a while to mail-order it in) ?

> A pretty good reasonable assumption is that statistics does not work
> in complexe numbers unless explicitly stated.

Consider also the case for fitting higher dimensional REAL data.
Admittedly, this is a contrived example but you may have an experiment
with one input variable and three output variables. You know from theory
that the output variables are all related to each other by some function
but you need to find particular parameters in that function:

> d <- data.frame( x=1:10 )
> d$y <- eval( expression( outer( x, 1:3 )), d )
> d$y <- d$y + rnorm( length( d$y ))
> library( nls )
> w <- nls( y ~ I( outer( x, P1 * 1:3 ) + P2 ), d, start=list( P1=3, P2=3 ))
Error in qr.qty(QR, resid) : qr and y must have the same number of rows

If the nls system supported these `wide' problems properly then complex
number support could just be recast as a special case of the same.
>From a mathematical point of view, this sort of regression can be recast
into a normal regression but with more sample points (and each sample
point repeated N times) and with a bogus input variable that gives the
column number. Thus, the same fundamental algorithm must work for both
cases if the structure of the data is suitably munged around.

I've been looking through the code in nls.c and nls.R quite closely and
I can see a lot of stuff which seems to be catering to cases of various
dimensional objects. I cannot fully comprehend what it does and I cannot
get a high dimensional case (such as the above) to work properly.

The crux of it seems to be that qr() needs a matrix to work on
so it uses as.matrix() to coerce its input. If you feed it a vector
you get a column matrix (sounds good), if you feed it a matrix you
have no problem, if you feed it a high dimensional tensor you get
a column vector (which surprised me), and qr() goes with that:

> dim( qr( 1:10 )$qr )
[1] 10  1
> dim( qr( outer( 1:3, 1:5 ))$qr )
[1] 3 5
> dim( qr( outer( 1:2, outer( 1:3, 1:4 )))$qr )
[1] 24  1

The result is that the nls model comes to grief when it feeds
the three dimensional tensor into qr() then what it gets back
(i.e. a long, thin column matrix) does not match up with anything.

>  That's why in R
> 
> > sqrt(-1)
> [1] NaN
> Warning message: 
> NaNs produced in: sqrt(-1) 

But sqrt() DOES support complex numbers... so the assumption is that
most users will not expect complex numbers to crop up in their calculations
unless they were deliberately put there, but once deliberatly put there,
they work properly.

I would be happy to see two examples added to the documentation for sqrt,
those being sqrt(-1) and sqrt(-1+0i).

	- Tel



-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._