#---- data and parameters nr = 1000000 ns = 10000 nP = 2010 r0 = s0 = 1 r1 = 409.1; s1 = 9.1 r2 = 25.2; s2 = 193.1 #---- JAGS model ------------------------------ model = "tmp_model.bug" # name of the model file ('temporary') write(" model { nP ~ sum(nP.I, nP.NI) nP.I ~ dbin(pi1, n.I) nP.NI ~ dbin(pi2, n.NI) pi1 ~ dbeta(r1, s1) pi2 ~ dbeta(r2, s2) n.I ~ dbin(p, ns) n.NI <- ns - n.I p ~ dbeta(r0,s0) } ", model) #---- call JAGS --------------------------------------------- library(rjags) data <- list(ns=ns, nP=nP, r0=r0, s0=s0, r1=r1, s1=s1, r2=r2, s2=s2) jm <- jags.model(model, data) update(jm, 10000) to.monitor <- c('p', 'n.I') chain <- coda.samples(jm, to.monitor, n.iter=nr) #---- show results print(summary(chain)) plot(chain, col='blue')