[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)

        stop("Too few unique irradiance values")

    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)

        stop("Too few unique irradiance values")

    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
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Thu Oct 27 02:16:16 2005

This archive was generated by hypermail 2.1.8 : Thu 27 Oct 2005 - 03:32:22 EST