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,
	print.level=0, typsiz=abs(p), ndigit=10, gradtol=0.00001,
	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,
		ndigit=ndigit, gradtol=gradtol, stepmax=stepmax,
		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,
		rss=z0$minimum,
		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),
		gradient=z0$gradient,
		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)
	cat("\nVariance estimate:",z$rss/z$df,"\n")
	if(z$np>1){
		cat("\nCorrelations:\n")
		dimnames(z$corr) <- list(seq(1,z$np),seq(1,z$np))
		print.default(z$corr, digits=4)}
	invisible(z)}
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._