Description Details Author(s) References See Also Examples
The package provides some tools related to using the Nash variant of Marquardt's algorithm for nonlinear least squares.
Package: | nlsr |
Type: | Package |
Version: | 1.0 |
Date: | 2012-03-05 |
License: | GPL-2 |
This package includes methods for solving nonlinear least squares problems specified by a modeling expression and given a starting vector of named paramters. Note: You must provide an expression of the form lhs ~ rhsexpression so that the residual expression rhsexpression - lhs can be computed. The expression can be enclosed in quotes, and this seems to give fewer difficulties with R. Data variables must already be defined, either within the parent environment or else in the dot-arguments. Other symbolic elements in the modeling expression must be standard functions or else parameters that are named in the start vector.
The main functions in nlsr
are:
Nash variant of the Marquardt procedure for nonlinear least squares,
with bounds constraints, using a residual and optionally Jacobian
described as R
functions.
Nash variant of the Marquardt procedure for nonlinear least squares,
with bounds constraints, using an expression to describe the residual via
an R
modeling expression. The Jacobian is computed via symbolic
differentiation.
Uses nlxb
to solve nonlinear least squares then calls nls()
to create an object of type nls.
returns a function with header function(prm)
, which
evaluates the residuals (and if jacobian is TRUE the
Jacobian matrix) of the model at prm
. The residuals are defined
to be the right hand side of modelformula
minus the left hand side.
returns a function with header function(prm)
, which
evaluates the sum of squared residuals (and if gradient is TRUE
the
gradient vector) of the model at prm
.
returns the expression used to calculate the vector of residuals (and possibly the Jacobian) used in the previous functions.
John C Nash and Duncan Murdoch
Maintainer: <nashjc@uottawa.ca>
Nash, J. C. (1979, 1990) _Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation._ Adam Hilger./Institute of Physics Publications
Nash, J. C. (2014) _Nonlinear Parameter Optimization Using R Tools._ Wiley
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 | rm(list=ls())
# require(nlsr)
traceval <- TRUE # traceval set TRUE to debug or give full history
# Data for Hobbs problem
ydat <- 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
tdat <- seq_along(ydat) # for testing
# A simple starting vector -- must have named parameters for nlxb, nls, wrapnlsr.
start1 <- c(b1=1, b2=1, b3=1)
startf1 <- c(b1=1, b2=1, b3=.1)
eunsc <- y ~ b1/(1+b2*exp(-b3*tt))
cat("LOCAL DATA IN DATA FRAMES\n")
weeddata1 <- data.frame(y=ydat, tt=tdat)
weeddata2 <- data.frame(y=1.5*ydat, tt=tdat)
anlxb1 <- try(nlxb(eunsc, start=start1, trace=traceval, data=weeddata1))
print(anlxb1)
# illustrate predict
newdta <- colMeans(weeddata1)
newdta["tt"]<-25 # This only works for 1D example -- CAUTION
predn <- predict(anlxb1, as.list(newdta))
print(predn)
anlxb2 <- try(nlxb(eunsc, start=start1, trace=traceval, data=weeddata2))
print(anlxb2)
escal <- y ~ 100*b1/(1+10*b2*exp(-0.1*b3*tt))
suneasy <- c(b1=200, b2=50, b3=0.3)
ssceasy <- c(b1=2, b2=5, b3=3)
st1scal <- c(b1=100, b2=10, b3=0.1)
shobbs.res <- function(x){ # scaled Hobbs weeds problem -- residual
# This variant uses looping
if(length(x) != 3) stop("hobbs.res -- parameter vector n!=3")
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)
tt <- 1:12
res <- 100.0*x[1]/(1+x[2]*10.*exp(-0.1*x[3]*tt)) - y
}
shobbs.jac <- function(x) { # scaled Hobbs weeds problem -- Jacobian
jj <- matrix(0.0, 12, 3)
tt <- 1:12
yy <- exp(-0.1*x[3]*tt)
zz <- 100.0/(1+10.*x[2]*yy)
jj[tt,1] <- zz
jj[tt,2] <- -0.1*x[1]*zz*zz*yy
jj[tt,3] <- 0.01*x[1]*zz*zz*yy*x[2]*tt
attr(jj,"gradient")<-jj # needed for nlfb to function seamlessly with nlxb
return(jj)
}
cat("try nlfb\n")
st <- c(b1=1, b2=1, b3=1)
low <- -Inf
up <- Inf
## Not run:
ans1 <- nlfb(st, shobbs.res, shobbs.jac, trace=traceval)
ans1
cat("No jacobian function -- use internal approximation\n")
ans1n <- nlfb(st, shobbs.res, trace=TRUE, control=list(watch=TRUE)) # NO jacfn
ans1n
# tmp <- readline("Try with bounds at 2")
time2 <- system.time(ans2 <- nlfb(st, shobbs.res, shobbs.jac, upper=c(2,2,2),
trace=traceval))
ans2
time2
## End(Not run) # end dontrun
cat("BOUNDS")
st2s <- c(b1=1, b2=1, b3=1)
## Not run:
an1qb1 <- try(nlxb(escal, start=st2s, trace=traceval, data=weeddata1,
lower=c(0,0,0), upper=c(2, 6, 3), control=list(watch=FALSE)))
print(an1qb1)
tmp <- readline("next")
ans2 <- nlfb(st2s,shobbs.res, shobbs.jac, lower=c(0,0,0), upper=c(2, 6, 3),
trace=traceval, control=list(watch=FALSE))
print(ans2)
cat("BUT ... nls() seems to do better from the TRACE information\n")
anlsb <- nls(escal, start=st2s, trace=traceval, data=weeddata1, lower=c(0,0,0),
upper=c(2,6,3), algorithm='port')
cat("However, let us check the answer\n")
print(anlsb)
cat("BUT...crossprod(resid(anlsb))=",crossprod(resid(anlsb)),"\n")
## End(Not run) # end dontrun
tmp <- readline("next")
cat("Try wrapnlsr\n")
traceval <- TRUE
# Data for Hobbs problem
ydat <- 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
tdat <- seq_along(ydat) # for testing
start1 <- c(b1=1, b2=1, b3=1)
escal <- y ~ 100*b1/(1+10*b2*exp(-0.1*b3*tt))
up1 <- c(2,6,3)
up2 <- c(1, 5, 9)
weeddata1 <- data.frame(y=ydat, tt=tdat)
## Not run:
an1w <- try(wrapnlsr(escal, start=start1, trace=traceval, data=weeddata1))
print(an1w)
cat("BOUNDED wrapnlsr\n")
an1wb <- try(wrapnlsr(escal, start=start1, trace=traceval, data=weeddata1, upper=up1))
print(an1wb)
cat("BOUNDED wrapnlsr\n")
an2wb <- try(wrapnlsr(escal, start=start1, trace=traceval, data=weeddata1, upper=up2))
print(an2wb)
cat("TRY MASKS ONLY\n")
an1xm3 <- try(nlxb(escal, start1, trace=traceval, data=weeddata1,
masked=c("b3")))
printsum(an1xm3)
#an1fm3 <- try(nlfb(start1, shobbs.res, shobbs.jac, trace=traceval,
# data=weeddata1, maskidx=c(3)))
an1fm3 <- try(nlfb(start1, shobbs.res, shobbs.jac, trace=traceval,
data=weeddata1, maskidx=c(3)))printsum(an1fm3)
an1xm1 <- try(nlxb
(escal, start1, trace=traceval, data=weeddata1,
masked=c("b1")))
printsum(an1xm1)
#an1fm1 <- try(nlfb(start1, shobbs.res, shobbs.jac, trace=traceval,
an1fm1 <- try(nlfb(start1, shobbs.res, shobbs.jac, trace=traceval,
data=weeddata1, maskidx=c(1)))
printsum(an1fm1)
## End(Not run) # end dontrun
# Need to check when all parameters masked.
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.