Este template é baseado no livro do Kruschke Doing Bayesian data analysis:

É necessário a versão actual do OpenBUGS

library(BRugs)
## Welcome to BRugs connected to OpenBUGS version 3.2.3
#------------------------------------------------------------------------------
# THE MODEL.

# Specify the model in BUGS language, but save it as a string in R:
modelString = "
# BUGS model specification begins ...
model {
    # Likelihood:
    for ( i in 1:nFlips ) {
        y[i] ~ dbern( theta )
    }
    # Prior distribution:
    theta ~ dbeta( priorA , priorB )
    priorA <- 1
    priorB <- 1
}
# ... BUGS model specification ends.
" # close quote to end modelString

# Write the modelString to a file, using R commands:
writeLines(modelString,con="model.txt")
# Use BRugs to send the model.txt file to BUGS, which checks the model syntax:
modelCheck( "model.txt" )
## model is syntactically correct
#------------------------------------------------------------------------------
# THE DATA.

# Specify the data in R, using a list format compatible with BUGS:
dataList = list(
    nFlips = 14 ,
    y = c( 1,1,1,1,1,1,1,1,1,1,1,0,0,0 )
)

# Use BRugs commands to put the data into a file and ship the file to BUGS:
modelData( bugsData( dataList ) )
## data loaded
#------------------------------------------------------------------------------
# INTIALIZE THE CHAIN.

modelCompile()  # BRugs command tells BUGS to compile the model.
## model compiled
modelGenInits() # BRugs command tells BUGS to randomly initialize a chain.
## initial values generated, model initialized
#------------------------------------------------------------------------------
# RUN THE CHAINS.

# BRugs tells BUGS to keep a record of the sampled "theta" values:
samplesSet( "theta" )
## monitor set for variable 'theta'
# R command defines a new variable that specifies an arbitrary chain length:
chainLength = 10000
# BRugs tells BUGS to generate a MCMC chain:
modelUpdate( chainLength )
## 10000 updates took 0 s
#------------------------------------------------------------------------------
# EXAMINE THE RESULTS.

thetaSample = samplesSample( "theta" ) # BRugs asks BUGS for the sample values.
thetaSummary = samplesStats( "theta" ) # BRugs asks BUGS for summary statistics.
thetaSummary
##         mean     sd MC_error val2.5pc median val97.5pc start sample
## theta 0.7503 0.1043 0.001108   0.5226 0.7607     0.918     1  10000
hist(thetaSample)

samplesHistory("theta", mfrow = c(1, 1))

samplesDensity("theta", mfrow = c(1, 1))

samplesAutoC("theta", 1, mfrow = c(1, 1))

Exercício 1

modelString = "
model {
  for (i in 1:N) {
     r[i] ~ dbin(p[i], n[i])
     p[i] ~ dbeta(1,1)
  }
}
"
writeLines(modelString,con="modelo_ex1.txt")
modelCheck( "modelo_ex1.txt" )
## model is syntactically correct
########################################################
# Load Data
dataList = list(
  N=12, 
  n=c(47,148,119,810,211,196,148,215,207,97,256,360), 
  r=c(0,18,8,46,8,13,9,31,14,8,29,24)
)
modelData( bugsData( dataList ) )
## data loaded
modelCompile( numChains = 2 ) # we will use two chains
## model compiled
########################################################
# Load Initialization 
initData = list (      # in this case we have two chains
  c(0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5),
  c(0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1)
)

genInitFactory <- function()  {
  i <- 0
  function() {
    i <<- i + 1
    list(             # \
      p=initData[[i]] # | adapt to your specific exercise    
    )                 # /
  }
}

genInit <- genInitFactory()

for (chainIdx in 1:length(initData)) {
  modelInits( bugsInits( genInit ) )
}
## Initializing chain 1: 
## initial values loaded and chain initialized but another chain contain uninitialized variables
## Initializing chain 2: 
## model is initialized
########################################################
burn.in = 1000               # burn the first iterations
modelUpdate( burn.in )
## 1000 updates took 0 s
# Keep track of:
samplesSet( "p" )
## monitor set for variable 'p'
# Run the chain:
chainLength = 1000
modelUpdate( chainLength, thin = 10 )
## 1000 updates took 0 s
########################################################
# Check results:
N=12
pSample = NULL
pSummary = NULL
for (i in 1:N) {
  pSample  <- rbind( pSample,  samplesSample( paste0("p[",i,"]") ) )
  pSummary <- rbind( pSummary, samplesStats(  paste0("p[",i,"]") ) )
}
pSummary
##          mean       sd  MC_error  val2.5pc  median val97.5pc start sample
## p[1]  0.02049 0.019700 0.0004593 0.0005119 0.01488   0.07276  1001   2000
## p[2]  0.12590 0.027240 0.0005839 0.0777600 0.12480   0.18150  1001   2000
## p[3]  0.07466 0.023010 0.0005596 0.0357800 0.07264   0.12670  1001   2000
## p[4]  0.05773 0.007985 0.0001821 0.0432400 0.05767   0.07440  1001   2000
## p[5]  0.04233 0.013900 0.0003100 0.0192500 0.04099   0.07277  1001   2000
## p[6]  0.07077 0.018290 0.0003933 0.0389200 0.06959   0.10990  1001   2000
## p[7]  0.06707 0.020310 0.0003697 0.0326500 0.06538   0.11290  1001   2000
## p[8]  0.14750 0.024020 0.0005492 0.1014000 0.14630   0.19840  1001   2000
## p[9]  0.07204 0.017960 0.0003779 0.0399900 0.07086   0.11070  1001   2000
## p[10] 0.09047 0.028250 0.0006204 0.0423800 0.08834   0.15500  1001   2000
## p[11] 0.11620 0.019940 0.0004339 0.0788200 0.11620   0.15600  1001   2000
## p[12] 0.06935 0.013200 0.0002657 0.0459000 0.06868   0.09824  1001   2000
samplesHistory("*", mfrow = c(4, 3))

