library(Matrix)
library(numDeriv)

dat<-read.table('../DATA/min.dat', skip=3, header=FALSE)

nlogL<-function(theta){
  k<-exp(theta[1:3])
  sigma<-exp(theta[4])
  A<-rbind(
       c(-k[1], k[2]),
       c( k[1], -(k[2]+k[3]))
     )
  x0<-c(0,100)
  sol<-function(t)100-sum(expm(A*t)%*%x0)
  pred<-sapply(dat[,1],sol)
  -sum(dnorm(dat[,2],mean=pred,sd=sigma, log=TRUE))
}

time<-unname(system.time({
    fit<-nlminb(c(-2,-2,-2,-2),nlogL,hessian=TRUE)
    hes<-hessian(nlogL,fit$par)
  }
)["elapsed"])

coef=fit$par
sd=sqrt(diag(solve(hes)))
ci.q=coef+sd%o%c(-2,2)
obj=nlogL(coef)
maxgrad=max(abs(grad(nlogL,fit$par)))
ev <- eigen(hes)$value
eratio <- min(ev)/max(ev)

results <- list(obj=obj,                     ## negative log L
                coef=coef,                   ## coefficients
                sd=sd,                       ## standard deviations 
                confint.quad=ci.q,           ## Nx2 - matrix 95% ci 
                confint.profile=NULL,        
                confint.mcmc=NULL,           ## 0.025 0.975 quantiles
                time=c(fit=time,profile=NULL),  
                terminfo=c(maxgrad=maxgrad,eratio=eratio), ## MLE termination info
                convinfo=NULL       ## MCMC convergence info
                )

save("results",file="fit.RData")