#!/bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #!/bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create the files:
# ./nlsfit.d
# ./source.S
# This archive created: Wed Nov 21 13:32:06 1990
export PATH; PATH=/bin:$PATH
if test -f './nlsfit.d'
then
echo shar: over-writing existing file "'./nlsfit.d'"
fi
cat << \SHAR_EOF > './nlsfit.d'
.BG
.FN nlsfit
.TL
nonlinear regression
.CS
nlsfit(model, theta=eval(model$theta), maxiter=50, minfactor=1/2^10,
tolerance=0.001, verbose=T)
.AG model
a function that takes arguments "theta" and "derivs" where "theta" is
the parameter vector for the model and "derivs" is a logical value that
indicates if the gradient matrix should be calculated. It "derivs" is
false, a vector of residuals should be returned. If "derivs" is true,
a structure with components "residuals" and "gradient" should be
returned. The "gradient" component is the gradient of the predicted
responses with respect to the parameter vector. Thus it is the
negative of the gradient of the residual vector.
.AG theta
Starting estimates for the parameters. If the model function has a
default value of the parameters, that is used as the default here.
.AG maxiter
Maximum number of iterations in the Gauss-Newton algorithm. Default is 50.
.AG minfactor
Minimum factor by which to multiply the Gauss-Newton increment when
trying to get a decrease in the sum of squares. Default is 1/2^10 or
approximately 0.001. When the stepfactor becomes less than this, it
usually indicates that the gradient has been incorrectly specified.
.AG tolerance
Tolerance level for the orthogonality convergence criterion. Defaults
to 0.001.
.AG verbose
Logical value indicating if verbose output on the iterations should be
written to the terminal. Defaults to T.
.RT
If the iterations are successful in locating a local minimum of the
sum of squares of the residuals, a list is returned with components
corresponding to those returned from the lsfit function. This means
the ls.summary function can be used to produce a summary based on the
linear approximation to the nonlinear regression model. The
components are:
.RC model
argument as passed in to nlsfit.
.RC coef
estimates of the parameters.
.RC qr
the QR decomposition of the gradient matrix.
.RC residuals
the residual vector at the parameter estimates.
.RC criterion
the value of the orthogonality convergence criterion at the optimum.
.PP
This function is a slightly modified copy of the function given in
Appendix A3.1.1 of Bates and Watts, "Nonlinear Regression Analysis and
Its Applications", Wiley, 1988. That book describes the Gauss-Newton
algorithm used here and the orthogonality convergence criterion.
.EX
> # Using the Puromycin example enclosed with this function
> Pur.out <- nlsfit(Puromycin,verb=F)
> Pur.out$coef
[1] 212.68294013 0.06412003
> Pur.summ <- ls.summary(Pur.out)
> Pur.summ$std.err
[,1]
[1,] 6.947111233
[2,] 0.008280816
> Pur.summ$corr
[,1] [,2]
[1,] 1.0000000 0.7650823
[2,] 0.7650823 1.0000000
.SA
lsfit, qr
.WR
SHAR_EOF
if test -f './source.S'
then
echo shar: over-writing existing file "'./source.S'"
fi
cat << \SHAR_EOF > './source.S'
"nlsfit"<-
function(model, theta = eval(as.list(model)$theta), maxiter = 50, minfactor = 1/2^10,
tolerance = 0.001, qrtol = 1e-7, verbose = T)
{
resid <- model(theta, derivs = F)
newssq <- sum(resid^2)
P <- length(theta)
N <- length(resid)
ndof <- N - P
mult <- sqrt(ndof/P)
iteration <- 0
stepfactor <- 1
repeat {
iteration <- iteration + 1
if(iteration > maxiter)
stop("Maximum number of iterations exceeded")
oldssq <- newssq
resgrad <- model(theta, derivs = T)
qrstr <- qr(resgrad$gradient, tol = qrtol)
if(qrstr$rank!=P)
stop("Singular gradient matrix")
incr <- qr.coef(qrstr, resgrad$residual)
qty <- qr.qty(qrstr, resgrad$residual)
converge <- mult * sqrt(sum(qty[1:P]^2)/sum(qty[ - (1:P)]^2))
if(verbose)
cat(iteration, "<", converge, ">", incr, fill = T)
if(converge < tolerance)
break
repeat {
trial <- theta + stepfactor * incr
if(stepfactor < minfactor)
stop("Step factor reduced below minimum")
newssq <- sum(model(trial, derivs = F)^2)
if(verbose)
cat(" ", stepfactor, "(", newssq, ")", trial,
fill = T)
if(newssq <= oldssq)
break
stepfactor <- stepfactor/2
}
stepfactor <- min(1, 2 * stepfactor)
theta <- trial
}
list(model = model, coef = theta, qr = qrstr, residuals = resgrad$resid,
criterion = converge)
}
"Puromycin"<-
function(theta = c(205, 0.08), derivs = F)
{
conc <- c(0.02, 0.02, 0.06, 0.06, 0.11, 0.11, 0.22, 0.22,
0.5600000000000001, 0.5600000000000001, 1.1, 1.1)
rate <- c(76, 47, 97, 107, 123, 139, 159, 152, 191, 201, 207, 200)
denom <- theta[2] + conc
res <- rate - ((theta[1] * conc)/denom)
if(!derivs)
return(res)
grad <- matrix(c(conc/denom, (( - theta[1]) * conc)/denom^2), ncol = 2)
list(residual = res, gradient = grad)
}
SHAR_EOF
# End of shell archive
exit 0