AST 501 Bayesian Statistics

Chapter 01

Fig 1.1 (left)

n<-20 
x<-0:n 
del<-0.25 

plot(range(x-del), c(0,0.4), xlab = "number infected in the sample", ylab = "probability", type = "n")

points(x-del, dbinom(x, n, 0.05), type = "h", col = gray(0.75), lwd=2) 
points(x, dbinom(x, n, 0.1), type = "h", col = gray(0.5), lwd=2)
points(x+del, dbinom(x, n, 0.05), type = "h", col = gray(0), lwd=2)

legend(10, 0.4, legend = c(
  expression(paste(theta," = 0.05", sep = "")),
  expression(paste(theta," = 0.1", sep = "")),
  expression(paste(theta," = 0.2", sep = ""))),
  lwd = c(2,2,2),
  col = gray(c(0.75,0.5,0)),
  bty = "n"
  )

library(ggplot2)
n<-20
x<- 0:20
del<-0.2
y1<-dbinom(x = x,size = 20,prob = .05)
y2<-dbinom(x = x,size = 20,prob = .1)
y3<-dbinom(x = x,size = 20,prob = .2)


ggplot()+
  geom_linerange(mapping = aes(x = x-del,y = y1,ymin = 0,ymax = y1,colour = "Y1"))+
  geom_linerange(mapping = aes(x = x,y = y2,ymin = 0,ymax = y2,colour = "Y2"))+
  geom_linerange(mapping = aes(x = x+del,y = y3,ymin = 0,ymax = y3,colour = "Y3"))

plot(x-del, y1,type = "h" ,col="red1")
points(x,y2,type= "h",col="yellow")
points(x,y3,type= "h",col="blue")

Figure 1.1 (Right)

We plot the prior and posterior on the same axis where p(θ) ~ beta(2,20) and p(θ|y) ~ beta(2,40). If the prior p(θ) ~ beta(a,b) and the sampling model p(y|θ) ~binomial(n,θ) then the posterior is p(θ|y)~beta(a+y,b+n−y)

We set the prior parameters, observed data, and generate theta for the x axis for plotting the prior.

a<-2
b<-20
y<-0
theta<-seq(0, 1, length.out=500)

plot(theta, dbeta(theta, a, b), 
     type="l",
     xlab=expression(paste("percentage infected in the population", theta,
                           sep = " ")),
     ylab = "probability",
     lwd = 2,
     col = "gray",
     ylim = c(0,16))

lines(theta, dbeta(theta, a+y, b+n-y), lwd = 2)
legend(0.6, 14, legend = c(
  expression(paste(italic("p"),"(",theta,")", sep = "")),
  expression(paste(italic("p"),"(",theta,"|",italic("y"),")", sep = ""))),
  lwd = c(2,2),
  col = c("grey", "black"),
  bty="n")

## Lecture 2

pbeta(0.2,2,20)-pbeta(0.05,2,20)

## [1] 0.6593258

pbeta(0.1,2,40)

## [1] 0.9260956

mean<- 2/(2+40)
mean

## [1] 0.04761905

x<-seq(0,1,length=10000)
y<-dbeta(x,2,20) # prior
plot(x,y)

z<-dbeta(x,2,40) # posterior
plot(x,z)

ggplot()+
  geom_line(aes(x = x,y = y,colour="prior"))+
  geom_line(aes(x = x,y = z,colour="posterior"))

Figure 1.2 (Left)

library("graphics")

 w<- seq(0, 25, length.out = 20)
theta_0<- seq(0, 0.5, length.out = 20)

y<-0
a<-2
b<-20
z <- matrix(0, nrow = length(w), ncol = length(theta_0))
for (i in 1: length(w)){
  for( j in 1: length(theta_0)){
    z[i,j]<- (w[i]*theta_0[j])/(w[i]+20)
  }
}



contour(x = w,
        y = theta_0,
        z)

Figure 1.2 (Right)

library("graphics")

 w<- seq(0, 25, length.out = 20)
theta_0<- seq(0, 0.5, length.out = 20)

y<-0
a<-2
b<-20
z <- matrix(0, nrow = length(w), ncol = length(theta_0))

for (i in 1: length(w)){
  for( j in 1: length(theta_0)){
    z[i,j]<- pbeta(q = 0.1,shape1 = (w[i]*theta_0[j]),shape2 = (w[i]*(1-theta_0[j])+20))
  }
}



contour(x = w,
        y = theta_0,
        z)

Chapter 02

Figure 2.1 (Left)

x<-0:10
plot(x,dpois(x,2.1), type="h",lwd=1, xlab=expression(italic(y)), ylab= expression(paste(italic("p"),"(",italic("y"),"|",theta==2.1,")"))) 

Figure 2.1 (Right)

x<-0:100
plot(x,dpois(x,21), type="h",lwd=1, xlab=expression(italic(y)),ylab=
  expression(paste(italic("p"),"(",italic("y"),"|",theta==21,")")) )

Figure 2.2 (Left)

mu<-10.75
sig<- .8
x<- seq(7.9,13.9,length=500)

plot(x,pnorm(x,mu,sig),type="l",ylab=expression(paste(italic("F"),"(",italic("y"),")")),xlab=
    expression(italic(y)),lwd=1)
abline(h=c(0,.5,1),col="gray")

Figure 2.2 (Right)

