############################################################################
#						 		           #
#  Program to perform the 3-step GLS algorithm:                            #
#						 		           #
#       (i) Get OLS estimate                                               #
#      (ii) Form estimated weights                                         #
#     (iii) Do WLS with weights from (ii)                                  #
#						 		           #
#  The user must define 3 functions below:                                 #
#								           #
#	f(x,b)	the nonlinear regression model depending on a (n x k)      #
#		matrix of predictor variables and a (p x 1) regression     #
#		parameter b                                                #
#								           #
#	fb(x,b)	the (n x p) matrix of derivatives of f wrt b               #
#								           #
#	g(b,theta,z)  the variance model with variance parameters theta    #
#								           #
#  The user must also provide the following:                               #
#								           #
#	y	(n x 1) data vector                                        #
#								           #
#	x	(n x k) precictor values                                   #
#								           #
#	bstart	starting value for b                                       #
#								           #
#	theta	the fixed, known value for theta                           #
#								           #
#	tol	convergence criterion for WLS --if 2 iterates of b         #
#		are within tol of each other relative to the current       #
#		value, declare that the algorithm has converged            #
#								           #
#	imax	maximum number of iterations to perform in WLS             #
#								           #
#       cmax    number of GLS iterations to perform                        #
#								           #
############################################################################


#  Functions defining f, fb, g

f <- function(x,b){
        dose <- 5.50
        ka <- exp(b[1])
        cl <- exp(b[2])
         v <- exp(b[3])
        ke <- cl/v
        kk <- ka-ke
        f1 <- dose*ka*(exp(-ke*x)-exp(-ka*x))/(v*kk)
	f1
}

fb <- function(x,b){
	n <- length(x)
	p <- length(b)
        dose <- 5.50
        ka <- exp(b[1])
        cl <- exp(b[2])
         v <- exp(b[3])
        ke <- cl/v
        kk <- ka-ke
        f1 <- dose*ka*(exp(-ke*x)-exp(-ka*x))/(v*kk)
        tt <- dose*x*ka/(v*kk)

	fb1 <- matrix(0,n,p)
	fb1[,1] <- -f1*(ke/kk)+tt*ka*exp(-ka*x)
	fb1[,2] <- f1*(ke/kk)-tt*ke*exp(-ke*x)
	fb1[,3] <- -f1*(ka/kk)+tt*ke*exp(-ke*x)
	fb1
}
	
#  power variance function

g <- function(b,theta,z){	
	g1 <- exp(theta*log(f(z,b)))
	g1
}

#  data
	
thedat <- matrix(scan("theo10.dat"),ncol=2,byrow=T)
x <- thedat[,1]
y <- thedat[,2]

#  starting value for beta
	
bstart <- log(c(0.75,0.05,0.50))

#  values of theta to consider
	
tvec <- c(0.0,0.5,1.0)

#  convergence tolerance and max number of iterations for WLS routine
	
tol <- 1e-8
imax <- 500

#  number of iterations of GLS algorithm to perform

cmax <- 20

#  name of output file
	
outfile <- "theo10.glsfit"

cat("FITS OF SUBJECT 10 THEOPHYLLINE DATA BY 3-STEP GLS FOR DIFFERENT THETA VALUES", file=outfile,"\n","\n",append=F)

########################################################################

#  generic code to perform wls and compute sigma

n <- length(x)
p <- length(bstart)

#  function specifying convergence criterion
	
converged <- function(tol,db,iter,imax){
	bmax <- max(abs(db))
	gmax <- iter==imax
	bmax <- bmax<tol
	converged <- gmax | bmax
	converged
}

#  function to perform WLS with fixed weights
	
wls <- function(y,x,b,tol,imax,f,fb,w){

#  initialize

	db <- 1
	iter <- 1

#  start iterating - stop when converged
	
	while(! converged(tol,db,iter,imax)){

#  current gradient matrix and weights
	
	  Xb <- fb(x,b)

#  construct X'WX
	
	  XbWXb <- t(Xb) %*% diag(w) %*% Xb

#  constructed vaariable z

	  z <- Xb %*% b + (y-f(x,b))
	  XbWz <- t(Xb) %*% diag(w) %*% z

#  get new iterate

          bnew <- solve(XbWXb) %*% XbWz

#  compute relative change

	  db <- (bnew-b)/b

	  b <- bnew
	  iter <- iter+1
  }


#  list of results returned by function
	
	list(beta=b,iend=iter-1)

}

#  initial fit by OLS -- call wls function with w vector of all ones

cat("Initial fit of model by OLS",file=outfile,"\n","\n",append=T)

w <- rep(1,n)

result <- wls(y,x,bstart,tol,imax,f,fb,w)

bols <- result$beta
cat("OLS estimate of beta",bols,file=outfile,"\n",append=T)

#  estimates of sigma^2 and sigma using OLS
	
s2 <- sum((y-f(x,bols))^2)/(n-p)

cat("\n","\n",file=outfile,append=T)
cat("Estimate of sigma squared",s2,file=outfile,"\n",append=T)
cat("Estimate of sigma",sqrt(s2),file=outfile,"\n","\n","\n",append=T)

#  for each non-zero value of theta, iterate steps (ii) and (iii) cmax times
	
for (j in 2:length(tvec)){

	theta <- tvec[j]
        cat("3-step GLS fit of model with power",theta,file=outfile,
	"\n","\n",append=T)

	bgls <- bols

#  GLS iteration loop
	
	for (k in 1:cmax){       

#  step (ii) -- form weights

	  w <- 1/g(bgls,theta,x)^2
	
#  step (iii) -- WLS with these estimated weights	
	
          result <- wls(y,x,bstart,tol,imax,f,fb,w)
	  bgls <- result$beta
	  cat("Estimate of beta at",k,"iterations",bgls,file=outfile,
	  "\n",append=T)

	}

#  compute final estimates of sigma^2, sigma

	w <- 1/g(bgls,theta,x)^2

	s2 <- sum(w*(y-f(x,bgls))^2)/(n-p)

        cat("\n","\n",file=outfile,append=T)
	cat("Estimate of sigma squared",s2,file=outfile,"\n",append=T)
	cat("Estimate of sigma",sqrt(s2),file=outfile,"\n","\n","\n",append=T)

}