Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
ISYE 6420 Homework
Problem 1 (a) The original data summarizes the number of breakdowns at each voltage level and we can convert the outcome to be binary. As shown in our lecture note, we have yi = { 1 if zi > 0, 0 if zi ≤ 0, zi|β ∼ N(β0 + β1xi, 1), θi = p(yi = 1) = p(zi > 0) = Φ(β0 + β1xi). Assuming noninformative prior p(β) ∝ 1, the full conditional distributions for zi and β can be derived as follows: β|z ∼ N((X ′X)−1X ′z, (X ′X)−1), p(zi|β, y) = { p(zi|β, zi > 0) if yi = 1 p(zi|β, zi ≤ 0) if yi = 0. The Gibbs sampler is implemented by the following R code: ################################### library(numDeriv) library(MASS) library(coda) library(mcmc) library(adaptMCMC) library(mined) #Problem#1 voltage=c(1065,1071,1075,1083,1089,1094,1100,1107,1111,1120,1128,1135) trials=rep(100,12) breakdown=c(2,3,5,11,10,21,29,48,56,88,98,99) voltagec=voltage-mean(voltage)#center the voltage a=glm(cbind(breakdown,trials-breakdown)~voltagec,family=binomial(link="probit")) summary(a) m=10000 beta=matrix(0,nrow=m,ncol=2) z=numeric(100) 1 zbar=matrix(0,nrow=m,ncol=12) beta[1,]=a$coef library(MASS) X=as.matrix(cbind(1,voltagec)) Sinv=solve(t(X)%*%X) library(truncnorm) for(i in 2:m) { for(j in 1:12) { mu=sum(X[j,]*beta[i-1,]) z1=rtruncnorm(breakdown[j],a=0,b=Inf,mean=mu,sd=1) z0=rtruncnorm(100-breakdown[j],a=-Inf,b=0,mean=mu,sd=1) z=c(z1,z0) zbar[i,j]=mean(z) } beta[i,]=mvrnorm(1,Sinv%*%t(X)%*%zbar[i,],1/100*Sinv) } library(coda) effectiveSize(beta) cumuplot(mcmc(beta,start=1,end=m,thin=10)) summary(beta) #1(a) par(mfrow=c(1,2)) plot(density(beta[,1]),main=expression(beta[0]),xlab="") plot(density(beta[,2]),main=expression(beta[1]),xlab="") par(mfrow=c(1,1)) #1(b) ld50=-beta[,1]/beta[,2]+mean(voltage) plot(density(ld50),main="LD50") ################################### For m = 20, 000, the posterior densities for β0 and β1 are shown in Figure 1. (b) Since θ = 0.5 ⇐⇒ β0 + β1x = 0 ⇐⇒ x = −β0β1 , the LD50 point corresponds to − β0 β1 . The distribution of LD50 is plotted in Figure 2. Problem 2 2 Figure 1: Posterior densities for β0 and β1. −0.7 −0.4 0 2 4 6 8 β0 D en si ty 0.055 0.070 0 20 40 60 80 10 0 β1 D en si ty If we fit a logistic regression model to the above data, the posterior distribution of β is p(β|y) ∝ 12∏ i=1 θyii (1− θi)100−yi = 12∏ i=1 ( ex T i β 1 + ex T i β )yi( 1 1 + ex T i β )100−yi . Note that yi here is the number of breakdowns instead of 0-1. Both mcmc and adaptMCMC packages use log of unnormalized posterior as input, and we can obtain the following results via the R code provided. ################################################################################### a=glm(cbind(breakdown,trials-breakdown)~voltagec,family=binomial(link="logit")) summary(a) logh=function(beta){ val=sum(breakdown*(X%*%beta))-100*sum(log(1+exp(X%*%beta))) return(val) } nlogh=function(theta) -logh(theta) betahat=optim(a$coef,nlogh)$par Sigmahat=solve(hessian(nlogh,betahat)) 3 Figure 2: Distribution of LD50 using Gibbs sampling. 1102 1103 1104 1105 1106 1107 1108 1109 0. 0 0. 1 0. 2 0. 3 0. 4 0. 5 Density plot for LD50 N = 18001 Bandwidth = 0.09821 D en si ty #mcmc eig=eigen(Sigmahat) sqrtSigma=eig$vec%*%diag(sqrt(eig$val))%*%t(eig$vec) s=2.4/sqrt(2)*sqrtSigma m=10000 out=metrop(logh,initial=betahat,nbatch=m,scale=s) out$accept beta=out$batch effectiveSize(beta) cumuplot(mcmc(beta,start=100,end=m,thin=10)) 1-rejectionRate(mcmc(beta)) ld50=-beta[,1]/beta[,2]+mean(voltage) plot(density(ld50),main="LD50") #adaptMCMC s=(2.4/sqrt(2))^2*Sigmahat#specify covariance matrix of the jumping distribution out=MCMC(logh,n=m,init=betahat,scale=s,adapt=T,acc.rate=.4) out$accept beta=out$samples effectiveSize(beta) 4 cumuplot(mcmc(beta,thin=10)) 1-rejectionRate(mcmc(beta)) ld50=-beta[,1]/beta[,2]+mean(voltage) plot(density(ld50),main="LD50") #mined #first find approximate ranges for the parameters summary(a) a$coef[1]+3*c(-1,1)*.0945 a$coef[2]+3*c(-1,1)*.0064 l1=-1.2 u1=-.6 l2=.09 u2=.14 X=as.matrix(cbind(1,voltagec)) logh=function(x){ beta=x beta[1]=l1+(u1-l1)*x[1] beta[2]=l2+(u2-l2)*x[2] val=sum(breakdown*(X%*%beta))-100*sum(log(1+exp(X%*%beta))) return(val) } ini=Lattice(109,2) D=mined(ini,logh)$points beta=matrix(0,nrow=109,ncol=2) beta[,1]=l1+(u1-l1)*D[,1] beta[,2]=l2+(u2-l2)*D[,2] ld50=-beta[,1]/beta[,2]+mean(voltage) plot(density(ld50),main="LD50")