# CHAPTER 4: R COMMANDS # 19/12/2006 ################################################################################ # GAUSSIAN RANDOM WALK MH FOR THE STANDARD NORMAL DISTRIBUTION test=hm1(10000,1,0.001) par(mfrow=c(3,1)) plot(test,type="l",xlab="Iterations",ylab="MH chain") hist(test[8001:10000],nclass=50,prob=TRUE,main="",xlab="x") curve(dnorm,-3,3,add=TRUE) acf(test,lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hm1.eps") test=hm1(10000,1,1000) par(mfrow=c(3,1)) plot(test,type="l",xlab="Iterations",ylab="MH chain") hist(test[8001:10000],nclass=50,prob=TRUE,main="",xlab="x") curve(dnorm,-3,3,add=TRUE) acf(test,lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hm2.eps") test=hm1(10000,1,1) par(mfrow=c(3,1)) plot(test,type="l",xlab="Iterations",ylab="MH chain") hist(test[8001:10000],nclass=50,prob=TRUE,main="",xlab="x") curve(dnorm,-3,3,add=TRUE) acf(test,lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hm3.eps") ################################################################################ # bank DATASET IMPLEMENTATION bank=read.table(file="bank") bank=as.matrix(bank) y=bank[,5] X=bank[,1:4] par(mfrow=c(1,2)) plot(y~X[,4],xlab="Bottom margin width (mm)",ylab="Status") boxplot(X[,4]~y,ylab="Bottom margin width (mm)",xlab="Status") dev.copy2eps(file="bottom.eps") dev.off() model1=glm(y~-1+X,family=binomial(link="probit")) summary(model1) model2=glm(y~-1+X,family=binomial(link="logit")) summary(model2) source("#4.R") ################################################################################ # bank DATASET: PROBIT MODEL # bank DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER FLAT PRIOR flatprobit=hmflatprobit(10000,y,X,1) par(mfrow=c(4,3),mar=1+c(1.5,1.5,1.5,1.5)) plot(flatprobit[,1],type="l",xlab="Iterations",ylab=expression(beta[1])) hist(flatprobit[1001:10000,1],nclass=50,prob=TRUE,main="",xlab=expression(beta[1])) acf(flatprobit[1001:10000,1],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatprobit[,2],type="l",xlab="Iterations",ylab=expression(beta[2])) hist(flatprobit[1001:10000,2],nclass=50,prob=TRUE,main="",xlab=expression(beta[2])) acf(flatprobit[1001:10000,2],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatprobit[,3],type="l",xlab="Iterations",ylab=expression(beta[3])) hist(flatprobit[1001:10000,3],nclass=50,prob=TRUE,main="",xlab=expression(beta[3])) acf(flatprobit[1001:10000,3],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatprobit[,4],type="l",xlab="Iterations",ylab=expression(beta[4])) hist(flatprobit[1001:10000,4],nclass=50,prob=TRUE,main="",xlab=expression(beta[4])) acf(flatprobit[1001:10000,4],lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hmflatprobit.eps") mean(flatprobit[1001:10000,1]) mean(flatprobit[1001:10000,2]) mean(flatprobit[1001:10000,3]) mean(flatprobit[1001:10000,4]) pnorm(-1.2193*214.9+0.9540*130.1+0.9795*129.9+1.1481*9.5) # bank DATASET: GIBBS SAMPLING UNDER FLAT PRIOR gibbsflatprobit=gibbsprobit(10000,y,X) par(mfrow=c(4,3),mar=1+c(1.5,1.5,1.5,1.5)) plot(gibbsflatprobit[,1],type="l",xlab="Iterations",ylab=expression(beta[1])) hist(gibbsflatprobit[1001:10000,1],nclass=50,prob=TRUE,main="",xlab=expression(beta[1])) acf(gibbsflatprobit[1001:10000,1],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(gibbsflatprobit[,2],type="l",xlab="Iterations",ylab=expression(beta[2])) hist(gibbsflatprobit[1001:10000,2],nclass=50,prob=TRUE,main="",xlab=expression(beta[2])) acf(gibbsflatprobit[1001:10000,2],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(gibbsflatprobit[,3],type="l",xlab="Iterations",ylab=expression(beta[3])) hist(gibbsflatprobit[1001:10000,3],nclass=50,prob=TRUE,main="",xlab=expression(beta[3])) acf(gibbsflatprobit[1001:10000,3],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(gibbsflatprobit[,4],type="l",xlab="Iterations",ylab=expression(beta[4])) hist(gibbsflatprobit[1001:10000,4],nclass=50,prob=TRUE,main="",xlab=expression(beta[4])) acf(gibbsflatprobit[1001:10000,4],lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="gibbsflatprobit.eps") mean(gibbsflatprobit[1001:10000,1]) mean(gibbsflatprobit[1001:10000,2]) mean(gibbsflatprobit[1001:10000,3]) mean(gibbsflatprobit[1001:10000,4]) # bank DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER NON-INFORMATIVE PRIOR noinfprobit=hmnoinfprobit(10000,y,X,1) par(mfrow=c(4,3),mar=1+c(1.5,1.5,1.5,1.5)) plot(noinfprobit[,1],type="l",xlab="Iterations",ylab=expression(beta[1])) hist(noinfprobit[1001:10000,1],nclass=50,prob=TRUE,main="",xlab=expression(beta[1])) acf(noinfprobit[1001:10000,1],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinfprobit[,2],type="l",xlab="Iterations",ylab=expression(beta[2])) hist(noinfprobit[1001:10000,2],nclass=50,prob=TRUE,main="",xlab=expression(beta[2])) acf(noinfprobit[1001:10000,2],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinfprobit[,3],type="l",xlab="Iterations",ylab=expression(beta[3])) hist(noinfprobit[1001:10000,3],nclass=50,prob=TRUE,main="",xlab=expression(beta[3])) acf(noinfprobit[1001:10000,3],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinfprobit[,4],type="l",xlab="Iterations",ylab=expression(beta[4])) hist(noinfprobit[1001:10000,4],nclass=50,prob=TRUE,main="",xlab=expression(beta[4])) acf(noinfprobit[1001:10000,4],lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hmnoinfprobit.eps") mean(noinfprobit[1001:10000,1]) mean(noinfprobit[1001:10000,2]) mean(noinfprobit[1001:10000,3]) mean(noinfprobit[1001:10000,4]) # bank DATASET: BAYES FACTORS CALCULATION library(mnormt) mkprob=apply(noinfprobit,2,mean) vkprob=var(noinfprobit) simk=rmnorm(100000,mkprob,2*vkprob) usk=probitnoinflpost(simk,y,X)-dmnorm(simk,mkprob,2*vkprob,log=TRUE) noinfprobit0=hmnoinfprobit(10000,y,X[,3:4],1) mk0=apply(noinfprobit0,2,mean) vk0=var(noinfprobit0) simk0=rmnorm(100000,mk0,2*vk0) usk0=probitnoinflpost(simk0,y,X[,3:4])-dmnorm(simk0,mk0,2*vk0,log=TRUE) bf0probit=mean(exp(usk))/mean(exp(usk0)) rm(mk0,vk0,simk0,usk0) rm(noinfprobit0) noinfprobit1=hmnoinfprobit(10000,y,X[,2:4],1) mk1=apply(noinfprobit1,2,mean) vk1=var(noinfprobit1) simk1=rmnorm(100000,mk1,2*vk1) usk1=probitnoinflpost(simk1,y,X[,2:4])-dmnorm(simk1,mk1,2*vk1,log=TRUE) bf1probit=mean(exp(usk))/mean(exp(usk1)) rm(mk1,vk1,simk1,usk1) rm(noinfprobit1) noinfprobit2=hmnoinfprobit(10000,y,cbind(X[,1],X[,3:4]),1) mk2=apply(noinfprobit2,2,mean) vk2=var(noinfprobit2) simk2=rmnorm(100000,mk2,2*vk2) usk2=probitnoinflpost(simk2,y,cbind(X[,1],X[,3:4]))-dmnorm(simk2,mk2,2*vk2,log=TRUE) bf2probit=mean(exp(usk))/mean(exp(usk2)) rm(mk2,vk2,simk2,usk2) rm(noinfprobit2) noinfprobit3=hmnoinfprobit(10000,y,cbind(X[,1:2],X[,4]),1) mk3=apply(noinfprobit3,2,mean) vk3=var(noinfprobit3) simk3=rmnorm(100000,mk3,2*vk3) usk3=probitnoinflpost(simk3,y,cbind(X[,1:2],X[,4]))-dmnorm(simk3,mk3,2*vk3,log=TRUE) bf3probit=mean(exp(usk))/mean(exp(usk3)) rm(mk3,vk3,simk3,usk3) rm(noinfprobit3) noinfprobit4=hmnoinfprobit(10000,y,X[,1:3],1) mk4=apply(noinfprobit4,2,mean) vk4=var(noinfprobit4) simk4=rmnorm(100000,mk4,2*vk4) usk4=probitnoinflpost(simk4,y,X[,1:3])-dmnorm(simk4,mk4,2*vk4,log=TRUE) bf4probit=mean(exp(usk))/mean(exp(usk4)) rm(mk4,vk4,simk4,usk4) rm(noinfprobit4) rm(simk,usk) log10(bf0probit) mkprob diag(vkprob) log10(c(bf1probit,bf2probit,bf3probit,bf4probit)) ################################################################################ # bank DATASET: LOGIT MODEL # bank DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER FLAT PRIOR flatlogit=hmflatlogit(10000,y,X,1) par(mfrow=c(4,3),mar=1+c(1.5,1.5,1.5,1.5)) plot(flatlogit[,1],type="l",xlab="Iterations",ylab=expression(beta[1])) hist(flatlogit[1001:10000,1],nclass=50,prob=TRUE,main="",xlab=expression(beta[1])) acf(flatlogit[1001:10000,1],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatlogit[,2],type="l",xlab="Iterations",ylab=expression(beta[2])) hist(flatlogit[1001:10000,2],nclass=50,prob=TRUE,main="",xlab=expression(beta[2])) acf(flatlogit[1001:10000,2],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatlogit[,3],type="l",xlab="Iterations",ylab=expression(beta[3])) hist(flatlogit[1001:10000,3],nclass=50,prob=TRUE,main="",xlab=expression(beta[3])) acf(flatlogit[1001:10000,3],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(flatlogit[,4],type="l",xlab="Iterations",ylab=expression(beta[4])) hist(flatlogit[1001:10000,4],nclass=50,prob=TRUE,main="",xlab=expression(beta[4])) acf(flatlogit[1001:10000,4],lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hmflatlogit.eps") mean(flatlogit[1001:10000,1]) mean(flatlogit[1001:10000,2]) mean(flatlogit[1001:10000,3]) mean(flatlogit[1001:10000,4]) exp(-2.5888*214.9+1.9967*130.1+2.1260*129.9+2.1879*9.5)/(1+exp(-2.5888*214.9+1.9967*130.1+2.1260*129.9+2.1879*9.5)) # bank DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER NON-INFORMATIVE PRIOR noinflogit=hmnoinflogit(10000,y,X,1) par(mfrow=c(4,3),mar=1+c(1.5,1.5,1.5,1.5)) plot(noinflogit[,1],type="l",xlab="Iterations",ylab=expression(beta[1])) hist(noinflogit[1001:10000,1],nclass=50,prob=TRUE,main="",xlab=expression(beta[1])) acf(noinflogit[1001:10000,1],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinflogit[,2],type="l",xlab="Iterations",ylab=expression(beta[2])) hist(noinflogit[1001:10000,2],nclass=50,prob=TRUE,main="",xlab=expression(beta[2])) acf(noinflogit[1001:10000,2],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinflogit[,3],type="l",xlab="Iterations",ylab=expression(beta[3])) hist(noinflogit[1001:10000,3],nclass=50,prob=TRUE,main="",xlab=expression(beta[3])) acf(noinflogit[1001:10000,3],lag=1000,main="",ylab="Autocorrelation",ci=F) plot(noinflogit[,4],type="l",xlab="Iterations",ylab=expression(beta[4])) hist(noinflogit[1001:10000,4],nclass=50,prob=TRUE,main="",xlab=expression(beta[4])) acf(noinflogit[1001:10000,4],lag=1000,main="",ylab="Autocorrelation",ci=F) dev.copy2eps(file="hmnoinflogit.eps") mean(noinflogit[1001:10000,1]) mean(noinflogit[1001:10000,2]) mean(noinflogit[1001:10000,3]) mean(noinflogit[1001:10000,4]) # bank DATASET: BAYES FACTORS CALCULATION library(mnormt) mklog=apply(noinflogit,2,mean) vklog=var(noinflogit) simk=rmnorm(100000,mklog,2*vklog) usk=logitnoinflpost(simk,y,X[,1:4])-dmnorm(simk,mklog,2*vklog,log=TRUE) noinflogit0=hmnoinflogit(10000,y,X[,3:4],1) mk0=apply(noinflogit0,2,mean) vk0=var(noinflogit0) simk0=rmnorm(100000,mk0,2*vk0) usk0=logitnoinflpost(simk0,y,X[,3:4])-dmnorm(simk0,mk0,2*vk0,log=TRUE) bf0logit=mean(exp(usk))/mean(exp(usk0)) rm(mk0,vk0,simk0,usk0) rm(noinflogit0) noinflogit1=hmnoinflogit(10000,y,X[,2:4],1) mk1=apply(noinflogit1,2,mean) vk1=var(noinflogit1) simk1=rmnorm(100000,mk1,2*vk1) usk1=logitnoinflpost(simk1,y,X[,2:4])-dmnorm(simk1,mk1,2*vk1,log=TRUE) bf1logit=mean(exp(usk))/mean(exp(usk1)) rm(mk1,vk1,simk1,usk1) rm(noinflogit1) noinflogit2=hmnoinflogit(10000,y,cbind(X[,1],X[,3:4]),1) mk2=apply(noinflogit2,2,mean) vk2=var(noinflogit2) simk2=rmnorm(100000,mk2,2*vk2) usk2=logitnoinflpost(simk2,y,cbind(X[,1],X[,3:4]))-dmnorm(simk2,mk2,2*vk2,log=TRUE) bf2logit=mean(exp(usk))/mean(exp(usk2)) rm(mk2,vk2,simk2,usk2) rm(noinflogit2) noinflogit3=hmnoinflogit(10000,y,cbind(X[,1:2],X[,4]),1) mk3=apply(noinflogit3,2,mean) vk3=var(noinflogit3) simk3=rmnorm(100000,mk3,2*vk3) usk3=logitnoinflpost(simk3,y,cbind(X[,1:2],X[,4]))-dmnorm(simk3,mk3,2*vk3,log=TRUE) bf3logit=mean(exp(usk))/mean(exp(usk3)) rm(mk3,vk3,simk3,usk3) rm(noinflogit3) noinflogit4=hmnoinfprobit(10000,y,X[,1:3],1) mk4=apply(noinflogit4,2,mean) vk4=var(noinflogit4) simk4=rmnorm(100000,mk4,2*vk4) usk4=logitnoinflpost(simk4,y,X[,1:3])-dmnorm(simk4,mk4,2*vk4,log=TRUE) bf4logit=mean(exp(usk))/mean(exp(usk4)) rm(mk4,vk4,simk4,usk4) rm(noinflogit4) rm(simk,usk) log10(bf0logit) mklog diag(vklog) log10(c(bf1logit,bf2logit,bf3logit,bf4logit)) ################################################################################ # airquality DATASET IMPLEMENTATION airquality=read.table(file="airquality") airqual=na.omit(airquality) ozone=cut(airqual\$Ozone,c(min(airqual\$Ozone),median(airqual\$Ozone),max(airqual\$Ozone)),include.lowest=TRUE) month=as.factor(airqual\$Month) tempe=cut(airqual\$Temp,c(min(airqual\$Temp),median(airqual\$Temp),max(airqual\$Temp)),include.lowest=TRUE) counts=table(ozone,tempe,month) is.array(counts) ftable(ozone,tempe,month) counts=as.vector(counts) ozo=gl(2,1,20) temp=gl(2,2,20) mon=gl(5,4,20) model1=glm(counts~ozo+temp+mon+ozo*temp+ozo*mon+temp*mon,family=poisson()) anova(model1) summary(model1) x1=rep(1,20) lulu=rep(0,20) x2=lulu x2[ozo==2]=1 x3=lulu x3[temp==2]=1 x4=lulu x4[mon==2]=1 x5=lulu x5[mon==3]=1 x6=lulu x6[mon==4]=1 x7=lulu x7[mon==5]=1 x8=lulu x8[ozo==2 & temp==2]=1 x9=lulu x9[ozo==2 & mon==2]=1 x10=lulu x10[ozo==2 & mon==3]=1 x11=lulu x11[ozo==2 & mon==4]=1 x12=lulu x12[ozo==2 & mon==5]=1 x13=lulu x13[temp==2 & mon==2]=1 x14=lulu x14[temp==2 & mon==3]=1 x15=lulu x15[temp==2 & mon==4]=1 x16=lulu x16[temp==2 & mon==5]=1 X=cbind(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16) model2=glm(counts~-1+X,family=poisson()) summary(model2) ################################################################################ # airquality DATASET: LOG-LINEAR MODEL # airquality DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER FLAT PRIOR flatloglin=hmflatloglin(10000,counts,X,0.5) par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) plot(flatloglin[,i],type="l",ylab="",xlab="Iterations") dev.copy2eps(file="trajflatloglin.eps") par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) hist(flatloglin[1001:10000,i],nclass=30,main="",ylab="",xlab="") dev.copy2eps(file="histflatloglin.eps") par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) acf(flatloglin[1001:10000,i],main="") dev.copy2eps(file="autoflatloglin.eps") apply(flatloglin[1001:10000,],2,mean) # airquality DATASET: GAUSSIAN RANDOM WALK MH SAMPLER UNDER NON-INFORMATIVE PRIOR noinfloglin=hmnoinfloglin(10000,counts,X,0.5) par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) plot(noinfloglin[,i],type="l",ylab="",xlab="Iterations") dev.copy2eps(file="trajnoinfloglin.eps") par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) hist(noinfloglin[1001:10000,i],nclass=30,main="",ylab="",xlab="") dev.copy2eps(file="histnoinfloglin.eps") par(mfrow=c(4,4),mar=1+c(1.5,1.5,1.5,1.5)) for (i in 1:16) acf(noinfloglin[1001:10000,i],main="") dev.copy2eps(file="autonoinfloglin.eps") apply(noinfloglin[1001:10000,],2,mean) # airquality DATASET: BAYES FACTORS CALCULATION library(mnormt) mklog=apply(noinfloglin,2,mean) vklog=var(noinfloglin) simk=rmnorm(100000,mklog,2*vklog) usk=loglinnoinflpost(simk,counts,X)-dmnorm(simk,mklog,2*vklog,log=TRUE) noinfloglin1=hmnoinfloglin(10000,counts,cbind(X[,1:7],X[,9:16]),0.5) mk1=apply(noinfloglin1,2,mean) vk1=var(noinfloglin1) simk1=rmnorm(100000,mk1,2*vk1) usk1=loglinnoinflpost(simk1,counts,cbind(X[,1:7],X[,9:16]))-dmnorm(simk1,mk1,2*vk1,log=TRUE) bf1loglin=mean(exp(usk))/mean(exp(usk1)) rm(mk1,vk1,simk1,usk1) rm(noinfloglin1) noinfloglin2=hmnoinfloglin(10000,counts,cbind(X[,1:8],X[,13:16]),0.5) mk2=apply(noinfloglin2,2,mean) vk2=var(noinfloglin2) simk2=rmnorm(100000,mk2,2*vk2) usk2=loglinnoinflpost(simk2,counts,cbind(X[,1:8],X[,13:16]))-dmnorm(simk2,mk2,2*vk2,log=TRUE) bf2loglin=mean(exp(usk))/mean(exp(usk2)) rm(mk2,vk2,simk2,usk2) rm(noinfloglin2) noinfloglin3=hmnoinfloglin(10000,counts,X[,1:12],0.5) mk3=apply(noinfloglin3,2,mean) vk3=var(noinfloglin3) simk3=rmnorm(100000,mk3,2*vk3) usk3=loglinnoinflpost(simk3,counts,X[,1:12])-dmnorm(simk3,mk3,2*vk3,log=TRUE) bf3loglin=mean(exp(usk))/mean(exp(usk3)) rm(mk3,vk3,simk3,usk3) rm(noinfloglin3) mklog diag(vklog) log10(c(bf1loglin,bf2loglin,bf3loglin)) ################################################################################