pause <- function() { cat ("\n >> press Enter to continue\n"); scan() }

call.jags <- function(model, data, nr) {
  jm <- jags.model(model, data)
  update(jm, 100)
  chain <- coda.samples(jm, 'p', n.iter=nr)
  print(summary(chain))
  chain.df <- as.data.frame( as.mcmc(chain) )
  return(chain.df$p)
}

library(rjags)
model = "tmp_model.bug"    # name of the model file ('temporary')
write("
model {
  n ~ dbin(p, N)
  p ~ dbeta(r0,s0)   
}
", model)

nr = 100000
N = 10; n = 3 
r0 = s0 = 1    # flat prior

# First call to JAGS
data <- list(N=N, n=n, r0=r0, s0=s0)
p <- call.jags(model, data, nr)
pause()

h <- hist(p, freq=FALSE, nc=100, col='cyan', xlim=c(0,1), ylim=c(0,4))
p.m <- sum(h$mids*h$counts)/sum(h$counts)
p.s <- sqrt(sum(h$mids^2*h$counts)/sum(h$counts) - p.m^2)
cat(sprintf(">>> JAGS:  %.4f +- %.4f\n", p.m, p.s))
pause()

# overimpose the posterion got from a Beta conjugate
curve(dbeta(x,r0+n,s0+(N-n)), col='blue', add=TRUE)
pause()

# Not-flat prior used to reshape the JAGS result
curve(dbeta(x,r,s), col='magenta', add=TRUE)
pause()

# reweighing
w <- dbeta(h$mids, r, s)
f.p.w <- h$density*w
f.p.w <- f.p.w/sum(f.p.w)/(h$mids[2]-h$mids[1])
points(h$mids, f.p.w, col='blue', ty='l', lwd=2)     
p.w.m <- sum(h$mids*f.p.w)/sum(f.p.w)
p.w.s <- sqrt(sum(h$mids^2*f.p.w)/sum(f.p.w) - p.w.m^2)
cat(sprintf(">>> Reweighed:  %.4f +- %.4f\n\n", p.w.m, p.w.s))
pause()

# second call to JAGS, with the new prior
data <- list(N=N, n=n, r0=r, s0=s)
p1 <- call.jags(model, data, nr)
h1 <-  hist(p1, nc=100, plot=FALSE)
points(h1$mids, h1$density, col='red', ty='l', lwd=1.5)  
p1.m <- sum(h1$mids*h1$counts)/sum(h1$counts)
p1.s <- sqrt(sum(h1$mids^2*h1$counts)/sum(h1$counts) - p1.m^2)
cat(sprintf(">>> JAGS(%d,%d):  %.4f +- %.4f\n", r,s, p1.m, p1.s))
pause()

# New Beta obtained by the well known updating rule
curve(dbeta(x,r+n,s+(N-n)), col='green', add=TRUE)