First Steps with OpenBUGS

Frequentist Approach

Let’s start with a simple example of a random experiment. Say we toss a coin 11 times and observe 3 heads and 8 tails. We then can estimate the probability of obtaining heads as follows.

The Binomial Model

Let \(y\) be the random vector representing the number of heads received when we toss a coin 11 times. Let \(\theta\) be the probability of obtaining heads. It follows that: \[y\mid \theta\;\sim\; \mathcal{B}(n,\theta)\] where \(\forall\,k=0,\ldots,n\) \[\mathbb{P}(y=k)={{n}\choose{k}} \theta^k (1-\theta)^{n-k}\]

The Maximum Likelihood Estimator

The Likelihood function of \(\theta\) is

\[l(y\mid \theta)={{n}\choose{y}} \theta^y (1-\theta)^{n-y}\] Taking the logarithm, we find: \[\log\left( l(y\mid \theta)\right)=y \log \theta +(n-y) \log (1-\theta)\]

Hence the Maximum Likelihood Estimator of \(\theta\) is defined by \[\widehat{\theta}=\mbox{argmax}_\theta \log\left( l(y\mid \theta)\right)\] Solving, we find: \[\widehat{\theta} =\displaystyle\frac{y}{n}\]

We can also compute the confidence interval of \(\theta\). It’s an asymptotic CI \[\left(\widehat{\theta}\pm z_{\alpha/2} \sqrt{\displaystyle\frac{1}{n}\widehat{\theta}(1-\widehat{\theta})}\right).\] where \(\alpha\) is the level of confidence.

The Bayesian approach

We consider a Beta distribution as a prior distribution:

\[\pi(\theta)\propto \theta^{a-1}(1-\theta)^{b-1}.\] Then the posteriori distribution on \(\theta\) is \[\begin{array}{rcl} \pi(\theta\mid y)&\propto &l(y\mid \theta) \pi(\theta) \\ & \propto & \theta^{y+a-1}(1-\theta)^{n-y+b-1}. e\end{array}\] Then \(y\mid \theta\sim \mbox{Beta}(y+a,n-y+b)\).

Let’s now show, using this example how can we use OpenBUGS to run the MCMC with probability distribution as the posterior distribution of \(\theta\).

Introducing OpenBUGS

OpenBUGS

  • This software is part of the BUGS project (Bayesian inference Using Gibbs Sampler) that aims to facilitate the use of MCMC methods for Bayesian Statistics. It’s developed by MRC Biostatistics of Cambridge University.

  • OpenBUGS is a free software.

  • OpenBUGS has a simple GUI with predefined models that can be constructed using DoodleBUGS.

  • We can also use OpenBUGS from R using R2OpenBUGS package.

  • OpenBUGS is based then on the Gibbs Sampler Algorithm.

1st code with OpenBUGS: Binomial Model

The model

model{
for(i in 1:n){
y[i] ~ dbern(p)
}
p~dbeta(a,b)
}

The data

list(n=11,y=c(0,0,0,0,0,0,1,0,0,1,0),a=2,b=3)

Initials

list(p=.2)
list(p=.9)
list(p=.4)

How to execute the code under OpenBUGS

  1. Check the model
  2. Load the data
  3. Specify the number of Markov Chains.
  4. Compile the model
  5. Initialize the model
  6. Generated the MC
  7. Specify which parameters to keep
  8. Check the convergence of the MC with coda

2nd Example

  • A fair coin is flipped 8 times. Can you estimate the probability of obtaining heads at most twice?

  • The model \[Y\sim\mbox{Binomial}(0.5,8)\] and we had to estimate \(\mathbb{P}(Y\leq 2)\)

  • The model under OpenBUGS

model{
  Y ~ dbin(0.5,8)  
  P <- step(2.5-Y)
}
  • step is a predefined function in BUGS that means P is equal to 1 if \(2.5-Y\geq 0\) and 0 otherwise.

OpenBUGS can be used to simulate one probability distribution

Example 1:

  • we would like a simulate a random sample with a Student-t probability distribution with mean 10, precision 2 and with degree of freedom 4:

  • The OpenBUGS code

model{
y ~ dt(10,2,4)
}

Example 2

  • Assume that we want to simulate a random sample from the random variable \(Y=(2Z+1)^3\) where \(Z\) follows a standard Normal distribution. We would like to estimate \(\mathbb{P}(Y>10)\)

  • The OpenBUGS code

model{
  z ~ dnorm(0,1)
  y <- pow(2*z+1,3)
  p <- step(y-10)
}

