#########################################################################
#                                                                       #
#  Program to implement normal theory ML using profiling and            #
#  the "trick" to estimate beta and theta (and sigma)                   #
#                                                                       #
#  Details on the nls() function may be found in Chapter 10             #
#  of the book "Statistical Models in S" edited by J.M. Chambers        #
#  and T. Hastie (1993, Chapman and Hall)                               #
#                                                                       #
#  Variance proportional to a power 2*theta of the mean                 #
#                                                                       #
#########################################################################

#########################################################################
#                                                                       #
#  Define the mean function f                                           #
#                                                                       #
#########################################################################

meanfunc <- function(t,W0,W1,beta0,beta1,tau0,tau1,v){

     log2 <- log(2)
     t0 <- 20

     f10 <- W0*(1-exp(-log2*(t/tau0)^beta0))
     joinpt <- W0*(1-exp(-log2*(t0/tau0)^beta0))
     f11 <- (t>t0)*(joinpt + W1*(1-exp(-log2*((t-t0)*(t>t0)/tau1)^beta1)))
     f1 <-  (1-v)*f10 + v*( (t<=t0)*f10 + (t>t0)* f11 ) 

#  I have removed the derivatives, because R will get confused.  The 
#  mean function is part of the "trick function."  When the "fake 
#  regression problem is run, the nls function needs the derivative of
#  the trick function.  "meanfunc" will carry around the derivatives
#  below as an attribute, and this will cause nls to screw up computation
#  of derivatives for "trkfunc," as it will attempt to make use of these
#  derivatives in the wrong way.  
        
   f1
}

#########################################################################
#                                                                       #
#  Define the variance function g                                       #
#                                                                       #
#########################################################################
  
varfunc <- function(x,b1,b2,b3,b4,b5,b6,v,theta){
     g1 <- meanfunc(x,b1,b2,b3,b4,b5,b6,v)^theta
     g1
}

#########################################################################
#                                                                       #
#  Define the mean function for the "trick"                             #
#                                                                       #
#########################################################################
  
trkfunc <- function(y,x,b1,b2,b3,b4,b5,b6,v,theta){

#  define the geometric mean of the variance function

   gdot <- exp( mean( log(varfunc(x,b1,b2,b3,b4,b5,b6,v,theta)) ) )
   trk <- (y-meanfunc(x,b1,b2,b3,b4,b5,b6,v))*gdot/varfunc(x,b1,b2,b3,b4,b5,b6,v,theta)	

#  too lazy to find analytic derivatives -- let nls compute them numerically!

   trk
}

#########################################################################
#                                                                       #
#  The data                                                             #
#                                                                       #
#########################################################################

thedata <- scan("dissolution.dat")
thedata <- matrix(thedata,ncol=3,byrow=T)

outfile <- "diss_ml.Rout"

v <- thedata[,1]
x <- thedata[,2]
y <- thedata[,3]
n <- length(x)

p <- 6

#########################################################################
#                                                                       #
#  Define list of starting values for beta and theta                    #
#                                                                       #
#########################################################################

allstart <- list(b1=60, b2=30, b3=1, b4=1, b5=10, b6=10, theta=1)

#########################################################################
#                                                                       #
#  Create the data frame for nls() with "dummy" responses               #
#                                                                       #
#########################################################################

dummy <- rep(0,n)
thedat <- data.frame(dummy,x,y,v)

#########################################################################
#                                                                       #
#  Call nls() with the trick function and dummy data                    #
#                                                                       #
#########################################################################
	
mlefit <- nls(dummy ~ trkfunc(y,x,b1,b2,b3,b4,b5,b6,v,theta),thedat,allstart,trace=T)
	
#########################################################################
#                                                                       #
#  Extract the estimates from the object mlefit and compute             #
#  sigma                                                                #
#                                                                       #
#########################################################################

bt <- coef(mlefit)
bmle <- bt[1:6]
theta <- bt[7]
mu <- meanfunc(x,bmle[1],bmle[2],bmle[3],bmle[4],bmle[5],bmle[6],v)
resid <- y-mu
sigmle <- sqrt(sum((resid/varfunc(x,bmle[1],bmle[2],bmle[3],bmle[4],bmle[5],bmle[6],v,theta))^2)/n)
	
#########################################################################
#                                                                       #
#  Print out the results to a file -- we round the results to 6         #
#  decimal places to the right of the decimal                           #
#                                                                       #
#########################################################################

cat("FIT OF THE DISSOLUTION DATA BY NORMAL ML USING THE TRICK",
	file=outfile,"\n","\n","\n",append=F)
	
cat("beta ",round(bmle,6),file=outfile,"\n","\n","\n",append=T)
cat("theta ",round(theta,6),file=outfile,"\n","\n","\n",append=T)
cat("sigma ",round(sigmle,6),file=outfile,"\n","\n","\n",append=T)