Fitting a 4-parameter logistic curve using the Levenburg-Marquardt algorithm

lma.logfit.R

require('minpack.lm');
###### START FUNCTION DEFINITIONS ###########

##Coded by M. Furia June 2007

lma.dologfit = function(x,y){
	compute.logistic.x = function(y,min,max,infl,slope){
        	return(infl - (1/slope)*log10((max - min)/(y - min)- 1));
	}
	# define the logistic function
	f = function(x,min,max,infl,slope){
		expr = expression(min + (max-min)/(1 + 10^((infl-x)*slope)));
		return(eval(expr));
	}
	j = function(x,min,max,infl,slope){
		expr = expression(min + (max-min)/(1 + 10^((infl-x)*slope)));
		retVal = c(eval(D(expr,'min')),
			eval(D(expr,'max')),
			eval(D(expr,'infl')),
			eval(D(expr,'slope'))
		)
		return(retVal);
	}
	# the function for the LMA algorithm to minimize
	fcn = function(p,y,x,fcall,jcall){
		y - do.call("fcall",c(list(x=x),as.list(p)));
	}
	# Jacobian
	fcn.jac = function(p,y,x,fcall,jcall){	
		-do.call("jcall",c(list(x=x),as.list(p)));
	}
	# compute the initial guesses
	guess= c('min' = min(y), 'max'=max(y), 'infl'=median(x), 'slope'=1);
	
	#control = c('maxfev' = 10000);
	# do the fit
	out = nls.lm(par=as.list(guess),fn = fcn, jac =fcn.jac,fcall=f, jcall=j, x=x,y=y)#,control=control);

	# compute the x and y residuals and predicted values	

	out$pred.y = f(x,out$par$min, out$par$max, out$par$infl, out$par$slope);
	out$resid.y = out$pred.y - y;
	sse = sum(out$resid.y^2);
	sst = sum((y - mean(y))^2);
	out$rSquare.y = 1 - (sse/sst);

	##out$pred.x = compute.logistic.x(y,out$par$min, out$par$max, out$par$infl, out$par$slope)
	##out$resid.x = out$pred.x - x;

	return(out);
}

lma.dolinfit = function(x,y){
	compute.linear.x = function(y,slope,intercept){
        	return((y-intercept)/slope);
	}
	# define the linear function
	f = function(x,slope,intercept){
		expr = expression(x*slope + intercept);
		return(eval(expr));
	}
	j = function(x,slope,intercept){
		expr = expression(x*slope + intercept);
		retVal = c(eval(D(expr,'slope')),
			eval(D(expr,'intercept'))
		)
		return(retVal);
	}
	# the function for the LMA algorithm to minimize
	fcn = function(p,y,x,fcall,jcall){
		y - do.call("fcall",c(list(x=x),as.list(p)));
	}
	# Jacobian
	fcn.jac = function(p,y,x,fcall,jcall){	
		-do.call("jcall",c(list(x=x),as.list(p)));
	}
	# compute the initial guesses
	guess= c('intercept'=0, 'slope'=1);

	#control = c('maxfev' = 10000);
	# do the fit
	out = nls.lm(par=as.list(guess),fn = fcn,fcall=f, jcall=j, x=x,y=y)#, control=control);

	out$pred.y = f(x,out$par$slope, out$par$intercept);
	out$resid.y = out$pred.y - y;
	

	# compute the x and y residuals and predicted values	
	out$pred.x = compute.linear.x(y,out$par$slope, out$par$intercept);
	out$resid.x = out$pred.x - x;

	out$expr = out$par$min + (out$par$max-out$par$min)/(1 + 10^((out$par$infl-x)*out$par$slope))

	return(out);
}

f.log = function(x,min,max,infl,slope){
	expr = expression(min + (max-min)/(1 + 10^((infl-x)*slope)));
	return(eval(expr));
}
f.lin = function(x,slope,intercept){
	expr = expression(x*slope + intercept);
	return(eval(expr));
}
compute.logistic.x = function(y,min,max,infl,slope){
       	return(infl - (1/slope)*log10((max - min)/(y - min)- 1));
}
compute.linear.x = function(y,slope,intercept){
       	return((y-intercept)/slope);
}