#----------------------------------------
# inferring p, predicting n1  -->  Gibbs
#
# GdA, May 2020
#---------------------------------------

# initialize the history vectors
N = 10000; p = x1 = rep(0, N) 

# observed nodes
n0 = 20; x0 = 10; n1 = 10

# initial p and n1 
p[1] = rbeta(1,1,1)           # uniform
x1[1] = rbinom(1, n1, p[1])

for (i in 2:N) {
   p[i] = rbeta(1, x0+x1[i-1]+1, n0+n1-x0-x1[i-1]+1)
   x1[i] = rbinom(1, n1, p[i])
}

# 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))