# [R] self starting function for nonlinear least squares.

From: Bill Shipley <bill.shipley_at_usherbrooke.ca>
Date: Thu 27 Oct 2005 - 02:10:47 EST

Following on my posting of this morning, concerning a problem that I am having constructing a self-starting function for use with nls (and eventually with nlsList and nlme), the following is the self-starting function called NRhyperbola:

> NRhyperbola

function (Irr,theta,Am,alpha,Rd)

{

# Am is the maximum gross photosynthetic rate

# Rd is the dark resiration rate (positive value)

# theta is the shape parameter

# alpha is the initial quantum yeild

#

(1/(2*theta))*

(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2 -4*alpha*theta*Am*Irr))- Rd

}

attr(,"initial")

function (mCall, LHS, data)

{

xy <- sortedXyData(mCall[["Irr"]], LHS, data)

if (nrow(xy) < 3)

fit <- smooth.spline(xy[, "x"], xy[, "y"])

alpha <- predict(fit, x = 0, deriv = 1)\$y

if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)- 0)

if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0

# to correct for minimal values greater that zero

Rd <- (predict(fit, 0)\$y - delta.for.zero)

# to correct for maximal values

alpha <- predict(fit, x = max(xy[,"x"],na.rm=T), deriv = 1)\$y

if(alpha>0){

Am <- predict(fit,
max(xy[,"x"],na.rm=T))\$y+alpha*(3000-max(xy[,"x"],na.rm=T))

Am <- Am - Rd

}

theta1 <- rep(NA, 4)

light <- c(100, 200, 500, 800)

for (i in 1:4) {

y.pred <- predict(fit, light[i])\$y

theta1[i] <- ((y.pred + Rd) * (alpha * light[i] + Am) -

alpha * Am * light[i])/(y.pred + Rd)^2

}

theta <- mean(theta1)

Rd<- -1*Rd

value <- c(theta, Am, alpha, Rd)

names(value) <- mCall[c("theta", "Am", "alpha", "Rd")]

value

}

attr(,"class")

[1] "selfStart"

This function was created, following pp. 345-346 of Pinheiro & Bates (Mixed-effects models in S and S-PLUS) by first creating NRhyperbola:

NRhyperbola<- function (Irr,theta,Am,alpha,Rd)

{

# Am is the maximum gross photosynthetic rate

# Rd is the dark resiration rate (positive value)

# theta is the shape parameter

# alpha is the initial quantum yeild

#

(1/(2*theta))*

(alpha*Irr + Am - sqrt((alpha*Irr + Am)^2 -4*alpha*theta*Am*Irr))- Rd

}

then creating a function

NRhyperbolaInit<-

function (mCall, LHS, data)

{

xy <- sortedXyData(mCall[["Irr"]], LHS, data)

if (nrow(xy) < 3)

fit <- smooth.spline(xy[, "x"], xy[, "y"])

alpha <- predict(fit, x = 0, deriv = 1)\$y

if(min(xy[,"x"],na.rm=T)>0)delta.for.zero<-alpha*(min(xy[,"x"],na.rm=T)- 0)

if(min(xy[,"x"],na.rm=T)<=0)delta.for.zero<-0

# to correct for minimal values greater that zero

Rd <- (predict(fit, 0)\$y - delta.for.zero)

# to correct for maximal values

alpha <- predict(fit, x = max(xy[,"x"],na.rm=T), deriv = 1)\$y

if(alpha>0){

Am <- predict(fit,
max(xy[,"x"],na.rm=T))\$y+alpha*(3000-max(xy[,"x"],na.rm=T))

Am <- Am - Rd

}

theta1 <- rep(NA, 4)

light <- c(100, 200, 500, 800)

for (i in 1:4) {

y.pred <- predict(fit, light[i])\$y

theta1[i] <- ((y.pred + Rd) * (alpha * light[i] + Am) -

alpha * Am * light[i])/(y.pred + Rd)^2

}

theta <- mean(theta1)

Rd<- -1*Rd

value <- c(theta, Am, alpha, Rd)

names(value) <- mCall[c("theta", "Am", "alpha", "Rd")]

value

}

and then using "selfStart"

NRhyperbola<-selfStart(NRhyperbola,initial=NRhyperbolaInit)

Bill Shipley

[[alternative HTML version deleted]]

R-help@stat.math.ethz.ch mailing list