Re: [R] nls

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
Date: Sat 23 Sep 2006 - 10:39:03 GMT

>>>>> "SpG" == Spencer Graves <spencer.graves@pdf.com> >>>>> on Sat, 23 Sep 2006 00:52:30 -0700 writes:

    SpG> I used "debug" to walk through your example line by line,  I found 
    SpG> that the error message was misleading.  By making 
    SpG> as.vector(semivariance) and as.vector(h)  columns of a data.frame, I got 
    SpG> it to work.  My revised code appears below. 

    SpG> Thanks for providing a self-contained and relatively simple 
    SpG> example. Without it, I could have done little more that suggest     SpG> 'traceback' (which provided zero information when I tried it) and 'debug'.

Thank you, Spencer, for picking this up. Your solution is probably fine for Mzabalazo Ngwenya.

However your note about the unhelpful error message made me dig a bit further.
My current diagnosis: nls() is at first doing a *multivariate* model in this case which model.frame() happily does. Later in nlsModel(); qr() does not deal with such multivariateness and hence produces a 'QR' result which does not "fit" the other parts of the computations.

I think we should decide if it is easy enough to make nls() work with multivariate models and do it if it's "easy" and otherwise disallow it consciously -- both via a better error message and the help page.

Well, for R 2.4.0 it should be possible to do the latter and I think we should...
This is now becoming an "R-devel" rather than "R-help" topic..

Martin

    SpG> Hope this helps.
    SpG> Spencer Graves

    SpG> fit.gaus<-function(coordinates,values,guess.c0,guess.c1,guess.a)
    SpG> {
    SpG> long<-rep(coordinates[,1],each=length(coordinates[,1]))
       
    SpG> lag.long<-t(matrix(long,nrow=length(coordinates[,1]),byrow=TRUE))
    SpG> dif.long <-(lag.long-t(lag.long))^2
    SpG> lat <-rep(coordinates[,2],each=length(coordinates[,2]))
    SpG> lag.lat<-t(matrix(lat,nrow=length(coordinates[,2]),byrow=TRUE))
    SpG> dif.lat <-(lag.lat-t(lag.lat))^2     SpG> h <-sqrt(dif.long+dif.lat)
        
                                       
    SpG> if( length(values[1,])>1)
    SpG> {
    SpG> y.m <-apply(values,1,sum,na.rm=TRUE)
    SpG> y.m <-as.matrix(y.m)
    SpG> y.mod <-(1/length(values[1,]))*(y.m)
    SpG> }
    SpG> else
    SpG> {

    SpG> y.mod <-as.matrix(values)
    SpG> }
    SpG> semi <-rep(y.mod,each=length(y.mod))
    SpG> mat1<-t(matrix(semi,nrow=length(y.mod),byrow=TRUE))
    SpG> mat2<-t(mat1)
    SpG> semivariance <-(1/2)*(mat1-mat2)^2

    SpG> #        model <-semivariance ~c0+c1*(1-exp(-(h^2)/a^2))
    SpG> DF <- data.frame(smvar=as.vector(semivariance),
    SpG> h.=as.vector(h) )
    SpG> mdl <- smvar~c0+c1*(1-exp(-(h.^2)/a^2))
        
    SpG> #        parameters <-nls(model,start =
    SpG> #list(c0=guess.c0,c1=guess.c1,a=guess.a),trace=TRUE)
    SpG> parameters <-nls(mdl,start =
    SpG> list(c0=guess.c0,c1=guess.c1,a=guess.a),trace=TRUE,
    SpG> data=DF )
    SpG> results <-summary(parameters)

    SpG> print(results)
    SpG> }
    SpG> ################################
    SpG> mzabalazo ngwenya wrote:
>> Hello everyone !
>>
>> I am trying to write a short program to estimate semivariogram
>> parameters. But I keep running into a problem when using the nls
>> function.
>>
>> Could you please shed some light. I have put a sample of one of the
>> codes and ran a short example so you see what I mean.
>>
>> -----------------
>>
>>
>>
>> fit.gaus<-function(coordinates,values,guess.c0,guess.c1,guess.a)
>> {
>> long<-rep(coordinates[,1],each=length(coordinates[,1]))
>>
>> lag.long<-t(matrix(long,nrow=length(coordinates[,1]),byrow=TRUE))
>> dif.long <-(lag.long-t(lag.long))^2
>> lat <-rep(coordinates[,2],each=length(coordinates[,2]))
>> lag.lat<-t(matrix(lat,nrow=length(coordinates[,2]),byrow=TRUE))
>> dif.lat <-(lag.lat-t(lag.lat))^2
>> h <-sqrt(dif.long+dif.lat)
>>
>>
>> if( length(values[1,])>1)
>> {
>> y.m <-apply(values,1,sum,na.rm=TRUE)
>> y.m <-as.matrix(y.m)
>> y.mod <-(1/length(values[1,]))*(y.m)
>> }
>> else
>> {
>> y.mod <-as.matrix(values)
>> }
>>
>> semi <-rep(y.mod,each=length(y.mod))
>> mat1<-t(matrix(semi,nrow=length(y.mod),byrow=TRUE))
>> mat2<-t(mat1)
>> semivariance <-(1/2)*(mat1-mat2)^2
>>
>> model <-semivariance ~c0+c1*(1-exp(-(h^2)/a^2))
>> parameters <-nls(model,start =
>> list(c0=guess.c0,c1=guess.c1,a=guess.a),trace=TRUE)
>> results <-summary(parameters)
>> print(results)
>> }
>> --------------------------
>>
>>
    >>> don <-matrix(c(2,3,9,6,5,2,7,9,5,3),5,2)
    >>> don
    >>> 

>> [,1] [,2]
>> [1,] 2 2
>> [2,] 3 7
>> [3,] 9 9
>> [4,] 6 5
>> [5,] 5 3
>>
    >>> data <-matrix(c(3,4,2,4,6))
    >>> data
    >>> 

>> [,1]
>> [1,] 3
>> [2,] 4
>> [3,] 2
>> [4,] 4
>> [5,] 6
>>
    >>> fit.gaus(don,data,2,3,5)
    >>> 

>> [,1] [,2] [,3] [,4] [,5]
>> [1,] 0.000000 5.099020 9.899495 5.000000 3.162278
>> [2,] 5.099020 0.000000 6.324555 3.605551 4.472136
>> [3,] 9.899495 6.324555 0.000000 5.000000 7.211103
>> [4,] 5.000000 3.605551 5.000000 0.000000 2.236068
>> [5,] 3.162278 4.472136 7.211103 2.236068 0.000000
>> 178.9113 : 2 3 5
>> Error in qr.qty(QR, resid) : 'qr' and 'y' must have the same number of
>> rows
>>
>>
>> ______________________________________________
>> R-help@stat.math.ethz.ch mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-help
>> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
    SpG> ______________________________________________
    SpG> R-help@stat.math.ethz.ch mailing list
    SpG> https://stat.ethz.ch/mailman/listinfo/r-help
    SpG> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html     SpG> and provide commented, minimal, self-contained, reproducible code.

R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat Sep 23 20:59:43 2006

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.1.8, at Sat 23 Sep 2006 - 11:30:06 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.