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!)

JAGS users: runjags has been updated

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

The package runjags by Matt Denwood has been recently updated. The runjags package has many nice facilities for running JAGS, especially a seamless 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!

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.)

 
 
#
# 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)
 
#

Created by Pretty R at inside-R.org

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.

From Maxwell Joseph's post.

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

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.