## Monday, September 16, 2013

### Visualizing a normal distribution as outcomes from a spinner

Mike Meredith has a new blog post showing an animation that converts a spinner to a normal distribution. It shows hows segments of the spinner map to regions under the normal distribution. Give it a look!

(Mike Meredith is also the fellow who created the BEST package for R. See this previous post about that!)

## Sunday, September 15, 2013

### JAGS users: runjags has been updated

Image credit: http://www.potracksgalore.com/ accessories/hanging-chains/premierpotrackchains.cfm |

*interface for running parallel chains on multi-core computers*and thereby greatly reducing the time needed to achieve a desired MCMC sample size from JAGS. Give it a look!

## Thursday, September 12, 2013

### Visualization of Metropolis with rejected proposals

In response to a previous post, Tinu Schneider has embellished the Metropolis visualization of Maxwell Joseph so it includes explicit plotting of rejected proposals (shown as grey segments). Here are the animated gif and R code. Many thanks to Tinu for sending this! (The gif, below, may take a few moments to load.)

Created by Pretty R at inside-R.org

# # Animating the Metropolis Algorithm # ------------------------------------------- # source: # http://mbjoseph.github.io/blog/2013/09/08/metropolis/ # Prepare ---------------------------------------- require(sm) # needs 'rpanel' dir.create("metropolis_ex") # Setup # -------------------------------------- # population level parameters mu <- 6 sigma <- 3 # data-values n <- 50 # metropolis iter <- 10000 chains <- 3 sigmaJump <- 0.2 # plotting seq1 <- seq( 1, 200, by= 5) seq2 <- seq(200, 450, by= 10) seq3 <- seq(450, 750, by= 50) seq4 <- seq(750, iter, by=150) xlims <- c(4, 8) ylims <- c(1, 5) zlims <- c(0, 0.7) widthPNG <- 1200 heightPNG <- 750 # Helpers # -------------------------------------- # collect some data (e.g. a sample of heights) gendata <- function(n, mu, sigma) { rnorm(n, mu, sigma) } # log-likelihood function ll <- function(x, muhat, sigmahat){ sum(dnorm(x, muhat, sigmahat, log=TRUE)) } # prior densities pmu <- function(mu){ dnorm(mu, 0, 100, log=TRUE) } psigma <- function(sigma){ dunif(sigma, 0, 10, log=TRUE) } # posterior density function (log scale) post <- function(x, mu, sigma){ ll(x, mu, sigma) + pmu(mu) + psigma(sigma) } geninits <- function(){ list(mu = runif(1, 4, 10), sigma = runif(1, 2, 6)) } jump <- function(x, dist = 0.1){ # must be symmetric x + rnorm(1, 0, dist) } # compute # ---------------------------------------- # init proposed <- array(dim = c(chains, 2, iter)) posterior <- array(dim = c(chains, 2, iter)) accepted <- array(dim = c(chains, iter-1)) # the data x <- gendata(n, mu, sigma) for (c in 1:chains){ props <- array(dim=c(2, iter)) theta.post <- array(dim=c(2, iter)) inits <- geninits() theta.post[1, 1] <- inits$mu theta.post[2, 1] <- inits$sigma for (t in 2:iter){ # theta_star = proposed next values for parameters jumpMu <- jump(theta.post[1, t-1], sigmaJump) jumpSigma <- jump(theta.post[2, t-1], sigmaJump) theta_star <- c(jumpMu, jumpSigma) pstar <- post(x, mu = theta_star[1], sigma = theta_star[2]) pprev <- post(x, mu = theta.post[1, t-1], sigma = theta.post[2, t-1]) lr <- pstar - pprev r <- exp(lr) # theta_star is accepted if posterior density is higher w/ theta_star # if posterior density is not higher, it is accepted with probability r # else theta does not change from time t-1 to t accept <- rbinom(1, 1, prob = min(r, 1)) accepted[c, t - 1] <- accept if (accept == 1){ theta.post[, t] <- theta_star } else { theta.post[, t] <- theta.post[, t-1] } props[, t] <- theta_star } proposed[c, , ] <- props posterior[c, , ] <- theta.post } # Plot # ---------------------------------------------- oldPath <- getwd() setwd("metropolis_ex") sequence <- c(seq1, seq2, seq3, seq4) cols <- c("blue", "purple", "red") png(file = "metrop%03d.png", width=widthPNG, height=heightPNG) for (i in sequence){ par(font.main=1, cex.main=1.2) layout(matrix(1:2, nrow=1), widths=c(1.0, 1.5)) # prepare empty plot plot(0, type="n", xlim=xlims, ylim=ylims, xlab="mu", ylab="sigma", main="Markov chains") for (j in 1:chains) { # first all the gray lines lines(proposed[j, 1, 1:i], proposed[j, 2, 1:i], col="gray") } for (j in 1:chains) { # now the accepteds lines(posterior[j, 1, 1:i], posterior[j, 2, 1:i], col=cols[j]) } text(x=6, y=1.0, paste("Iteration ", i), cex=1.2) sm.density(x=cbind(c(posterior[, 1, 1:i]), c(posterior[, 2, 1:i])), xlab="mu", ylab="sigma", zlab="", xlim=xlims, ylim=ylims, zlim=zlims, col="white") title(main="Posterior density") } dev.off() # combine .png's into one .gif # system('convert2gif.bat') # had to write a .bat-file # system('convert -delay 20 *.png "metropolis.gif"') # doesent work #file.remove(list.files(pattern=".png")) setwd(oldPath) #

## Tuesday, September 10, 2013

### Cool visualization of Metropolis sampling

A cool visualization of Metropolis sampling has been posted by Maxwell Joseph here. I've pasted in his animated gif below, as a teaser for you to check his post (and I figured it would be okay because his work focuses on "parasitic biodiversity" ;-). Thanks to Joseph Young for pointing me to this.

What would also be cool is an animation of Metropolis that shows the proposals at every step, even if rejected.

From Maxwell Joseph's post. |

## Friday, September 6, 2013

### Bayesian methods article wins best paper award in Organizational Research Methods

The article, "The Time Has Come: Bayesian Methods for Data Analysis in the Organizational Sciences," won the 2012 Best Paper Award for the journal, Organizational Research Methods. The publisher has made the article freely available. You can see the announcement and find a link to the article here.
The accompanying programs can be found here .
If you find that article interesting, you would probably also like the article, "Bayesian estimation supersedes the

*t*test," described here.
Subscribe to:
Posts (Atom)