library(R2jags)
load.module('msm') # the module containing the matrix exponential 

dat <- read.table('../DATA/min.dat', skip=3, header=FALSE)
time <- dat[,1]
M <- dat[,2]
noObs <- length(time)

set.seed(12345)
     
jags.data <- list("noObs","time","M")
jags.params <- c("logK1","logK2","logK3","sigma")
jags.inits <- function(){
  list("logK1"=-2+.2*rnorm(1),"logK2"=-2+.2*rnorm(1),"logK3"=-2+.2*rnorm(1), 
       "sigma"=exp(-2+.2*rnorm(1)), .RNG.name="base::Wichmann-Hill", .RNG.seed=round(runif(1)*1000000))
}

time<-unname(system.time(       
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params, 
                  n.iter=10000, n.thin=1,n.burnin=0,model.file="model.txt")
)["elapsed"])         

mcmc <- as.mcmc(jagsfit)
mcmcall <- do.call(rbind,mcmc[-c(1:1000),-1])
mcmcall <- cbind(mcmcall[,1:3],logSigma=log(mcmcall[,4]))

coef <- apply(mcmcall,2,mean)
sd <- apply(mcmcall,2,sd)
ci.q <- coef+sd%o%c(-2,2)
ci.mc <- t(apply(mcmcall,2,quantile, probs=c(0.025,0.975)))

mcmc<-as.mcmc(jagsfit)
mcmcall<-do.call(cbind,mcmc)

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

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