Chapter 16. BayesianMethods#
Figure 16.1. Posterior for Bernoulli Model#
x <- seq(0,1,.001)
alpha <- beta <- 5
b1 <- dbeta(x,alpha,beta)
b2 <- dbeta(x,8+alpha,2+beta)
b3 <- dbeta(x,24+alpha,6+beta)
leg1 <- "Prior"
leg2 <- "Posterior, n=10"
leg3 <- "Posterior, n=30"
m2 <- (8+alpha)/(10+alpha+beta)
bm2 <- dbeta(m2,8+alpha,2+beta)
m3 <- (24+alpha)/(30+alpha+beta)
bm3 <- dbeta(m3,24+alpha,6+beta)
plot(x,b1,type="l",lty=1,xaxs="i",yaxs="i",xlim=c(0,1),ylim=c(0,6),ylab="",xlab="",yaxt="n",bty="n",cex.axis=.75)
lines(x,b2)
lines(x,b3)
arrows(m2,bm2,m2,0,angle=20,length=.1)
arrows(m3,bm3,m3,0,angle=20,length=.1)
points(.8,0.05,pch=19,col="black",cex=.6)
text(.94,.5,"MLE",cex=.75)
arrows(.89,.45,.82,.1,angle=20,length=.1)
text(.15,1.5,leg1,cex=.75)
text(.43,3.5,leg2,cex=.75)
text(.5,5,leg3,cex=.75)
arrows(.2,1.4,.27,1.1,angle=20,length=.1)
arrows(.56,3.4,.6,3.2,angle=20,length=.1)
arrows(.63,4.85,.685,4.6,angle=20,length=.1)
Figure 16.2. Posterior for Mean in Normal Model#
x <- seq(-2,2,.001)
mu <- 0
lambda <- 3
alpha <- 2
beta <- 2
n2 <- 10
n3 <- 30
xbar <- 1
sigma2 <- 2
xbar1 <- mu
xbar2 <- (n2*xbar + lambda*mu)/(n2+lambda)
xbar3 <- (n3*xbar + lambda*mu)/(n3+lambda)
dof1 <- 2*alpha
dof2 <- n2 - 1 + 2*alpha
dof3 <- n3 - 1 + 2*alpha
s1 <- 2*beta
s2 <- n2*sigma2+2*beta
s3 <- n3*sigma2+2*beta
scale1 <- sqrt(beta/alpha/lambda)
scale2 <- sqrt(s2/(n2+lambda)/dof2)
scale3 <- sqrt(s3/(n3+lambda)/dof3)
b1 <- dt((x-xbar1)/scale1,dof1)/scale1
b2 <- dt((x-xbar2)/scale2,dof2)/scale2
b3 <- dt((x-xbar3)/scale3,dof3)/scale3
bm2 <- dt(0,dof2)/scale2
bm3 <- dt(0,dof3)/scale3
plot(x,b1,type="l",lty=1,xaxs="i",yaxs="i",xlim=c(-2,2),ylim=c(0,1.7),ylab="",xlab="",yaxt="n",bty="n",cex.axis=.75)
lines(x,b2)
lines(x,b3)
points(xbar,0.02,pch=19,col="black",cex=.6)
text(1.15,.25,"MLE",cex=.75)
arrows(1.15,.2,1.03,.05,angle=20,length=.1)
arrows(xbar2,bm2,xbar2,0,angle=20,length=.1)
arrows(xbar3,bm3,xbar3,0,angle=20,length=.1)
text(-1.5,.3,leg1,cex=.75)
text(-.5,1,leg2,cex=.75)
text(0,1.5,leg3,cex=.75)
arrows(-1.3,.28,-1,.2,angle=20,length=.1)
arrows(0.05,.98,.51,.89,angle=20,length=.1)
arrows(.53,1.49,.76,1.42,angle=20,length=.1)
Figure 16.3. Posterior for Precision in NormalModel#
x <- seq(0,2,.001)
b1 <- dchisq(x*s1,dof1)*s1
b2 <- dchisq(x*s2,dof2)*s2
b3 <- dchisq(x*s3,dof3)*s3
v2 <- dof2/s2
v3 <- dof3/s3
vm2 <- dchisq(v2*s2,dof2)*s2
vm3 <- dchisq(v3*s3,dof3)*s3
plot(x,b1,type="l",lty=1,xaxs="i",yaxs="i",xlim=c(0,2),ylim=c(0,3.5),ylab="",xlab="",yaxt="n",bty="n",cex.axis=.75)
lines(x,b2)
lines(x,b3)
arrows(v2,vm2,v2,0,angle=20,length=.1)
arrows(v3,vm3,v3,0,angle=20,length=.1)
points(1/sigma2,0.03,pch=19,col="black",cex=.6)
text(.37,.35,"MLE",cex=.75)
arrows(.40,.24,.48,.05,angle=20,length=.1)
text(1.75,.65,leg1,cex=.75)
text(1.2,1.1,leg2,cex=.75)
text(1.1,2,leg3,cex=.75)
arrows(1.7,.55,1.65,.3,angle=20,length=.1)
arrows(.93,1.1,.785,.9,angle=20,length=.1)
arrows(.83,2,.65,1.9,angle=20,length=.1)
Figure 16.4. Credible Set for NormalMean#
x <- seq(-2,2.5,.001)
mu <- 0
lambda <- 3
alpha <- 2
beta <- 2
n2 <- 10
xbar <- 1
sigma2 <- 2
xbar2 <- (n2*xbar + lambda*mu)/(n2+lambda)
dof <- n2 - 1
dof2 <- n2 - 1 + 2*alpha
s1 <- 2*beta
s2 <- n2*sigma2+2*beta
v2 <- dof2/s2
scale2 <- 1/sqrt((n2+lambda)*v2)
b2 <- dt((x-xbar2)/scale2,dof2)/scale2
bm2 <- dt(0,dof2)/scale2
q2 <- qt(.975,dof2)
cr1 <- xbar2 - q2/sqrt((n2+lambda)*v2)
cr2 <- xbar2 + q2/sqrt((n2+lambda)*v2)
br1 <- dt((cr1-xbar2)/scale2,dof2)/scale2
br2 <- dt((cr2-xbar2)/scale2,dof2)/scale2
q <- qt(.975,dof)
cL <- xbar - q*sqrt(sigma2/dof)
cU <- xbar + q*sqrt(sigma2/dof)
br <- -br1/2
plot(x,b2,type="l",lty=1,xaxs="i",yaxs="i",xlim=c(-.5,2.2),ylim=c(-0.07,1.1),ylab="",xlab="",yaxt="n",xaxt="n",bty="n")
points(xbar,0.01,pch=19,col="black",cex=.6)
text(0.5,.22,"Bayes",cex=.75)
text(1.27,.22,"MLE",cex=.75)
arrows(1.16,.19,1.02,.03,angle=20,length=.1)
arrows(0.63,.19,.75,.03,angle=20,length=.1)
points(xbar2,0.01,pch=19,col="black",cex=.6)
lines(c(cr1,cr2),c(br1,br2),lwd=1.5)
arrows(cr1,br1,cr1,0,angle=20,length=.1)
arrows(cr2,br2,cr2,0,angle=20,length=.1)
arrows(.35,.2,.2,br1,angle=20,length=.1)
arrows(1.36,.19,1.4,br,angle=20,length=.1)
lines(c(cL,cU),c(br,br),lwd=1.5)
arrows(cL,br,cL,0,angle=20,length=.1)
arrows(cU,br,cU,0,angle=20,length=.1)
lines(c(-1,3),c(0,0))
text(1.4,.72,expression(paste(pi,"(",theta,"|X)")),cex=.8)
Figure 16.5. Credible Set for Normal Precision#
x <- seq(0,2,.001)
b2 <- dchisq(x*s2,dof2)*s2
v2 <- dof2/s2
vm2 <- dchisq(v2*s2,dof2)*s2
mle <- 1/sigma2
c1 = qchisq(.025,dof)/(n2*sigma2)
c2 = qchisq(.975,dof)/(n2*sigma2)
br <- -.1
credible <- function(L){
F <- pchisq(L*s2,dof2)
Q <- qchisq(.95+F,dof2)/s2
f <- dchisq(L*s2,dof2)*s2 - dchisq(Q*s2,dof2)*s2
return(f)
}
LU <- qchisq(.05,dof2)/s2
cred <- uniroot(credible,c(0,LU))
L <- cred$root
U <- qchisq(.95+pchisq(L*s2,dof2),dof2)/s2
qL <- dchisq(L*s2,dof2)*s2
qU <- dchisq(U*s2,dof2)*s2
plot(x,b2,type="l",lty=1,xaxs="i",yaxs="i",xlim=c(0,1.5),ylim=c(-.12,2.1),ylab="",xlab="",xaxt="n",yaxt="n",bty="n")
lines(c(0,2),c(0,0))
points(mle,0.01,pch=19,col="black",cex=.6)
points(v2,0.01,pch=19,col="black",cex=.6)
text(.33,.4,"MLE",cex=.75)
arrows(.38,.34,.48,.03,angle=20,length=.1)
arrows(.27,.34,.24,br,angle=20,length=.1)
text(.7,.4,"Bayes",cex=.75)
arrows(.62,.34,.55,.03,angle=20,length=.1)
arrows(.78,.36,.82,qL,angle=20,length=.1)
lines(c(c1,c2),c(br,br),lwd=1.5)
arrows(c1,br,c1,0,angle=20,length=.1)
arrows(c2,br,c2,0,angle=20,length=.1)
lines(c(L,U),c(qL,qU),lwd=1.5)
arrows(L,qL,L,0,angle=20,length=.1)
arrows(U,qU,U,0,angle=20,length=.1)
text(0.8,1.4,expression(paste(pi,"(",theta,"|X)")),cex=.8)
