#----------------------------------------
# inferring p, predicting n1  -->  JAGS
# (with same output of inf_p_pred_gibbs.R etc)
# GdA, May 2020
#---------------------------------------


model = "tmp_model.bug"
write("
model{
    x0 ~ dbin(p, n0);
    x1 ~ dbin(p, n1);
     p ~ dbeta(1, 1);   
}
", model)

library(rjags)

data = list(n0=20, x0=10, n1=10)   

jm <- jags.model(model, data)
chain <- coda.samples(jm, c("p", "x1"), n.iter=10000)


p  <- as.vector(chain[[1]][,1])
x1 <- as.vector(chain[[1]][,2])

# analysis and plots of the chain
par(mfcol=c(2,2))
hist(p, nc=50, col='cyan', xlim=c(0,1))
barplot(table(x1)/N, col='cyan', xlab='x1')    #,xlim=c(0,n1))
plot(x1, p, col='blue', cex=1)
plot(NULL,  xaxt = 'n', yaxt = 'n', bty = 'n', pch = '',
     ylab = '', xlab = '', xlim=c(0,10), ylim=c(0,10))
text(4,9, sprintf("mean(p) = %.3f\n sigma(p) =  %.3f",
                  mean(p), sd(p)), cex=1.5,  col='blue')
text(4,5, sprintf("mean(x1) = %.3f\n sigma(x1) =  %.3f",
                  mean(x1), sd(x1)), cex=1.5, col='blue')
text(4,1, sprintf("rho(p,x1) = %.3f", cor(p,x1)), cex=1.5, col='blue')
par(mfcol=c(1,1))