t <- c(0.166,0.333,0.5,0.666,1,1.5,2,2.5,3,4,5,6,8,10,12,24,30,48)
conc <- c(10.1,14.8,19.9,22.1,20.8,20.3,19.7,18.9,17.3,16.1,15,14.2,
13.2,12.3,10.8,6.5,4.6,1.7)

plot(t,conc)

plot(t,conc,log="x")

resid <- function(t,y,theta){
  yhat <- theta[3]*(exp(-1*theta[1]*t)-exp(-1*theta[2]*t))
  return(y-yhat)
}

a.matrix <- function(t,theta){
  a1 <- -1*theta[3]*t*exp(-1*theta[1]*t)
  a2 <- theta[3]*t*exp(-1*theta[2]*t)
  a3 <- exp(-1*theta[1]*t)-exp(-1*theta[2]*t)
  return(matrix(cbind(a1,a2,a3),ncol=3))
}

theta0 <- c(1,10,22)

r0 <- resid(t,conc,theta0)
sse <- sum(r0^2)
sse

a0 <- a.matrix(t,theta0)


theta1 <- theta0+solve(t(a0)%*%a0,t(a0)%*%r0)
theta1

r1 <- resid(t,conc,theta1)
sse <- sum(r1^2)
sse


theta1 <- theta0+0.5*solve(t(a0)%*%a0,t(a0)%*%r0)
theta1

r1 <- resid(t,conc,theta1)
sse <- sum(r1^2)
sse



theta1 <- theta0+0.2*solve(t(a0)%*%a0,t(a0)%*%r0)
theta1

r1 <- resid(t,conc,theta1)
sse <- sum(r1^2)
sse


a1 <- a.matrix(t,theta1)


theta2 <- theta1+0.2*solve(t(a1)%*%a1,t(a1)%*%r1)
theta2

r2 <- resid(t,conc,theta2)
sse <- sum(r2^2)
sse


a2 <- a.matrix(t,theta2)


theta3 <- theta2+0.2*solve(t(a2)%*%a2,t(a2)%*%r2)
theta3

r3 <- resid(t,conc,theta3)
sse <- sum(r3^2)
sse


a3 <- a.matrix(t,theta3)


theta4 <- theta3+solve(t(a3)%*%a3,t(a3)%*%r3)
theta4

r4 <- resid(t,conc,theta4)
sse <- sum(r4^2)
sse


a4 <- a.matrix(t,theta4)


theta5 <- theta4+solve(t(a4)%*%a4,t(a4)%*%r4)
theta5

r5 <- resid(t,conc,theta4)
sse <- sum(r5^2)
sse