## Warning: package 'knitr' was built under R version 3.5.3
\[ \begin{array}{l} f(y|\theta) = \binom{N}{y} \theta^y (1-\theta)^{N-y} \\ \pi(\theta) = \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1} \end{array} \]
\(\Longrightarrow\)
\[ \begin{array}{lcl} \pi(\theta|y) &=& {\displaystyle}\frac{{\displaystyle}f(y|\theta) \pi(\theta)}{{\displaystyle}f(y)} = {\displaystyle}\frac{f(y|\theta) \pi(\theta)}{{\displaystyle}\int_0^1 f(y|\theta) \pi(\theta) d\theta} \\ &=& {\displaystyle}\frac{{\displaystyle}\binom{N}{y} \theta^y (1-\theta)^{N-y} \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1}} {{\displaystyle}\int_0^1 \binom{N}{y} \theta^y (1-\theta)^{N-y} \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1} d\theta} \\ &=& {\displaystyle}\frac{1}{{\displaystyle}B(y+a,N+y-b)} \ \theta^{y+a-1} (1-\theta)^{N-y+b-1} \end{array} \]
par(mfrow=c(2,2))
N1<-10; y1<-2
N2<-50; y2<-10
theta <- seq(0.001,0.999,0.001)
# non-informative prior Beta(1,1)
a <- 1; b<-1
plot(theta, dbeta(theta, a, b), type="l", ylab="pdf", ylim=c(0,8), main="(A) Prior Beta(1,1)")
points(theta, dbeta(theta, y1+a, N1-y1+b), type="l", lty=2, col="blue")
points(theta, dbeta(theta, y2+a, N2-y2+b), type="l", lty=3, col="red")
# Informative prior Beta(2,1)
a <- 2; b<-1
plot(theta, dbeta(theta, a, b), type="l", ylab="pdf",ylim=c(0,8), main="(B) Prior Beta(2,1)")
points(theta, dbeta(theta, y1+a, N1-y1+b), type="l", lty=2, col="blue")
points(theta, dbeta(theta, y2+a, N2-y2+b), type="l", lty=3, col="red")
# Informative prior Beta(2,3)
a <- 2; b<-3
plot(theta, dbeta(theta, a, b), type="l", ylab="pdf",ylim=c(0,8), main="(C) Prior Beta(2,3)")
points(theta, dbeta(theta, y1+a, N1-y1+b), type="l", lty=2, col="blue")
points(theta, dbeta(theta, y2+a, N2-y2+b), type="l", lty=3, col="red")
# Non-informative prior Beta(1/2,1/2)
a <- b <- 1/2
plot(theta, dbeta(theta, a, b), type="l", ylab="pdf",ylim=c(0,8), main="(D) Prior Beta(1/2,1/2)")
points(theta, dbeta(theta, y1+a, N1-y1+b), type="l", lty=2, col="blue")
points(theta, dbeta(theta, y2+a, N2-y2+b), type="l", lty=3, col="red")
CI_1 <- function(a, b, N, y, level){
percentile <- c((1-level)/2, 1-(1-level)/2)
CI <- qbeta(percentile, y+a, N-y+b)
CI
}
CI_2 <- function(a, b, N, y, level, split=1000){
theta <- seq(0,1,length.out=split+1)[-c(1,(split+1))]
PDF <- dbeta(theta, y+a, N-y+b)/split
clist <- seq(0.00001,max(PDF),length.out=100)
SS <- numeric(100)
for(i in 1:100){ SS[i] <- sum(PDF[PDF > clist[i]]) }
Num <- which.min(abs(SS-level))
INT <- theta[PDF > clist[Num]]
CI <- c(min(INT), max(INT))
CI
}
#Example data
N1 <- 20
y1 <- 4
a <- b <- 1
#Posterior
par(mfrow=c(1,1))
plot(theta, dbeta(theta, y1+a, N1-y1+b), type="l", ylab="pdf")
theta.mean<-(a+y1)/(N1+a+b)
points(theta.mean, 0, pch=4, col="red")
theta.mode<-(a+y1-1)/(N1+a+b-2)
points(theta.mode, 0, pch=4, col="blue")
lines(c(theta.mode, theta.mode), c(0,dbeta(theta.mode, y1+a, N1-y1+b)), lty=2, col="blue")
CI_1.res <- CI_1(a=a, b=b, N=N1, y=y1, level=0.95)
CI_1.res
[1] 0.08217588 0.41906604
lines(CI_1.res, c(0,0), lwd=2, col="green")
CI_2.res <- CI_2(a=a, b=b, N=N1, y=y1, level=0.95)
CI_2.res
[1] 0.070 0.398
lines(CI_2.res, c(0.1,0.1), lwd=2, col="blue")
Simulation-based posterior calculaiton
- Grid method - Sampling Importance-resampling (SIR) method
- Markov Chain Monte carlo (MCMC) methods
etc.
For \(\theta_1 < \theta_2 < \cdots < \theta_k\), evaluate the likelihood and prior, and then compute the following values for obataining the posterior numerically. \[ \pi(\theta_i | y) = {\displaystyle}\frac{f(y|\theta_i) \pi(\theta_i)}{{\displaystyle}\sum_{j=1}^k f(y|\theta_j) \pi(\theta_j)} \]
N1 <- 10
y1 <- 2
a <- b <- 1
Increment <- 0.001
theta <- seq(0.001,0.999,Increment)
like <- dbinom(y1, N1, theta)
prior <- dbeta(theta, a, b)
posterior1 <- like*prior/sum(Increment*like*prior)
Note: If \(g(\theta)=\pi(\theta)\), then the weight is just proportional to the likelihood as follows: \[ {\displaystyle}\omega_i = \frac{{\displaystyle}f(y|\theta_i) }{{\displaystyle}\sum_{j=1}^n f(y|\theta_j) }. \]
#proposal = prior Beta(1,1)
MM <- 10^5
NN <- 10^3
theta.sample <- rbeta(MM, a, b)
weight <- dbinom(y1, N1, theta.sample)
theta.resample <- sample(theta.sample, NN, prob=weight, replace=T)
plot(theta, dbeta(theta, a+y1, b+N1-y1), type="l", ylab="Denstiy")
points(theta, posterior1, type="l", lty=2, col="blue") #overlaid
hist(theta.resample, add=T, freq=F) #perfect!
Perhaps we will use WinBUGS, Jags, Stan etc…