### R code from vignette source 'hierarchical.Rnw' ### Encoding: UTF-8 ################################################### ### code chunk number 1: hierarchical.Rnw:176-177 ################################################### library(mvtnorm) ################################################### ### code chunk number 2: hier1 ################################################### n=c(47,148,119,810,211,196,148,215,207,97,256,360) y=c(0,18,8,46,8,13,9,31,14,8,29,24) logmpostab=function(par) { -5/2*log(par[1]+par[2])+sum(lgamma(par[1]+par[2])+ lgamma(par[1]+y)+lgamma(par[2]+n-y)-lgamma(par[1])- lgamma(par[2])-lgamma(par[1]+par[2]+n)) } out=optim(logmpostab,par=c(1,1),control=list(fnscale=-1), hessian=TRUE) propab=rmvnorm(10000,mean=out$par,sigma=solve(-out$hessian)) gt0=apply(propab>0,FUN=all,MAR=1) propab=propab[gt0,] theta=matrix(ncol=12,nrow=nrow(propab)) for (i in 1:12) theta[,i]=rbeta(nrow(propab),shape1=(propab[,1]+y[i]), shape2=(propab[,2]+n[i]-y[i])) iw=exp(apply(propab,FUN=logmpostab,MAR=1)- dmvnorm(propab,mean=out$par,sigma=solve(-out$hessian),log=TRUE)) s=sample(nrow(propab),rep=T,prob=iw) theta=theta[s,]; boxplot(theta); points(y/n,col=2,pch=16) abline(h=sum(y)/sum(n)) ################################################### ### code chunk number 3: hier2 ################################################### betabinmwg=function(nMCMC,par.init,sigma.prop){ par=par.init post=matrix(ncol=(12+2),nrow=nMCMC) for (i in 1:nMCMC){ theta=rbeta(12,shape1=par[1]+y,shape2=par[2]+n-y) par.prop=rmvnorm(1,par,sigma.prop) if(all(par.prop>0)){ lnum=-5/2*log(sum(par.prop))+sum(dbeta(theta,shape1=par.prop[1], shape2=par.prop[2],log=TRUE)) lden=-5/2*log(sum(par))+sum(dbeta(theta,shape1=par[1], shape2=par[2],log=TRUE)) if(runif(1)