library(R2jags)

file = '../DATA/theta.dat';
data = scan(file,skip=5,sep=' ')
y=na.omit(data)
N=length(y);

     
jags.data <- list("N","y")
Q=0.010000; R=0.040000; K=900.000000; r0=0.100000; theta=0.500000;

n.chains = 1
n.iter = 20000
n.burnin = floor(n.iter/2)
n.thin = max(1, floor(n.chains * (n.iter - n.burnin)/1000))

##################################
## Gamma priors on precision
##################################

jags.params=c("K","logr0","logtheta","iR","iQ","x");

jags.inits <- function(){
  list("K"=800,"logr0"=log(0.2),"logtheta"=log(1),"iQ"=1/0.05,"iR"=1/0.1,"x"=y,.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(12345)
time<-system.time(       
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,n.chains=n.chains, 
                  n.iter=n.iter, n.thin=n.thin,n.burnin=n.burnin,model.file="bugmodel-GammaPrior.txt")
)         
time <- unname(time["elapsed"]);

mcmc <- as.mcmc(jagsfit)

mcmcall <- mcmc[,-2]
mcmcall <- cbind(mcmcall[,1],mcmcall[,2],mcmcall[,3],mcmcall[,4],mcmcall[,5])

## get parameter estimates (and transform inverse variances back to log SDs)
coef <- apply(mcmcall,2,mean)
coef[2:3]=log(sqrt(1/coef[2:3]));
##coef[4:5]=exp(coef[2:3]);

coef.median<-apply(mcmcall,2,median)
coef.median[2:3]=log(sqrt(1/coef.median[2:3]));
##coef.median[4:5]=exp(coef.median[4:5]);

dens<-apply(mcmcall,2,density)
coef.mode<-sapply(dens,function(x) x$x[which.max(x$y)] )
coef.mode[2:3]=log(sqrt(1/coef.mode[2:3]));
##coef.mode[4:5]=exp(coef.mode[4:5]);

ci.mc <- t(apply(mcmcall,2,quantile, probs=c(0.025,0.975)))
## OBS remember that upper and lower are reversed!
ci.mc[2:3,2:1]=log(sqrt(1/ci.mc[2:3,]));
##ci.mc[4:5,]=exp(ci.mc[4:5,]);

##reorder the stuff (K,Q,R,r0,theta)=>(theta r0 K Q R)
myorder=c(5,4,1,2,3);
coef=coef[myorder];
coef.median=coef.median[myorder];
coef.mode=coef.mode[myorder];
ci.mc=ci.mc[myorder,];

results <- list(
                GP=list(
                obj=NULL,                     ## negative log L
                coef=coef,                   ## coefficients
                coef.median=coef.median,
                coef.mode=coef.mode,
                sd=NA,  
                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=list(effsize=effectiveSize(mcmc),
                  geweke=geweke.diag(mcmc))       ## MCMC convergence info
                ), UP=list() )

#####################################################################
## Repeat with uniform prior on log scale for variance parameters
#####################################################################

jags.params=c("K","logr0","logtheta","stdQ","stdR","x");

jags.inits <- function(){
  list("K"=800,"logr0"=log(0.2),"logtheta"=log(1),"stdQ"=sqrt(0.05),"stdR"=sqrt(0.1),"x"=y,.RNG.name="base::Wichmann-Hill", .RNG.seed=123)
}

set.seed(12345)
time<-system.time(       
  jagsfit <- jags(data=jags.data, inits=jags.inits, jags.params,n.chains=n.chains, 
                  n.iter=n.iter, n.thin=n.thin,n.burnin=n.burnin,model.file="bugmodel-UPrior.txt")
)         
time <- unname(time["elapsed"]);

mcmc <- as.mcmc(jagsfit)

mcmcall <- mcmc[,-2]
mcmcall <- cbind(mcmcall[,1],mcmcall[,2],mcmcall[,3],mcmcall[,4],mcmcall[,5])

## get parameter estimates (and transform back)
coef <- apply(mcmcall,2,mean)
##coef[2:3]=exp(coef[2:3]);
coef[4:5]=log(coef[4:5]);

coef.median<-apply(mcmcall,2,median)
##coef.median[2:3]=exp(coef.median[2:3]);
coef.median[4:5]=log(coef.median[4:5]);

dens<-apply(mcmcall,2,density)
coef.mode<-sapply(dens,function(x) x$x[which.max(x$y)] )
##coef.mode[2:3]=exp(coef.mode[2:3]);
coef.mode[4:5]=log(coef.mode[4:5]);

ci.mc <- t(apply(mcmcall,2,quantile, probs=c(0.025,0.975)))
##ci.mc[2:3,]=exp(ci.mc[2:3,]);
ci.mc[4:5,]=log(ci.mc[4:5,]);

##reorder the stuff (K,r0,theta,Q,R)=>(theta r0 K Q R)
myorder=c(3,2,1,4,5) 
coef=coef[myorder];
coef.median=coef.median[myorder];
coef.mode=coef.mode[myorder];
ci.mc=ci.mc[myorder,];

results[["UP"]]=list(
                obj=NULL,                     ## negative log L
                coef=coef,                   ## coefficients
                coef.median=coef.median,
                coef.mode=coef.mode,
                sd=NA,
                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=list(effsize=effectiveSize(mcmc),
                  geweke=geweke.diag(mcmc))
         );

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