#Example:analysis o biossay expriment #Gelman and Rubin book, pag 82 #DD_input() input_function(LL=100){ biossay _ data.frame(doses=c(-.863,-.296,-.053,.727), rats=c(5,5,5,5),deaths=c(0.5,1,3,4.5),freq=c(0.5,1,3,4.5)/5) DD _ data.frame(y=log((biossay$freq)/(1-biossay$freq)),x=biossay$doses) estimates _ lm(y~x,data=DD) alpha.hat _ summary(estimates)$coeff[1,1] std.alpha _ summary(estimates)$coeff[1,2] beta.hat _ summary(estimates)$coeff[2,1] std.beta _ summary(estimates)$coeff[2,2] alpha _ seq(-2,2,length=LL) beta _ seq(-5,10,length=LL) return(biossay,estimates,alpha,beta) } log.post_function(alpha,beta,data=DD$biossay){ ldens _ 0 for (i in 1:length(data$doses)){ theta _ 1/(1+exp(-alpha-beta*data$doses[i])) ldens _ ldens + data$deaths[i]*log(theta)+(data$rats[i]-data$deaths[i])*log(1-theta) } ldens } plot.joint.post_function(w=0,data,draws){ if (w == 1) postscript("/home/biostats/fdominic/course/jposterior.ps") contours _ seq(.05,.95,.1) logdens _outer(DD$alpha,DD$beta,log.post,data) dens_exp(logdens-max(logdens)) contour(DD$alpha,DD$beta,dens,levels=contours,xlab="alpha",ylab="beta",labex=0,cex=1) points(draws$aaa,draws$bbb) mtext("Posterior density",3,line=1,cex=1.2) par(oma=c(0,0,0,0)) par(mfrow=c(1,1)) if (w == 1) dev.off() } grid.value_function(data=DD$biossay,alpha,beta){ ll _ length(alpha) PP _ matrix(NA,ll,ll) for(i in 1:ll){ for(j in 1:ll){ PP[i,j]_exp(log.post(alpha[i],beta[j],data)) }} ccc _ sum(PP) PP _ PP/ccc return(PP) } sampling_function(M=100,PP,alpha,beta){ alpha.mar_apply(PP,1,sum) alpha.cdf <- cumsum(alpha.mar) post.alpha <- post.beta <- rep(0,M) for( m in 1:M){ uuu_runif(1,0,1) Fhat.alpha <- max( alpha.cdf[ alpha.cdf <= uuu]) post.alpha[m] <- alpha[(1:length(alpha.cdf))[alpha.cdf == Fhat.alpha]] junk <- length(alpha[alpha <= post.alpha[m]]) PP[junk, ] <- PP[junk,]/sum(PP[junk,]) beta.cond.cdf <- cumsum(PP[junk,]) uuu_runif(1,0,1) Fhat.beta <- max( beta.cond.cdf[ beta.cond.cdf < uuu]) post.beta[m] <- beta[(1:length(beta.cond.cdf))[beta.cond.cdf == Fhat.beta]] } return(post.alpha,post.beta) }