HW 2 Question 3

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

HW 2 Question 4

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")