3rd Example: How many times should we repair?

  • Suppose that the cost of a reparation of a machine follows a Gamma distribution with mean $100 and standard deviation $50. How many times could we repair it with a budget of $1000? What’s the probability that this number exceeds 12?

  • We first easily check that a Gamma(4,0.04) probability distribution has a mean equal to $100 and standard deviation $50.

  • Let’s write first the mathematical model. Let \(Y_i, i=1,\ldots,I\) is the sequence of i.i.d random variables representing the cost of each reparation.

  • If we make \(M\) reparations the total cost of them can be represented by \[ \displaystyle\sum_{i=1}^M Y_i <1000\]

  • Since the average of each reparation is $100, assume that the maximum of the reparations cannot exceed 20. Thus, we can choose \(I=20.\)

  • Here’s the OpenBUGS code:

model{
  for(i in 1:20) {
    Y[i] ~ dgamma(4,0.04)
  }
  cum[1] <- Y[1]
  for(i in 2:20) {
    cum[i] <- cum[i-1]+Y[i]
  }
  for(i in 1:20) {
    cum.step[i] <- i*step(1000- cum[i])
  }
  number <- rank(cum.step[],20) ## le maximum dans cum.step
  P <- step(number-12)
}
  • Notice that we’ve utilized the function rank. Indeed, rank can be used as follows: rank(v,s), where v is a vector and s is a real number. rank will compute the number of components in v less than or equal to s.

Using now R2OpenBUGS, Binomial Model

  • We flip two coins 86 times. For the first coin, we obtain heads 83 times, and for the second, we obtain heads 72 times. Can we compare the probability of obtaining heads between these two coins?

  • Let’s first write the statistical model of this experience:

\[Y_1\sim\mbox{Binom}(86,p_1)\] and \[Y_2\sim\mbox{Binom}(86,p_2)\] Where \(Y_1\) and \(Y_2\) are the number of heads for each coin when it’s flipped 86 times. The data then is \((y_1,y_2)=(83,72)\)

  • Then unknown parameters in our model is the couple \((p_1,p_2)\in[0,1]^2\). We would like to use a Bayesian approach to estimate the posterior distribution of the difference: \[\delta=p_1-p_2\]

  • To do that, we will consider a uniform probability distribution on \([0,1]\) for both \(p_1\) and \(p_2\).

  • We are then able to write the BUGS code

model{
  for(i in 1:2){
    r[i] ~ dbin(p[i],n[i])
  }
  for (i in 1:2){ p[i] ~ dunif(0,1)}
  delta <- p[1] - p[2]
}

  • Let’s now show how to use the package R2OpenBUGS to perform a Bayesian analysis of this problem.

  • • First things to ensure are that:

    1. OpenBUGS is installed in your computer
    2. Install also from the CRAN the R2OpenBUGS

  • You can now use this code to write your BUGS code in a txt file:

sink('exBinom.txt')
cat("
model{
  for(i in 1:2){
    y[i] ~ dbin(p[i],n[i])
  }
  for (i in 1:2){ p[i] ~ dunif(0,1)}
  delta <- p[1] - p[2]
  }
",fill=T)
sink()
  • Write an R function that can be used to initialize the Markov Chain
inits<-function(){
  inits=list(p=runif(2,0,1))
}
  • We should specify which parameters have to be estimated:
params <- c('p','delta')
  • The data
n=c(86,86)
y=c(83,72)
data<-list("n","y")
  • Now, you can run the Markov Chain and get a simulation of the posterior distribution of the couple \((p_1,p_2)\)
library(R2OpenBUGS)
filename<-'exBinom.txt'
outBinom <-bugs(donnees,inits,params,filename,codaPkg=F,
                n.thin =1, n.iter=1000,debug=F,
                n.chains = 1,working.directory=getwd())
OpenBUGS.pgm="/Users/dhafermalouche/.wine/drive_c/Program Files (x86)/OpenBUGS/OpenBUGS323/OpenBUGS.exe"
outBinom <-bugs(data,inits,params,filename,codaPkg=F,
                n.thin =1, n.iter=1000,debug=F,
                n.chains = 1,working.directory=getwd(),
                OpenBUGS.pgm=OpenBUGS.pgm, useWINE=T)

where OpenBUGS.pgm is a variable in R that defines the path of the OpenBUGS software.

  • The result is then
> outBinom
## Inference for Bugs model at "exBinom.txt", 
## Current: 1 chains, each with 1000 iterations (first 500 discarded)
## Cumulative: n.sims = 500 iterations saved
##          mean  sd 2.5% 25% 50%  75% 97.5%
## p[1]      1.0 0.0  0.9 0.9 1.0  1.0   1.0
## p[2]      0.8 0.0  0.7 0.8 0.8  0.9   0.9
## delta     0.1 0.0  0.0 0.1 0.1  0.2   0.2
## deviance  9.3 2.1  7.3 7.9 8.7 10.0  14.9
## 
## DIC info (using the rule, pD = Dbar-Dhat)
## pD = 1.7 and DIC = 11.0
## DIC is an estimate of expected predictive error (lower deviance is better).