library(rbugs,lib.loc="~/myRpackages")
library(coda,lib.loc="~/myRpackages")

options("digits"=22)

################
## Theta = 0.5
################

file = '../data/parest_thet05.txt';
Q=0.010000; R=0.040000; K=900.000000; r0=0.100000; theta=0.500000;
params=c(r0,K,sqrt(Q),sqrt(R),theta);

data = read.table(file,skip=2,sep=',',header=FALSE,colClasses=c("character","character"));

N = nrow(data);
y = as.numeric(data[,2]);
x = as.numeric(data[,1]);

wang.data=list("N","y");
parameters=c("K","logr0","logtheta","logQ","logR","x");

wang.bug=file.path("~/wang/BUGS/final.bug");

inits = list(list(K=800,logr0=log(0.2),logtheta=log(1.0),logQ=log(0.05),logR=log(0.1),x=y)) 

niter=100000;

system.time(wang.run<-rbugs(data=wang.data,inits,parameters,wang.bug,n.chains=1, n.iter=niter,bugs="~/OpenBUGS310/bin/OpenBUGS", #workin
                            ,bugsWorkingDir=".",OpenBugs=TRUE,verbose=TRUE,seed=66)) 


ll=list();
for(i in 1:length(wang.run))
    {
      ll[[i]]=as.mcmc(wang.run[[i]][,-c(1,3)]);
    }
tt=as.mcmc.list(ll)

## plot chains and posterior distr.
plot(tt[[1]])

## Geweke test for convergence
Z=geweke.diag(tt[[1]])
qqnorm(Z$z)
KS=ks.test(Z$z,"pnorm"); ## test null hypothesis that Z is normal distributed.
KS$p.value


# MAP estimat for K=900
dens=density(tt[[1]][,1])
dens$x[which.max(dens$y)]
quantile(tt[[1]][,1],c(0.025,0.975)) #credible interval
# MAP estimat for Q=0.01  
dens=density(tt[[1]][,2])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,2],c(0.025,0.975))) #credible interval
# MAP estimat for R=0.04
dens=density(tt[[1]][,3])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,3],c(0.025,0.975))) #credible interval
# MAP estimat for r0=0.1
dens=density(tt[[1]][,4])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,4],c(0.025,0.975))) #credible interval

## MAP estimat for theta=0.5
dens=density(tt[[1]][,5])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,5],c(0.025,0.975))) #credible interval


## Write results to file
out=summary(tt);
write.table(out$quantiles,"./bugs_results_thet05.csv");


################
## Theta = 1.5
################

file = '../data/parest_thet15.txt';
Q=0.010000; R=0.040000; K=900.000000; r0=0.100000; theta=1.500000;
params=c(r0,K,sqrt(Q),sqrt(R),theta);

data = read.table(file,skip=2,sep=',',header=FALSE,colClasses=c("character","character"));

N = nrow(data);
y = as.numeric(data[,2]);
x = as.numeric(data[,1]);

wang.data=list("N","y");
parameters=c("K","logr0","logtheta","logQ","logR","x");

wang.bug=file.path("~/wang/BUGS/final.bug");
wang.bug=file.path("final.bug");

inits = list(list(K=800,logr0=log(0.2),logtheta=log(1.0),logQ=log(0.05),logR=log(0.1),x=y)) 

niter=100000;

system.time(wang.run<-rbugs(data=wang.data,inits,parameters,wang.bug,n.chains=1, n.iter=niter,bugs="~/OpenBUGS310/bin/OpenBUGS", #workin
                            ,bugsWorkingDir=".",OpenBugs=TRUE,verbose=TRUE,seed=66)) 


ll=list();
for(i in 1:length(wang.run))
    {
      ll[[i]]=as.mcmc(wang.run[[i]][,-c(1,3)]);
    }
tt=as.mcmc.list(ll)

## plot chains and posterior distr.
## plot(tt[[1]])

## Geweke test for convergence
Z=geweke.diag(tt[[1]])
qqnorm(Z$z)
KS=ks.test(Z$z,"pnorm"); ## test null hypothesis that Z is normal distributed.
KS$p.value


# MAP estimat for K=900
dens=density(tt[[1]][,1])
dens$x[which.max(dens$y)]
quantile(tt[[1]][,1],c(0.025,0.975)) #credible interval
# MAP estimat for Q=0.01  
dens=density(tt[[1]][,2])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,2],c(0.025,0.975))) #credible interval
# MAP estimat for R=0.04
dens=density(tt[[1]][,3])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,3],c(0.025,0.975))) #credible interval
# MAP estimat for r0=0.1
dens=density(tt[[1]][,4])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,4],c(0.025,0.975))) #credible interval

## MAP estimat for theta=1.5
dens=density(tt[[1]][,5])
exp(dens$x[which.max(dens$y)])
exp(quantile(tt[[1]][,5],c(0.025,0.975))) #credible interval


## Write results to file
out=summary(tt);
write.table(out$quantiles,"./bugs_results_thet15.csv");