Hi,
I encountered with some difficulty in preparing the nonlin.function for the
gnm function. Please help me.
I want to use the gnm function in the "gnm" package to fit the following
spherical function
y = a + b (1.5*x - 0.5 *(x/c)) if x< c
y = a + b if x >= c
I understand that the fit.variogram function in the gstat function was
designed for fitting variogram, but it seems to be questionable because of
the normality assumption on the residuals. It is interested that whether
the fitted model would be changed as another distribution was assumed.
Seems to me, the gnm function in the gnm/R package could be used for the
this. However, I failed to customize a proper "nonlin" object for the
aforementioned formula after a few hours of try.
The following are my crappy code defining the nonlin object.
-------------------------------------------------------------
myfun<-function(x) {
list( predictors = list(nuggest = .1, psill = 1, range=5), variables =
list(substitute(x)),
term <- function(predLabels, varLabels) {
paste0("pred=", predLabels[1], "+", predLabels[2], "*(1.0 - exp(-",
varLabels[1], "/", predLabels[3], ")")
# paste0("pred[which(", varLabels[1], ">", predLabels[3], ")]=",
predLabels[3], "+", predLabels[2])
# paste0("pred")
}
)
}
class(myfun) <- "nonlin"
-------------------------------------------------------------
In addition, I have also tried the gnls function in the nlme package but
did not find any parameters in function would allow to specify a
distribution for the residuals.
I appreciate any help, suggestion, comment on my question.Thank you very
much.
Best regards,
Zhiqiu
The following code are largely irrelevant to my question, but it could be
used to generate a semivarogram table for the analysis.
#The Sph model in SP analysis
sph.fun<-function(dist, nugget, psill, range) {
res<-nugget+psill*(1.5*dist/range-.5*(dist/range)^3)
id.far=which(dist>=range)
if (NROW(id.far)>0) res[id.far]= sum(nugget, psill)
res
}
#simulate semivariogram data
library(gstat)
dist=1:100
dat=data.frame(np=300.0, dist=(dist)+0.001,
gamma=sph.fun(dist, nugget = 1, psill = 9, range =
20)+rnorm(NROW(dist))*1e-3,
dir.hor=0.0, dir.ver=0.0, id=as.factor(1))
v=dat;
class(v)<-"gstatVariogram"
#By defalut, WLS is used and the weight is w=np/dist^2
ft1=fit.variogram(v, vgm(nugget = 1, mod="Sph", psill = 1, range = 5))
#use the nls to fit the data
ft2=nls(gamma ~ sph.fun(dist, nugget, psill, range),
data = dat,
start = list(nugget=1, psill= 1, range = 5),
weight=np/dist^2,
algorithm="port",
lower=c(0, 0, 0),
trace=F)
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org 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.
