Skip to content

# jet08013/BDA

Fetching contributors…
Cannot retrieve contributors at this time
37 lines (35 sloc) 1.52 KB
 post.a <- function(mu,sd,y){ ldens <- 0 for (i in 1:length(y)) ldens <- ldens + log(dnorm(y[i],mu,sd)) ldens} post.b <- function(mu,sd,y){ ldens <- 0 for (i in 1:length(y)) ldens <- ldens + log(pnorm(y[i]+0.5,mu,sd) - pnorm(y[i]-0.5,mu,sd)) ldens} summ <- function(x){c(mean(x),sqrt(var(x)), quantile(x, c(.025,.25,.5,.75,.975)))} nsim <- 2000 y <- c(10,10,12,11,9) n <- length(y) ybar <- mean(y) s2 <- sum((y-mean(y))^2)/(n-1) mugrid <- seq(3,18,length=200) logsdgrid <- seq(-2,4,length=200) contours <- c(.0001,.001,.01,seq(.05,.95,.05)) logdens <- outer (mugrid, exp(logsdgrid), post.a, y) dens <- exp(logdens - max(logdens)) contour (mugrid, logsdgrid, dens, levels=contours, xlab="mu", ylab="log sigma", labex=0, cex=2) mtext ("Posterior density, ignoring rounding", 3) sd <- sqrt((n-1)*s2/rchisq(nsim,4)) mu <- rnorm(nsim,ybar,sd/sqrt(n)) print (rbind (summ(mu),summ(sd))) logdens <- outer (mugrid, exp(logsdgrid), post.b, y) dens <- exp(logdens - max(logdens)) contour (mugrid, logsdgrid, dens, levels=contours,xlab="mu", ylab="log sigma", labex=0, cex=2) mtext ("Posterior density, accounting for rounding", cex=2, 3) dens.mu <- apply(dens,1,sum) muindex <- sample (1:length(mugrid), nsim, replace=T, prob=dens.mu) mu <- mugrid[muindex] sd <- rep (NA,nsim) for (i in (1:nsim)) sd[i] <- exp (sample(logsdgrid, 1, prob=dens[muindex[i],])) print (rbind (summ(mu),summ(sd)))
You can’t perform that action at this time.