Code for Lecture 3

This is a very simple example of a R Markdown file, where you will find all the code used during the lecture today. Just run this into R using the packages rmarkdown and knitr.

## At some point we will need these libraries
## 
library(reshape2)
library(plyr)
library(xtable)

library(ggplot2)
library(GGally)
library(RColorBrewer)

The next chunk of code shows how to plot the binomial distributions for different values of the parameter \(\theta\).

set.seed(1234) ## to ensure reproducibility of RNGs

n=30 ## number of trials
p=c(0.3,0.5,0.7,0.9) ## values of theta

## create a function to plot the binomial distributions

graph <- function(n,p){
  x <- dbinom(0:n,size=n,prob=p)
  barplot(x,names.arg=0:n,
      main=substitute(paste(theta, "=", pp),list(pp = p)))
}

par(mfrow=c(2,2))
mapply(graph,30,p)

##       [,1] [,2] [,3] [,4]
##  [1,]  0.7  0.7  0.7  0.7
##  [2,]  1.9  1.9  1.9  1.9
##  [3,]  3.1  3.1  3.1  3.1
##  [4,]  4.3  4.3  4.3  4.3
##  [5,]  5.5  5.5  5.5  5.5
##  [6,]  6.7  6.7  6.7  6.7
##  [7,]  7.9  7.9  7.9  7.9
##  [8,]  9.1  9.1  9.1  9.1
##  [9,] 10.3 10.3 10.3 10.3
## [10,] 11.5 11.5 11.5 11.5
## [11,] 12.7 12.7 12.7 12.7
## [12,] 13.9 13.9 13.9 13.9
## [13,] 15.1 15.1 15.1 15.1
## [14,] 16.3 16.3 16.3 16.3
## [15,] 17.5 17.5 17.5 17.5
## [16,] 18.7 18.7 18.7 18.7
## [17,] 19.9 19.9 19.9 19.9
## [18,] 21.1 21.1 21.1 21.1
## [19,] 22.3 22.3 22.3 22.3
## [20,] 23.5 23.5 23.5 23.5
## [21,] 24.7 24.7 24.7 24.7
## [22,] 25.9 25.9 25.9 25.9
## [23,] 27.1 27.1 27.1 27.1
## [24,] 28.3 28.3 28.3 28.3
## [25,] 29.5 29.5 29.5 29.5
## [26,] 30.7 30.7 30.7 30.7
## [27,] 31.9 31.9 31.9 31.9
## [28,] 33.1 33.1 33.1 33.1
## [29,] 34.3 34.3 34.3 34.3
## [30,] 35.5 35.5 35.5 35.5
## [31,] 36.7 36.7 36.7 36.7

The following chunk of code shows how to obtain plots of the prior, likelihood and posterior distributions in the Beta-Binomial model, using a grid of values for \(\theta \in (0,1)\).

n      <- 30 ## number of trials
Y      <- 20 ## number of successes
a      <- 2 ## prior hyperp.
b      <- 5

grid   <- seq(0,1,.01)

like   <- dbinom(Y,n,grid)
like   <- like/sum(like) #standardized to scale the figure

prior  <- dbeta(grid,a,b)
prior  <- prior/sum(prior) #standardized to scale the figure

post   <- like*prior
post   <- post/sum(post) 

  plot(grid,like,type="l",lty=2, 
     xlab="theta",ylab="Density") 
  lines(grid,prior)
  lines(grid,post,lwd=2, col="blue")
  legend("topleft",c("Likelihood","Prior","Posterior"),
                  lwd=c(1,1,2),lty=c(2,1,1), col=c("black", "black", "blue"))

The following code shows how to plot the posterior distribution of the \(\theta\) parameter in the Beta-Binomial model using the exact posterior distribution (Beta with updated parameters)

post <- dbeta(grid,Y+a,n-Y+b)
plot(grid,post,type="l",
     xlab="theta",
     ylab="Posterior density")

The next code shows how to obtain the posterior predictive distribution (i.e. the predicted distribution of the successes in the next \(m\) trials in the Beta-Binomial model.

m=10 ## let's see what may happen in the future 10 matches
post_samples = rbeta(10000,Y+a,n-Y+b)
## Generate one predictive value from the Binomial from each theta from the posterior
pred_samples=sapply(post_samples, rbinom, size=10, n=1) ## size=n. trials
hist(pred_samples, col="blue", probability=TRUE, main="predictive distribution", 
     xlab="predictive samples")