Part (a) - See solutions document.
Part (b) - The posterior distribution only involves the data through the sum of \(y_{i1}\) and \(y_{i2}\) for each \(i\). So we first reduce the data to these sums. Then we set up a grid approximation on \(\phi_1, \phi_2\) and sample from it.
# setting seed so I can reproduce these exact results
set.seed(1675493)
# set number of manufacturers and data
J <- 9
y <- c(2.3, 4.6,32.3, 4.0, 7.2, 25.5, 1.7, 11.9, 13.8)
xgridlow <- -150
xgridhigh <- 100
xgridint <- 50
ygridlow <- -100
ygridhigh <- 200
ygridint <- 25
# logpost is a function to evaluate the log posterior on the
# phi_1 = log(alpha/beta), phi_2 = log(beta) parameterization
logpost <- function(u, v)
{ a <- exp(u) * exp(v)
b <- exp(v)
x <- J * (lgamma(a + 2) - lgamma(a) + a*log(b)) + u - 0.5*v
for(i in (1:J)) {x <- x - (a+2)*log(b + y[i])}
x }
# set up grid
x4 <- c(xgridlow:xgridhigh)/xgridint
y4 <- c(ygridlow:ygridhigh)/ygridint
# evaluate grid
z <- outer(x4,y4,logpost)
# exponentiate log posterior (set maximum to be one)
z2 <- exp(z - max(z))
# contour plot
contour(x4,y4,z2,levels=c(0.05,0.15,0.25,0.35,0.45,0.55,0.65,0.75,0.85,0.95))
# sample from posterior distn
z2 <- z2/sum(z2)
p.x <- apply(z2,1,sum)
newx <- sample(x4,1000,replace=T,prob=p.x)
# next line takes sampled value 'newx' and identifies appropriate row 'xid' in the matrix z2
xid <- newx*xgridint - xgridlow + 1
newy <- rep(0,1000)
for (i in (1:1000)) newy[i] <- sample(y4,1,prob=z2[xid[i],]/sum(z2[xid[i],]))
The next step is to transform the \(\phi_1, \phi_2\) samples into samples of \(\alpha,\beta\) and then simulate \(\theta\)’s from their conditional posterior distribution.
# transform sample to get samples of alpha and beta
newa <- exp(newx)*exp(newy)
newb <- exp(newy)
# simulate thetas from conditional posterior
output <- matrix(0,1000,11)
for (j in (1:9)) {
output[,j] <- rgamma(1000,newa+2,newb+y[j])
}
output[,10] <- newa
output[,11] <- newb
dimnames(output)[[2]] <- c("theta1","theta2","theta3","theta4","theta5","theta6","theta7","theta8","theta9","alpha","beta")
describe <- function(x){c(n=length(x),mn=mean(x),sd=sqrt(var(x)),quantile(x,probs=c(0.025,.25,.50,.75,.975)))}
summaries <- apply(output,2,describe)
t(round(summaries,4))
## n mn sd 2.5% 25% 50% 75% 97.5%
## theta1 1000 0.5357 0.4030 0.1146 0.2631 0.4149 0.6917 1.6326
## theta2 1000 0.3807 0.2291 0.0983 0.2229 0.3208 0.4907 0.9656
## theta3 1000 0.1036 0.0571 0.0210 0.0586 0.0940 0.1382 0.2335
## theta4 1000 0.4127 0.2729 0.0901 0.2257 0.3462 0.5212 1.0724
## theta5 1000 0.2951 0.1700 0.0731 0.1749 0.2583 0.3758 0.7253
## theta6 1000 0.1285 0.0704 0.0273 0.0779 0.1168 0.1698 0.2926
## theta7 1000 0.6143 0.4786 0.1122 0.2791 0.4703 0.8229 1.8439
## theta8 1000 0.2155 0.1132 0.0508 0.1352 0.1936 0.2762 0.4871
## theta9 1000 0.1921 0.0966 0.0430 0.1217 0.1785 0.2493 0.4251
## alpha 1000 5.5906 26.3530 0.3263 0.8694 1.5220 3.0649 26.6169
## beta 1000 26.2577 142.8046 0.3393 1.8221 4.0552 11.0232 142.8907
Part (c) Find the probability that \(\theta_3 < \theta_7\):
sum(output[,3] < output[,7])/1000
## [1] 0.956
Part (d) To study the predictive distribution we generate a simulated lifetime for a lightbulb using each posterior draw of \(\theta_1\).
predbulb1 <- rgamma(1000,shape=1,rate=output[,1])
describe(predbulb1)
## n mn sd 2.5% 25%
## 1.000000e+03 3.158189e+00 5.228337e+00 3.504412e-02 5.853016e-01
## 50% 75% 97.5%
## 1.641156e+00 3.738150e+00 1.418875e+01
Part(e) To study the predictive distribution for a new manufacturer we simulate a new \(\theta\) and then the lifetime of a single bulb from that manufacturer. We repeat this 1000 times.
predtheta <- rgamma(1000,shape=newa,rate=newb)
predbulbn <- rgamma(1000,shape=1,rate=predtheta)
describe(predbulbn)
## n mn sd 2.5% 25%
## 1.000000e+03 2.266404e+03 4.727529e+04 5.678715e-02 8.014405e-01
## 50% 75% 97.5%
## 2.430258e+00 7.184158e+00 9.036538e+01
Part (a) - Set up a grid approximation on the interval [-1,10]
y<-c(1,0,0,0,1,0,1,0,0,1)
x<-c(1.1,-2.3,-1.0,-5.3,2.5,-5.2,2.6,-6.0,2.6,0.3)
#
# create grid approximation
#
logpost<-function(beta){
p<-exp(beta*x)/(1+exp(beta*x))
sum(y*log(p) + (1-y)*log(1-p))}
grid<-c(-100:1000)/100
dens<-rep(0,1101)
for (i in (1:1101)) dens[i]<-logpost(grid[i])
dens<- exp(dens-max(dens))
dens<-dens/sum(dens)
plot(grid,dens,type="l",col="black")
Part (b) - Find mode and 2nd derivative at mode using Newton’s method
der1<-function(beta){
p<-exp(beta*x)/(1+exp(beta*x))
sum(x*(y-p))}
der2<-function(beta){
p<-exp(beta*x)/(1+exp(beta*x))
sum(-x*x*p*(1-p))}
# set starting value for Newton's method
beta<-1
# run for ten Newton steps
for (j in (1:10)) beta<-beta-der1(beta)/der2(beta)
cat("posterior mode:",beta," 2nd deriv:",der2(beta))
## posterior mode: 0.7075039 2nd deriv: -5.389671
Part (c) - Normal approximation is centered at the posterior mode with variance equal to 1/-der2(beta)
mn <- beta
var <- -1/der2(beta)
#plot normal density and grid approx
densnorm <- dnorm(grid,mn,sqrt(var))
densnorm <- densnorm/sum(densnorm)
matplot(grid,cbind(dens,densnorm),type="l",col=c("black","red"))
text(1.9,.009,"norm approx",col="red")
text(2.9,.001,"posterior",col="black")