############################################################################
#						 		           #
#  Program to perform iteratively reweighted least squares by              #
#  by the simple Gauss-Newton approach (with no fancy modifications).      #
#								           #
#  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 -- if 2 successive 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                    #
#								           #
############################################################################

#  Indomethacin PK data with biexponential model

#  function to determine starting values

self.start <- function(y,x){

    n <- length(y)	

#  2nd phase -- start with last 3 observations

    isitout <- F
    i <- 2
    thisb34 <- 0

    while (!isitout){

#  set estimate to previous

         thisb34.old <- thisb34
	
	 k <- i+1
	
#  fit straight line on log scale by OLS
	
         thisb34 <- lsfit(x[(n-i):n],log(y[(n-i):n]))$coef
         res <- lsfit(x[(n-i):n],log(y[(n-i):n]))$residuals
         sigma2 <- sum(res^2)/(k-2)

#  form prediction interval at next x

         x0 <- x[n-i-1]
	 logy0 <- thisb34[1]+thisb34[2]*x0
	 xbar <- mean(x[(n-i):n])
         sxx <- (k-1)*var(x[(n-i):n])
         end <- sqrt(sigma2*(1+ 1/k + (x0-xbar)^2/sxx))
	 nextlogy <- log(y[n-i-1])
	 t95 <- qt(0.95,k-2)
         lo <- logy0 - t95*end
	 up <- logy0 + t95*end

#  is the next observation contained in the 90% prediction interval?
	
	 isitout <- (nextlogy<=lo)|(nextlogy>=up)
	
	 i <- i+1

    }

#  upon exit, thisb34.old contains the starting values and k-1 contains
#  the number of observations used in estimating the 2nd phase

#  now form residuals from this fit to estimate the 1st phase
	
    res <- y[1:(n-k+1)] - exp(thisb34.old[1]+thisb34.old[2]*x[1:(n-k+1)])
    thisb12 <- lsfit(x[1:(n-k+1)],log(res))$coef

#  the final starting values

    bstart <-c(exp(thisb12[1]),-thisb12[2],exp(thisb34.old[1]),-thisb34.old[2])

#  return the start value and number of observations used to estimate
#  2nd phase
	
    list(bstart=bstart,num=k-1)
}
     

#  Functions defining f, fb, g

#  Regression model 

f <- function(x,b){
       f1 <- b[1]*exp(-b[2]*x) + b[3]*exp(-b[4]*x)
       f1
}

fb <- function(x,b){
	n <- length(x)
	p <- length(b)
	fb1 <- matrix(0,n,p)
	fb1[,1] <- exp(-b[2]*x)
	fb1[,2] <- -b[1]*x*exp(-b[2]*x)
	fb1[,3] <- exp(-b[4]*x)
	fb1[,4] <- -b[3]*x*exp(-b[4]*x)
	fb1
}

g <- function(b,theta,z){

#  power variance function
	
	g1 <- exp(theta*log(f(z,b)))
	g1
}

#  data

thedat <- matrix(scan("biexp.dat"),ncol=2,byrow=T)

x <- thedat[,1]
y <- thedat[,2]

#  get starting value for beta by "self-starting" function and number
#  of observations used to estimate second phase

start.value <- self.start(y,x)	
bstart <- start.value$bstart
num <- start.value$num

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

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

#  name of output file
	
outfile <- "self.out"

cat("FIT OF BIEXPONENTIAL DATA BY IRWLS USING SELF-STARTING METHOD",file=outfile,"\n","\n",append=F)

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

#  output value of starting value and number of observations used

cat("Start value = ",round(bstart,6),file=outfile,"\n",append=T) 
cat("and # obs used in 2nd phase = ",num,file=outfile,"\n","\n",append=T)

#  generic code to perform irwls 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 IRWLS
	
irwls <- function(y,x,b,theta,tol,imax,f,fb,g){

#  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)
	  w <- 1/g(b,theta,x)^2

#  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)

}

#  loop through values of theta and compute beta, sigma by IRWLS
	
for (j in 1:length(tvec)){

	theta <- tvec[j]
        cat("Fit of model with power",round(theta,6),file=outfile,"\n",
	"\n",append=T)

	result <- irwls(y,x,bstart,theta,tol,imax,f,fb,g)
	betahat <- result$beta
	cat("Estimate of beta",round(betahat,6),file=outfile,"\n",append=T)
	cat("Number of iterations",result$iend,file=outfile,"\n",append=T)

#  compute final estimate of sigma^2, sigma

	w <- 1/g(betahat,theta,x)^2
	s2 <- sum(w*(y-f(x,betahat))^2)/(n-p)

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

}