madisonModel=function() { 
  for (c in 1:TotalCycles) {  
    y[c] ~ dbern(prob[id[c],tx[c]]) 
  } 

  for (w in 1:Women) { 
       prob[w,1] ~  dbeta(a,b) 
       pr[w,2]   <- rr*prob[w,1] 
       prob[w,2] <- min(pr[w,2],0.99) 
  } 

 a  ~ dunif(0,100) 
 b  ~ dunif(0,100) 
 rr ~ dunif(0,10) 
}

library(BRugs)

setwd("c:/tmp")

writeModel(madisonModel,"madisonModel.txt")

modelCheck("madisonModel.txt")        # check model file

modelData("madisonData.txt")          # read data file
modelData("madisonSampleSize.txt")    # read data file

modelCompile(numChains=1)             # compile model with 1 chain

modelInits("madisonInitialValues.txt",1)  # read init data file
modelGenInits()  # read init data file

modelUpdate(50000)                    # burn in

samplesSet(c("a", "b","rr"))         # a b and rr should be monitored

modelUpdate(25000)                    # 1000 more iterations ....

samplesStats("*")                    # the summarized results

res<-samplesStats("*")                    # the summarized results


a.hat=res[1,1] ; b.hat=res[2,1];

p.0=seq(0.01,0.99,0.01)
plot(p.0,dbeta(p.0,a.hat,b.hat),,type="l" )

## some plots

samplesHistory("a") # plot the chain,
samplesHistory("b") # plot the chain,
samplesHistory("rr") # plot the chain,
samplesDensity("rr")              

# plot the densities,

par(mfrow=c(2,1))