Greetings R users, maybe there is someone who can help
me with this problem:
I define a function "optim.fun" and want as output the
sum of squared errors between predicted and measured
values, as follows:
optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
{
predR <-
(a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
abserr <- abs(ch2no-predR)
qnum <- quantile(abserr, probs=0.95, na.rm=T)
is.na(abserr) <- (abserr > qnum)
errsq <- sum(abserr^2, na.rm=T)
errsq
}
Then I want to optimize parameters a,b,d and E as to
minimize the function output with the following:
optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
ch2no, a, b, d, E), data=tower,
start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action na.exclude )
were nulldat=0
At this point I get the following error message:
Error in qr.default(.swts * attr(rhs, "gradient")) :
NA/NaN/Inf in foreign function call (arg 1)
Warning messages:
1: In if (na.rm) x <- x[!is.na(x)] else if
(any(is.na(x))) stop("missing values and NaN's not
allowed if 'na.rm' is FALSE") ... :
the condition has length > 1 and only the first
element will be used
(this warning is repeated 12 times)
Question: what does the error mean? What am I doing
wrong?
Thanks a bunch.
Nano
Jen, Germany
Max Planck for Biogeochemistry
____________________________________________________________________________________
[[elided Yahoo spam]]
The error message means that the gradient (first derivative of residual vector with respect to the parameter vector) is not possible to work with; calling the function qr on the gradient multiplied by the square root of the weight vector .swts (in your case all 1's) fails. If you want concrete advice it would be helpful to provide the commented, minimal, self-contained, reproducible code that the posting guide requests. what are the values of ST04, SM08b, ch2no, and tower? General comments: If your goal is to minimize sum( (observed - predicted)^2), the function you give nls to minimize (optim.fun in your case) should return the vector (observed - predicted), not the scalar sum( (observed - predicted)^2). You may want to see the nls.lm function in package minpack.lm for another way to deal with the problem. On Wed, 7 May 2008, Fernando Moyano wrote:> Greetings R users, maybe there is someone who can help > me with this problem: > > I define a function "optim.fun" and want as output the > sum of squared errors between predicted and measured > values, as follows: > > optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E) > { > predR <- > (a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13)))) > abserr <- abs(ch2no-predR) > qnum <- quantile(abserr, probs=0.95, na.rm=T) > > is.na(abserr) <- (abserr > qnum) > errsq <- sum(abserr^2, na.rm=T) > errsq > } > > Then I want to optimize parameters a,b,d and E as to > minimize the function output with the following: > > optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b, > ch2no, a, b, d, E), data=tower, > start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action > na.exclude ) > > were nulldat=0 > At this point I get the following error message: > > Error in qr.default(.swts * attr(rhs, "gradient")) : > NA/NaN/Inf in foreign function call (arg 1) > Warning messages: > 1: In if (na.rm) x <- x[!is.na(x)] else if > (any(is.na(x))) stop("missing values and NaN's not > allowed if 'na.rm' is FALSE") ... : > the condition has length > 1 and only the first > element will be used > (this warning is repeated 12 times) > > Question: what does the error mean? What am I doing > wrong? > Thanks a bunch. > Nano > Jen, Germany > Max Planck for Biogeochemistry > > > > > ____________________________________________________________________________________ > > [[elided Yahoo spam]] > > ______________________________________________ > R-help at 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. >
I solved the problem arising from using certain quantile values simply by
printing the iterations with the argument nprint. This gave me correct
estimates. I have no idea why.
----- Original Message ----
From: elnano <nanomoyano at yahoo.com>
To: r-help at r-project.org
Sent: Thursday, May 8, 2008 5:43:31 PM
Subject: Re: [R] function in nls argument
I've basically solved the problem using the nls.lm function from the
minpack.lm (thanks Katharine) with some modifications for ignoring residuals
above a given percentile. This is to avoid the strong influence of points
which push my modeled vs. measured values away from the 1:1 line.
I based it on the example given for nls.lm. Here it is:
R # soil respiration data
ST <- ST [!is.na(R)] # soil temeprature data. Had to remove na to make
nls.lm work
SM <- SM [!is.na(R)] # soil moisture data
R <- R [!is.na(R)]
q <- 0.95 # quantile
p <- c("a"=-0.003, "b"=0.13, "c"=0.50, E=400)
# model parameters with
estimated values
Rf <- function(ST, SM, a, b, c, E)
{
expr <-
expression((a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13)))))
eval(expr)
}
optim.f <- function(p, ST, SM, R, Rfcall)
{
res <- R - do.call("Rfcall", c(list(ST = ST, SM = SM),
as.list(p))) #
the nls.lm example divides this by sqrt(R), I don't know why. I removed
that.
abserr <- abs(res)
qnum <- quantile(abserr, probs=q, na.rm=T) # calculate the
"q"
quantile of the absolute errors
res[abserr > qnum]=0 # convert
residuals above qnum to 0
res
}
Rmodel<- nls.lm(par = p, fn = optim.f, Rfcall = Rf, ST = ST, SM = SM, R
= R)
summary(Rmodel)
The only error I still get is when using lower q values. A q of around 0.95
or less (depending on the dataset) gives me completely wrong parameter
estimates resulting in negative predicted values. Maybe someone has a
suggestion here.
Fernando
Katharine Mullen wrote:>
> The error message means that the gradient (first derivative of residual
> vector with respect to the parameter vector) is not possible to work with;
> calling the function qr on the gradient multiplied by the square root of
> the weight vector .swts (in your case all 1's) fails.
>
> If you want concrete advice it would be helpful to provide the commented,
> minimal, self-contained, reproducible code that the posting guide
> requests. what are the values of ST04, SM08b, ch2no, and tower?
>
> General comments: If your goal is to minimize sum( (observed -
> predicted)^2), the function you give nls to minimize (optim.fun in your
> case) should return the vector (observed - predicted), not the scalar sum(
> (observed - predicted)^2). You may want to see the nls.lm function in
> package minpack.lm for another way to deal with the problem.
>
> On Wed, 7 May 2008, Fernando Moyano wrote:
>
>> Greetings R users, maybe there is someone who can help
>> me with this problem:
>>
>> I define a function "optim.fun" and want as output the
>> sum of squared errors between predicted and measured
>> values, as follows:
>>
>> optim.fun <- function (ST04, SM08b, ch2no, a, b, d, E)
>> {
>> predR <-
>>
(a*SM08b^I(2)+b*SM08b+d)*exp(E*((1/(283.15-227.13))-(1/(ST04+273.15-227.13))))
>> abserr <- abs(ch2no-predR)
>> qnum <- quantile(abserr, probs=0.95, na.rm=T)
>>
>> is.na(abserr) <- (abserr > qnum)
>> errsq <- sum(abserr^2, na.rm=T)
>> errsq
>> }
>>
>> Then I want to optimize parameters a,b,d and E as to
>> minimize the function output with the following:
>>
>> optim.model<-nls(nulldat ~ optim.fun(ST04, SM08b,
>> ch2no, a, b, d, E), data=tower,
>> start=list(a=-0.003,b=0.13,d=0.50, E=400), na.action >>
na.exclude )
>>
>> were nulldat=0
>> At this point I get the following error message:
>>
>> Error in qr.default(.swts * attr(rhs, "gradient")) :
>> NA/NaN/Inf in foreign function call (arg 1)
>> Warning messages:
>> 1: In if (na.rm) x <- x[!is.na(x)] else if
>> (any(is.na(x))) stop("missing values and NaN's not
>> allowed if 'na.rm' is FALSE") ... :
>> the condition has length > 1 and only the first
>> element will be used
>> (this warning is repeated 12 times)
>>
>> Question: what does the error mean? What am I doing
>> wrong?
>> Thanks a bunch.
>> Nano
>> Jen, Germany
>> Max Planck for Biogeochemistry
>>
>>
>>
>>
>>
>>
____________________________________________________________________________________
>>
>> [[elided Yahoo spam]]
>>
>> ______________________________________________
>> R-help at 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.
>>
>
> ______________________________________________
> R-help at 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.
>
>
--
View this message in context:
http://www.nabble.com/function-in-nls-argument-tp17108100p17127514.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help at 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.
____________________________________________________________________________________
[[elided Yahoo spam]]
You indeed found a bug. I can reproduce it (which I should have tried to do on other examples in the first place!). Thanks for finding it. It will be fixed in version 1.1-0 which I will submit to CRAN soon. On Fri, 9 May 2008, elnano wrote:> > Find the data (data_nls.lm_moyano.txt) here: > ftp://ftp.bgc-jena.mpg.de/pub/outgoing/fmoyano > > > > Katharine Mullen wrote: > > > > Thanks for the details - it sounds like a bug. You can either send me the > > data in an email off-list or make it available on-line somewhere, so that > > I and other people can download it. > > > > > > ______________________________________________ > > R-help at 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. > > > > > > -- > View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17146812.html > Sent from the R help mailing list archive at Nabble.com. > > ______________________________________________ > R-help at 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. >