Segue um boxplot da simulaçao de cada um dos hospitais:

boxplot(pSample, use.cols=FALSE, xlab="Hospital", ylab="% morte")
abline(h=mean(pSummary$mean), col="red", lwd=1)

Exercício 2

modelString = "
# y[i] quantidade de oxigenio para o i-esimo paciente
# x[i] freq. cardiaca para o i-esimo paciente

model {
    # verosimilhança
    for (i in 1:N) {
        y[i] ~ dnorm(mu[i], tau)  
        mu[i] <- alpha + beta * (x[i] - xbar) # retirar a média ajuda nas auto-correlações
    }
    # distribuições à priori
    tau ~ dgamma(1.0E-6, 1.0E-6)
    sigma <- 1/tau
    alpha ~ dnorm(0.0, 1.0E-6)
    beta ~ dnorm(0.0, 1.0E-6)
    alpha0 <- alpha - xbar*beta   # repor o ajuste da média: esta é variável que temos de monitorar
}
"
writeLines(modelString,con="modelo_ex2.txt")
modelCheck( "modelo_ex2.txt" )
## model is syntactically correct
########################################################
# Load Data
dataList = list(
  x    = c(94,96,94,95,104,106,108,113,115,121,131),
  xbar = 107,   
  y    = c(0.47,0.75,0.83,0.98,1.18,1.29,1.4,1.6,1.75,1.9,2.23), 
  N    = 11
)
modelData( bugsData( dataList ) )
## data loaded
modelCompile( numChains = 2 ) # we will use two chains
## model compiled
########################################################
# Load Initialization 
alphas <- c(0.0, 10.0) # values for each chain
betas  <- c(0.0, 20.0) 
taus   <- c(1.0, 10.0)

genInitFactory <- function()  {
  i <- 0
  function() {
    i <<- i + 1
    list( 
      alpha = alphas[i],
      beta  = betas[i],
      tau   = taus[i]
    ) 
  }
}

genInit <- genInitFactory()

for (chainIdx in 1:length(initData)) {
  modelInits( bugsInits( genInit ) )
}
## Initializing chain 1: 
## initial values loaded and chain initialized but another chain contain uninitialized variables
## Initializing chain 2: 
## model is initialized
########################################################
burn.in = 1000               # burn the first iterations
modelUpdate( burn.in )
## 1000 updates took 0 s
# Keep track of:
samplesSet( c("alpha0", "beta", "sigma") )
## monitor set for variable 'alpha0'
## monitor set for variable 'beta'
## monitor set for variable 'sigma'
# Run the chain:
chainLength = 1000
modelUpdate( chainLength, thin = 10 )
## 1000 updates took 0 s
########################################################
# Check results:
sample.alpha <- samplesSample("alpha0")
stats.alpha  <- samplesStats("alpha0")
stats.alpha
##          mean     sd MC_error val2.5pc median val97.5pc start sample
## alpha0 -3.277 0.4088 0.007766   -4.091 -3.276    -2.471  1001   2000
sample.beta <- samplesSample("beta")
stats.beta  <- samplesStats("beta")
stats.beta
##         mean       sd MC_error val2.5pc median val97.5pc start sample
## beta 0.04283 0.003788 7.42e-05  0.03527 0.0428   0.05036  1001   2000
sample.sigma <- samplesSample("sigma")
stats.sigma  <- samplesStats("sigma")
stats.sigma
##         mean      sd  MC_error val2.5pc  median val97.5pc start sample
## sigma 0.0221 0.01442 0.0003127 0.008145 0.01844   0.05776  1001   2000
samplesHistory("*", mfrow = c(2, 2))

samplesAutoC("*", 1, mfrow = c(2, 2))
samplesAutoC("*", 2, mfrow = c(2, 2))