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

nr = 100000

#-------------  JAGS model --------------------------------------
library(rjags)
model = "tmp_model.bug"    # name of the model file ('temporary')
write("
model {
  nP.I   ~ dbin(pA, nP)          
  nV.A  ~ dbin(pA, nV)         # quellli che si sarebbero beccati
                                # il Covid e il vaccino non avesse funzionato
  pA ~ dbeta(1,1)
  nV.I ~ dbin(ffe, nV.A)      # ffe = 1 - eff 
  ffe ~ dbeta(1,1)
  eff <- 1 - ffe
}
", model)

# Moderna
nP   = 15000     # placebo     NON CRUCIALI [-> 150, 1500, 15000]
nV   = 15000     # vaccino       "     "
nP.I = 90        # infetti placebo 
nV.I = 5         # infetti vaccino
#------------  JAGS call via rjags ------------------------------
data <- list(nP=nP, nV=nV, nP.I=nP.I, nV.I=nV.I)
jm <- jags.model(model, data)
update(jm, 100)
to.monitor <-  c('pA', 'eff', 'nV.A')
chain <- coda.samples(jm, to.monitor, n.iter=nr)
chain.df.Moderna <- as.data.frame( as.mcmc(chain) )

print(summary(chain))
plot(chain)

pause()

hist(chain.df.Moderna$eff, nc=100, col='cyan', freq=FALSE,
     xlab='efficacy', main='Moderna')
abline(v=0.945, col='blue')

pause()


# Pfizer
nP   = 20000     # placebo   NON CRUCIALI
nV   = 20000     # vaccino    "    "
nP.I = 162       # infetti placebo 
nV.I = 8         # infetti vaccino
#------------  JAGS call via rjags ------------------------------
data <- list(nP=nP, nV=nV, nP.I=nP.I, nV.I=nV.I)
jm <- jags.model(model, data)
update(jm, 100)
to.monitor <-  c('pA', 'eff', 'nV.A')
chain <- coda.samples(jm, to.monitor, n.iter=nr)
chain.df.Pfizer <- as.data.frame( as.mcmc(chain) )

print(summary(chain))
plot(chain)

pause()

hist(chain.df.Pfizer$eff, nc=100, col='cyan', freq=FALSE,
     xlab='efficacy', main='Pfizer')
abline(v=0.95, col='blue')