#example 8.1 Gelman book, fit a linear regression of radom measurements #explanatory variables:basement indicator and three counties indicators #beta_bayesian.regression(nsim=1000) make.indicators_function(x){ ux_unique(x) mat1_matrix(x,nrow=length(x),ncol=length(ux)) mat2_matrix(ux,nrow=length(x),ncol=length(ux),byrow=T) (mat1==mat2)*1} bayesian.regression_function(nsim=1000){ y.1_c(5.0, 13.0, 7.2,6.8,12.8,5.8,9.5,6.0,3.8,14.3,1.8,6.9,4.7,9.5) y.2_c(0.9,12.9,2.6,3.5,26.6,1.5,13.0,8.8,19.5,2.5,9.0,13.1,3.6,6.9) y.3_c(14.3,6.9,7.6,9.8,2.6,43.5,4.9,3.5,4.8,5.6,3.5,3.9,6.7) basement.1_c(1,1,1,1,1,0,1,1,1,0,1,1,1,1) basement.2_c(0,1,1,0,1,1,1,1,1,0,1,1,1,0) basement.3_c(1,0,1,0,1,1,1,1,1,1,1,1,1) counties_rep(1:3,c(length(y.1),length(y.2),length(y.3))) y_c(y.1,y.2,y.3) x_cbind(c(basement.1,basement.2,basement.3),make.indicators(counties)) ls.out_lsfit(x,log(y),intercept=F) lsd_ls.diag(ls.out) n_nrow(x) k_ncol(x) sigmasqr_rep(NA,nsim) beta_matrix(NA,nsim,k) for ( i in 1:nsim){ sigmasqr[i]_lsd$std.dev*(n-k)/rchisq(1,n-k) PP_ t(x)%*%x VV_solve(PP) VV_.5*(VV+t(VV)) beta[i,]_simulate.multnorm(as.numeric(ls.out$coef),sigmasqr[i]*VV) } output_exp(cbind(beta[,2],beta[,1]+beta[,2],beta[,3],beta[,1]+beta[,3], beta[,4],beta[,1]+beta[,4],sigmasqr)) for(i in 1:ncol(output)) print(round(quantile(output[,i],c(.25,.5,.75)),1)) return(beta) } makeboxplot_function(w=1,beta){ if (w == 1) postscript("/home/biostats/fdominic/course/radon.ps") coeff_c("BE, no basement","BE, basement","CC, no basement","CC, basement", "GC, no basement","GC, basement") boxplot(exp(beta[2,]),exp(beta[,1]+beta[,2]), exp(beta[3,]),exp(beta[,1]+beta[,3]), exp(beta[4,]),exp(beta[,1]+beta[4,]), style.bxp = "old",ylim=c(0,25),names=coeff,srt=45) #mtext(coeff, side=1,at=seq(1,80,length=6),srt=90,adj=.5 ,outer=F,line=1.5,cex=1) par(oma=c(0,0,0,0)) par(mfrow=c(1,1)) if (w == 1) dev.off() }