# R-beta: nonlinear least squares

Jim Lindsey (jlindsey@luc.ac.be)
Mon, 16 Mar 1998 10:34:28 +0100 (MET)

```From: Jim Lindsey <jlindsey@luc.ac.be>
Message-Id: <9803160934.AA29275@alpha.luc.ac.be>
Subject: R-beta: nonlinear least squares
To: r-help@stat.math.ethz.ch (r-help)
Date: Mon, 16 Mar 1998 10:34:28 +0100 (MET)

Because interest in nonlinear least squares has been expressed on this
list, I enclose a short function to do it. It is a stripped down
version of a function in my nonlinear library (I hope I didn't
introduce any errors in the process.) Jim

--------------------------------------------------------------------------
nls <- function(y, mu=NULL, p=NULL, wt=1, delta=1,
stepmax=10*sqrt(p%*%p), steptol=0.00001, iterlim=100, fscale=1){
call <- sys.call()
if(missing(p))stop("Initial parameter estimates must be supplied")
else np <- length(p)
if(!is.function(mu))stop("A mean function must be supplied")
fn <- function(p) sum(wt*(y-mu(p))^2)
if(fscale==1)fscale <- fn(p)
if(is.na(fscale))
stop("Non-numerical function value: probably invalid initial values")
z0 <- nlm(fn, p=p, hessian=T, print.level=print.level, typsiz=typsiz,
steptol=steptol, iterlim=iterlim, fscale=fscale)
if(length(delta)==1)delta <- rep(delta,length(y))
if(length(wt)==1)wt <- rep(wt,length(y))
maxlike <- length(y)*(log(2*pi*z0\$minimum/sum(wt))+1)/2-sum(log(delta))
fitted.values <-  as.vector(mu(z0\$estimate))
residuals <-  y-fitted.values
if(np==1)cov <- 1/z0\$hessian
else {
a <- qr(z0\$hessian)
if(a\$rank==np)cov <- solve(z0\$hessian)
else cov <- matrix(NA,ncol=np,nrow=np)}
cov <- 2*cov*z0\$minimum/sum(wt)
se <- sqrt(diag(cov))
z1 <- list(
call=call,
delta=delta,
mu=mu,
prior.weights=wt,
maxlike=maxlike,
fitted.values=fitted.values,
residuals=residuals,
aic=maxlike+np+1,
df=sum(wt)-np,
coefficients=z0\$estimate,
np=np,
se=se,
cov=cov,
corr=cov/(se%o%se),
iterations=z0\$iterations,
code=z0\$code)
class(z1) <- "nls"
return(z1)}

residuals.nls <- function(z) z\$residuals
fitted.values.nls <- function(z) z\$fitted.values
coefficients.nls <- function(z) z\$coefficients
weights.nls <- function(z) z\$prior.weights
df.residual.nls <- function(z) z\$df
deviance.nls <- function(z) 2*z\$maxlike

print.nls <- function(z) {
cat("\nCall:\n",deparse(z\$call),"\n\n",sep="")
if(z\$code>2)cat("Warning: no convergence - error",z\$code,"\n\n")
t <- deparse(z\$mu)
cat("Mean function:",t[2:length(t)],sep="\n")
cat("\n-Log likelihood   ",z\$maxlike,"\n")
cat("Degrees of freedom",z\$df,"\n")
cat("AIC               ",z\$aic,"\n")
cat("Iterations        ",z\$iterations,"\n\n")
cat("Mean parameters:\n")
coef.table <- cbind(z\$coefficients[1:z\$np], z\$se[1:z\$np])
dimnames(coef.table) <- list(seq(1,z\$np), c("estimate", "se"))
print.default(coef.table, digits=4, print.gap=2)