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.

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 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.

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\).

`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.

`OpenBUGS`

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

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

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

`OpenBUGS`

- Check the model
- Load the data
- Specify the number of Markov Chains.
- Compile the model
- Initialize the model
- Generated the MC
- Specify which parameters to keep
- Check the convergence of the MC with
`coda`

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 distributionwe 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)
}
```

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)
}
```

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`

.

`R2OpenBUGS`

, Binomial ModelWe 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:

`OpenBUGS`

is installed in your computer- 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())
```

- If you’re using Mac OS, it’s now possible to use OpenBUGS. I advise to follow the directions in the following tutorial: https://oliviergimenez.github.io/post/run_openbugs_on_mac/. I’m using a Mac and It worked for me perfectly. You have then to use this code:

```
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).
```