require(deSolve); # for ODE solver
require(stats4);   # for maximum likelihood estimation 

########################################################## 
# Data: Paramecium abundance with predators absent  
########################################################## 
tvals=1:14; 
xvals=c(15.58,30.04,66.05,141.6,274.6,410.0,468.8,526.4,
472.5, 496.6,489.5,492.0,496.8,473.0);

########################################################## 
# Define the model
##########################################################
ThetaLogist=function(t,x,parms) {
    r=parms[1]; K=100*parms[2]; theta=parms[3]
    xdot=r*x*(1-(x/K)^theta);	
    return(list(xdot));
}

########################################################## 
# Least squares objective function 
########################################################## 
TLfit=function(p) {
	parms=p[1:3]; x0=p[4]; times=1:14; 
	Nhat=lsoda(x0,times,ThetaLogist,parms)[,2];
	mse=mean((Nhat-xvals)^2);
	return(mse);
}

########################################################## 
# Minimize the objective function 
########################################################## 
p0=c(1,5,1,15); # initial guess at parameters 	
fit=optim(p0,TLfit); # call the optimizer 	
fit2=optim(fit$par,TLfit);  # re-call optimizer to check convergence 
fit$par; fit2$par; 

########################################################## 
# Maximum Likelihood fit 
########################################################## 
TLNegLogLik=function(r,K,theta,x0) {
	parms=c(r,K,theta); times=1:14; 
	Nhat=lsoda(x0,times,ThetaLogist,parms)[,2]
    -sum(dnorm(xvals-Nhat,mean=0,sd=0.42*sqrt(Nhat),log=TRUE));  
}
fitTL.MLE=mle(TLNegLogLik,start=list(r=1,K=5,theta=1,x0=15),method="Nelder-Mead")
			
########################################################## 
# Solve ODE with MLE parameters 
########################################################## 
p=coef(fitTL.MLE); parms=p[1:3]; x0=p[4]; times=(2:28)/2; 
Nhat=lsoda(x0,times,ThetaLogist,parms,hmin=0.001)[,2];

########################################################## 
#plot data and solution 
########################################################## 
win.graph(8,4); 
par(cex.axis=1.35,cex.lab=1.35,mar=c(5,5,1,1))
plot(times,Nhat,type="l",xlab="Time (days)",ylab="Population count",ylim=c(0,max(xvals))); 
points(tvals,xvals,type="p",pch=16,cex=1.15); 
StdDev=0.42*sqrt(Nhat);
matpoints(times,cbind(Nhat-2*StdDev,Nhat+2*StdDev),type="l",lty=2,col="black");