From: Katharine Mullen <kate_at_few.vu.nl>

Date: Fri, 09 May 2008 13:25:42 +0200 (CEST)

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. Received on Fri 09 May 2008 - 11:39:27 GMT

Date: Fri, 09 May 2008 13:25:42 +0200 (CEST)

On Fri, 9 May 2008, elnano wrote:

*>
*

> Thank you Katharine. I am certain nprint is affecting my solution. Let me

*> know how I can send the data (~300Kb). The script I used it:
**>
**> ST1 <- ST04
**> SM1 <- SM08
**> SR1 <- SRch2
**>
**> ST <- ST1 [!is.na(SR1)]
**> SM <- SM1 [!is.na(SR1)]
**> SR <- SR1 [!is.na(SR1)]
**> q <- 0.90
**>
**> p <- c("a"=-0.003, "b"=0.13, "c"=0.50, "E"=400)
**>
**> SRf <- 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, SR, SRfcall)
**> {
**> res <- SR - do.call("SRfcall", c(list(ST = ST, SM = SM), as.list(p)))
**> abserr <- abs(res)
**> qnum <- quantile(abserr, probs=q, na.rm=TRUE)
**> res[abserr > qnum]=0
**> res
**> }
**>
**> SRmodel<- nls.lm(par = p, fn = optim.f, SRfcall = SRf, ST = ST, SM = SM,
**> SR = SR, control = nls.lm.control(nprint=1))
**>
**> summary(SRmodel)
**> SRpred <- do.call(SRf, c(list(ST = ST1, SM = SM1), SRmodel$par))
**>
**> plot(SR1~SRpred)
**> a=c(0,2,4,6)
**> b=c(0,2,4,6)
**> lines(a,b,col=4)
**>
**> Changing the nprint argument to 0 (or removing nprint) results in different
**> and completely wrong solutions. The same is true for this equivalent
**> simplified script from your example.
**>
**> ST1 <- ST04
**> SM1 <- SM08
**> SR1 <- SRch2
**>
**> ST <- ST1 [!is.na(SR1)]
**> SM <- SM1 [!is.na(SR1)]
**> SR <- SR1 [!is.na(SR1)]
**> q <- 0.9
**>
**> p <- list(a=-0.003, b=0.13, c=0.50, E=400)
**>
**> optim.f <- function(p, ST, SM, SR, q) {
**> res <- SR -
**> (p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))
**> abserr <- abs(res)
**> qnum <- quantile(abserr, probs=q, na.rm=TRUE)
**> res[abserr > qnum] <- 0
**> res
**> }
**>
**> SRmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, SR = SR, q = q,
**> control = nls.lm.control(nprint=1))
**>
**> summary(SRmodel)
**> SRfun <- function(ST, SM, a, b, c, E)
**> {(a*SM^I(2)+b*SM+c)*exp(E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))))}
**> SRpred <- do.call(SRfun, c(list(ST = ST1, SM = SM1), SRmodel$par))
**>
**> plot(SR1~SRpred)
**> a=seq(0,6,1)
**> b=seq(0,6,1)
**> lines(a,b,col=4)
**>
**> Here, however, I get the following additional error after summary(SRmodel):
**> Error in param/se : non-numeric argument to binary operator
**> although the solution is same as for the first script.
**>
**> Note that a q of 95 is OK, a q of 90 is not, but a q of 50 is OK again...
**>
**>
**>
**> To make your example reproducible you have to provide the data somehow;
**> I am pretty sure nprint doesn't effect the solution, but if it does this
**> would be a bug and I would appreciate a reproducible report.
**>
**> The example in nls.lm is a little complicated in order to show how to use
**> an analytical expression for the gradient (possible to provide in the
**> argument jac); since you don't need this, note that your residual function
**> can be simplified into something along the lines of
**>
**> p <- list(a=-0.003, b=0.13, c=0.50, E=400)
**>
**> optim.f <- function(p, ST, SM, R, q) {
**> res <- R
**> -(p$a*SM^I(2)+p$b*SM+p$c)*exp(p$E*((1/(283.15-227.13))-(1/(ST+273.15-227.13))
**> abserr <- abs(res)
**> qnum <- quantile(abserr, probs=q, na.rm=T)
**> res[abserr > qnum] <- 0
**> res
**> }
**>
**> Rmodel <- nls.lm(par = p, fn = optim.f, ST = ST, SM = SM, R = R, q = q)
**>
**> On Thu, 8 May 2008, Fernando Moyano wrote:
**>
**> > 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.
**>
**> --
**> View this message in context: http://www.nabble.com/function-in-nls-argument-tp17108100p17145801.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.
**>
*

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. Received on Fri 09 May 2008 - 11:39:27 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 09 May 2008 - 13:30:35 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.
*