library(R2jags)
source("OrangeTree_BUGS_funs.R")

OrangeTreeData.y <- t(read.table("../DATA/OrangeTree.dat",skip=7))
OrangeTreeData.x<-scan(file='../DATA/OrangeTree.dat',what=numeric(), nlines=1, skip=5)
data<-list(n = 7, K = 5, x = OrangeTreeData.x,
Y = OrangeTreeData.y)

set.seed(1001) ## set RNG seed for R
inits<-list(phi1=200, phi2=800, phi3=400,tauC = 0.1,
            ## use different generator despite JN's nostalgia?
            .RNG.name="base::Wichmann-Hill",  ## set RNG seed/type for JAGS
            .RNG.seed=round(runif(1)*1000000))

st1 <- system.time(tfit_jags <- OrangeTree_BUGS_fit())["elapsed"]

tfit_jags_mcmc <- as.mcmc(tfit_jags)
coefs <- tfit_jags$BUGSoutput$mean
coefs <- unlist(coefs[names(coefs)!="deviance"])
coefs <- coefs[c(2:4,1)] ## !! rearrange order -- was alpha??
CI <- HPDinterval(tfit_jags_mcmc)
CI <- CI[rownames(CI)!="deviance",]
CI <- CI[c(2:4,1),]  ## ditto, rearrange
results <- list(coef=coefs,
                confint.mcmc=CI,
                time=c(fit=st1,mcmc=st1),
                convinfo=list(effsize=effectiveSize(tfit_jags_mcmc),
                  geweke=geweke.diag(tfit_jags_mcmc))
)

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