The predictive distribution for a Binomial sampling model with Beta prior
Let \(R \sim \mbox{Binomial}(n,\theta)\), where \(n\) denotes the number of trials. Suppose that the prior of \(\theta\) is \(\mbox{Beta}(a,b)\). Then, the posterior predictive distribution of \(r^\star\) successes out of \(n^\star\) trials (given the data) is:
\[\begin{align*} \pi(r^\star|r) &= \int_0^1 \pi_{R^\star}(r^\star|\theta) \pi_{\Theta|r}(\theta|r)\, d\theta \\ &= \int_0^1 \begin{pmatrix}n^\star\\r^\star\end{pmatrix} \theta^{r^\star} (1-\theta)^{(n^\star-r^\star)} \mbox{Beta}(a+r,b+n-r) \, d\theta \\ & = \begin{pmatrix}{n^\star}\\r^\star\end{pmatrix} \frac{B(r^\star+a+r,b+n-r+n^\star-r^\star)}{B(a+r,b+n-r)}, \end{align*}\]where \(B(\cdot,\cdot)\) denotes the Beta function. This distribution is known as the Beta-Binomial distribution. The following R code implements this distribution and illustrates the shape of it for \(n=10,1000\) and \(r=3,300\), and several values of the hyperparameters \((a,b)\).
# Beta Binomial Predictive distribution function
BetaBinom <- Vectorize(function(rp){
log.val <- lchoose(np, rp) + lbeta(rp+a+r,b+n-r+np-rp) - lbeta(a+r,b+n-r)
return(exp(log.val))
})
#####################
# Small sample
#####################
# Example 1
n <- 10; r <- 3; a <- 1; b <- 1; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=1, b=1",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 2
n <- 10; r <- 3; a <- 1; b <- 0.001; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=1, b=0.001",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 3
n <- 10; r <- 3; a <- 4; b <- 4; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=4, b=4",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 4
n <- 10; r <- 3; a <- 7; b <- 2; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=7, b=2",cex.axis= 1.5,cex.lab=1.5,lwd=4)
#####################
# Large sample
#####################
# Example 5
n <- 1000; r <- 300; a <- 1; b <- 1; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=1, b=1",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 6
n <- 1000; r <- 300; a <- 1; b <- 0.001; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=1, b=0.001",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 7
n <- 1000; r <- 300; a <- 4; b <- 4; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=4, b=4",cex.axis= 1.5,cex.lab=1.5,lwd=4)
# Example 8
n <- 1000; r <- 300; a <- 7; b <- 2; np <- 10
plot(1:10,BetaBinom(1:10),type="b",xlab="r*",ylab="P(R=r*|Data)", main = "Posterior predictive: a=7, b=2",cex.axis= 1.5,cex.lab=1.5,lwd=4)