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

#  HCV load data

#  Functions defining f, fb, g

f <- function(x,b){
       
    v0 <- exp(b[1])
    cc <- exp(b[2])
    e <- exp(b[3])
    t0 <- 0.10      #  t0 is fixed
    tt <- (x>t0)*(x-t0)
    f1 <- v0*(1-e+e*exp(-cc*tt))*1e6
    f1
}

fb <- function(x,b){

    n <- length(x)
    p <- length(b)
    v0 <- exp(b[1])
    cc <- exp(b[2])
    e <- exp(b[3])
    t0 <- 0.10      #  t0 is fixed
    tt <- (x>t0)*(x-t0)
    f1 <- v0*(1-e+e*exp(-cc*tt))*1e6

    fb1 <- matrix(0,n,p)
    fb1[,1] <- f1
    fb1[,2] <- -v0*e*exp(-cc*tt)*cc*tt*1e6
    fb1[,3] <- v0*e*(exp(-cc*tt)-1)*1e6
    fb1
}

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

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

#  data

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

#  starting value for beta
	
bstart <- log(c(4.0,3.5,0.85))

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

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

#  name of output file
	
outfile <- "hcv.irwlsfit"

cat("FITS OF HCV LOAD DATA BY IRWLS FOR DIFFERENT THETA VALUES",
	file=outfile,"\n","\n",append=F)

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

#  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",theta,file=outfile,"\n","\n",append=T)

	result <- irwls(y,x,bstart,theta,tol,imax,f,fb,g)
	betahat <- result$beta
	cat("Estimate of beta",betahat,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",s2,file=outfile,"\n",append=T)
	cat("Estimate of sigma",sqrt(s2),file=outfile,"\n","\n","\n",append=T)

}