[R] nls

[Martin Maechler]

>>>>> "SpG" == Spencer Graves <spencer.graves at 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 at 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 at 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.