model2grfun: Generate a gradient function from a nonlinear model...

In nlmrt: Functions for Nonlinear Least Squares Solutions

Description

Given a nonlinear model expressed as an expression of the form lhs ~ formula_for_rhs and a start vector where parameters used in the model formula are named, attempts to build the the R function for the gradient of the residual sum of squares. As a side effect, a text file with the program code is generated.

Usage

   model2grfun(modelformula, pvec, funname="mygr", filename=NULL)

Arguments

modelformula	This is a modeling formula of the form (as in nls) lhsvar ~ rhsexpression for example, y ~ b1/(1+b2*exp(-b3*tt)) You may also give this as a string.
pvec	A named parameter vector. For our example, we could use start=c(b1=1, b2=2.345, b3=0.123) WARNING: the parameters in the output function will be used in the order presented in this vector. Names are NOT respected in the output function.
funname	The (optional) name for the function that is generated in the file named in the next argument. The default name is 'mygr'.
filename	The (optional) name of a file that is written containing the (text) program code for the function. If NULL, no file is written.

Details

None.

Value

An R function object that computes the gradient of the sum of squared residuals of a nonlinear model at a set of parameters.

Author(s)

John C Nash <[e-mail]>

References

Nash, J. C. (1979, 1990) _Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation._ Adam Hilger./Institute of Physics Publications

See Also

Function nls(), packages optim and optimx.

Examples

  cat("See also examples in nlmrt-package.Rd\n")
  require(numDeriv)
  y <- c(5.308, 7.24, 9.638, 12.866, 17.069, 23.192, 31.443, 38.558, 
      50.156, 62.948, 75.995, 91.972)  # for testing
  tt <- seq_along(y)  # for testing
  f <- y ~ b1/(1 + b2 * exp(-1 * b3 * tt))
  p <- c(b1 = 1, b2 = 1, b3 = 1)
  mygr <- model2grfun(f, p)
  myss <- model2ssfun(f, p) # for check
  cat("mygr:\n")
  print(mygr)
  
  ans <- mygr(p, tt = tt, y = y)
  print(ans)
  gnum <- grad(myss, p, tt = tt, y = y)
  cat("Max(abs(ans-gnum)) = ",max(abs(ans-gnum)),"\n")

  bnew <- c(b1 = 200, b2 = 50, b3 = 0.3)
  ans <- mygr(prm = bnew, tt = tt, y = y)
  print(ans)
  gnum <- grad(myss, bnew, tt = tt, y = y)
  cat("Max(abs(ans-gnum)) = ",max(abs(ans-gnum)),"\n")
  
  cat("Test with un-named vector\n")  # At 20120424 should fail
  bthree <- c(100, 40, 0.1)
  ans <- mygr(prm = bthree, tt = tt, y = y)
  print(ans)
  gnum <- grad(myss, bthree, tt = tt, y = y)
  cat("Max(abs(ans-gnum)) = ",max(abs(ans-gnum)),"\n")

Example output

See also examples in nlmrt-package.Rd
Loading required package: numDeriv
mygr:
function (prm, y = NULL, tt = NULL) 
{
    b1 <- prm[[1]]
    b2 <- prm[[2]]
    b3 <- prm[[3]]
    localdf <- data.frame(y, tt)
    jstruc <- with(localdf, eval({
        .expr4 <- exp(-1 * b3 * tt)
        .expr6 <- 1 + b2 * .expr4
        .expr11 <- .expr6^2
        .value <- b1/.expr6 - y
        .grad <- array(0, c(length(.value), 3), list(NULL, c("b1", 
            "b2", "b3")))
        .grad[, "b1"] <- 1/.expr6
        .grad[, "b2"] <- -(b1 * .expr4/.expr11)
        .grad[, "b3"] <- b1 * (b2 * (.expr4 * tt))/.expr11
        attr(.value, "gradient") <- .grad
        .value
    }))
    jacmat <- attr(jstruc, "gradient")
    resids <- as.numeric(eval(b1/(1 + b2 * exp(-1 * b3 * tt)) - 
        y))
    grj <- as.vector(2 * crossprod(jacmat, resids))
}
<environment: 0x561ea8b9c228>
[1] -824.042084    4.764888  -11.025384
Max(abs(ans-gnum)) =  2.85223e-07 
[1]    -17.8438     48.2491 -24559.8187
Max(abs(ans-gnum)) =  2.885245e-07 
Test with un-named vector
[1]    -45.1255    105.5918 -41527.5337
Max(abs(ans-gnum)) =  9.532741e-07 

nlmrt documentation built on May 1, 2019, 11:31 p.m.