plot(x,dnorm(x,mu,sig),type="l",ylab=expression(paste(italic("p"),"(",italic("y"),")")),
xlab=expression(italic(y)),lwd=1)
abline(v=mu,col="gray")

Chapter 03

Figure 3.4 (Top Left)

a<-1
b<-1
n<-5
y<-1
theta<-seq(0,1,length.out=500)
plot(theta, dbeta(theta, a, b), 
     type="l",
     xlab=expression(theta),
     ylab = "probability",
     lwd = 2,
     col = "gray",
     ylim = c(0,3))
lines(theta, dbeta(theta, a+y, b+n-y), lwd = 2)
legend(0.6, 2.5, legend = c(
  expression(paste(italic("p"),"(",theta,")")),
  expression(paste(italic("p"),"(",theta,"|",italic("y"),")"))),
  lwd = c(2,2),
  col = c("gray", "black"),
  bty = "n"
  )

a=1
b=1
n=5
y=1
x<-seq(0,1,length=100)
plot(x = x,y=dbeta(x = x,shape1 = a,shape2 = b),type = "l",ylim = c(0,3))+
  points(x = x,y=dbeta(x = x,shape1 = a+y,shape2 = b+n-y),type = "l")

## integer(0)

Figure 3.4 (Top Right)

a<-3
b<-2
n<-5
y<-1
theta<-seq(0,1,length.out=500)
plot(theta, dbeta(theta, a, b), 
     type="l",
     xlab=expression(theta),
     ylab = "probability",
     lwd = 2,
     col = "gray",
     ylim = c(0,3))
lines(theta, dbeta(theta, a+y, b+n-y), lwd = 2)
legend(0.6, 2.5, legend = c(
  expression(paste(italic("p"),"(",theta,")")),
  expression(paste(italic("p"),"(",theta,"|",italic("y"),")"))),
  lwd = c(2,2),
  col = c("gray", "black"),
  bty = "n"
  )

a=3
b=2
n=5
y=1
x<-seq(0,1,length=100)
plot(x = x,y=dbeta(x = x,shape1 = a,shape2 = b),type = "l",ylim = c(0,3))+
  points(x = x,y=dbeta(x = x,shape1 = a+y,shape2 = b+n-y),type = "l")

## integer(0)

Figure 3.4 (Bottom Left)

a<-1
b<-1
n<-100
y<-20
theta<-seq(0,1,length.out=500)
plot(theta, dbeta(theta, a, b), 
     type="l",
     xlab=expression(theta),
     ylab = "probability",
     lwd = 2,
     col = "gray",
     ylim = c(0,10))
lines(theta, dbeta(theta, a+y, b+n-y), lwd = 2)
legend(0.6, 10, legend = c(
  expression(paste(italic("p"),"(",theta,")")),
  expression(paste(italic("p"),"(",theta,"|",italic("y"),")"))),
  lwd = c(2,2),
  col = c("gray", "black"),
  bty = "n"
)

a=1
b=1
n=100
y=20
x<-seq(0,1,length=100)
plot(x = x,y=dbeta(x = x,shape1 = a,shape2 = b),type = "l",ylim = c(0,10))+
  points(x = x,y=dbeta(x = x,shape1 = a+y,shape2 = b+n-y),type = "l")

## integer(0)

Figure 3.4 (Bottom Right)

a<-3
b<-2
n<-100
y<-20
theta<-seq(0,1,length.out=500)
plot(theta, dbeta(theta, a, b), 
     type="l",
     xlab=expression(theta),
     ylab = "probability",
     lwd = 2,
     col = "gray",
     ylim = c(0,10))
lines(theta, dbeta(theta, a+y, b+n-y), lwd = 2)
legend(0.6, 10, legend = c(
  expression(paste(italic("p"),"(",theta,")")),
  expression(paste(italic("p"),"(",theta,"|",italic("y"),")"))),
  lwd = c(2,2),
  col = c("gray", "black"),
  bty = "n"
)

a=3
b=2
n=100
y=20
x<-seq(0,1,length=100)
plot(x = x,y=dbeta(x = x,shape1 = a,shape2 = b),type = "l",ylim = c(0,10))+
  points(x = x,y=dbeta(x = x,shape1 = a+y,shape2 = b+n-y),type = "l")

## integer(0)

Figure 3.10

If the sampling model is p(y|θ)~pois(θ), and the prior p(θ)~gamma(a,b), then the posterior distribution p(θ|y)~gamma(a+∑y,b+n)

a<-2
b<-1
n1<-111
sy1<-217
n2<-44
sy2<-66

theta<-seq(0,5, length.out=500)

plot(theta, dgamma(theta,a+sy1,b+n1), type = "l", xlab = expression(theta), ylab = "probability", 
     lwd=2, col="gray", ylim = c(0,3.5))
lines(theta, dgamma(theta, a+sy2, b+n2), lwd = 2)
lines(theta, dgamma(theta, a, b), lty= 2, lwd = 2)
legend(3.5, 3, legend = c(expression(paste(italic("p"),"(",theta[1],"|",italic("y"),")")),
                          expression(paste(italic("p"),"(",theta[2],"|",italic("y"),")")),
                          expression(paste(italic("p"),"(",theta,")"))),
       lty = c(1,1,2), lwd = c(2,2,2), col = c("grey", "black", "black"), bty="n")