HPD Plotting for Beta-Binomial model

gistfile1.r

require(plyr)
randHDR <- function(df = function(x) {dbeta(x, 1, 1)}, 
	   	    getRand = function(n) {rbeta(n, 1, 1)}, 
		    p = 0.05,
                    n = 1000) {
    # Density quantile approach of Hyndman (1996)
    # Compute HDR level
    xs <- getRand(n)
    ps <- sort(aaply(xs, 1, df))
    fj <- ps[ceiling(p*n)]

    # Find endpoints when crossing level fj
    inhdr = FALSE
    ep <- c()
    for (x in sort(xs)) {
        if (df(x) > fj) {
            if (!inhdr) {
	        ep <- c(ep, x)
	        inhdr <- TRUE
    	    }
        } else {
            if (inhdr) {
    	        ep <- c(ep, x)
	        inhdr <- FALSE
    	    }
        }
    }
    if (inhdr) {
        ep <- c(ep, x)
    }
    
    # If there are no crosses, this must be approximately uniform
    if (length(ep) != 0) {
        return(ep)
    } else {
	return(c(min(xs), max(xs)))
    }
}

plotPriorBeta <- function(flips, a0, b0) {
    a <- a0
    b <- b0
    N = length(flips)
    plot(c(0, 1), c(1, N), col = 0)
    i = 1
    n = 80
    i <- i + 1
    for (flip in flips) {
        if (!flip) {
	    a <- a+1
	}
	b <- b+1
        ints = randHDR(function(x) dbeta(x, a, b), function(n) rbeta(n, a, b), n = n)
        a_ply(t(array(ints, c(2, length(ints)/2))), 1, 
             function(xs) {lines(x = xs, y = c(i, i), col = "red", lwd = 2)})
        points(x = c(a/(a+b)), y = i, col = "black", cex = 0.3, lwd = 3)
    i <- i + 1
    }
    ints = randHDR(function(x) dbeta(x, a, b), function(n) rbeta(n, a, b), n = n)
    a_ply(t(array(ints, c(2, length(ints)/2))), 1, 
          function(xs) {lines(x = xs, y = c(i, i), col = "black", lwd = 3)})
    points(x = c(a/(a+b)), y = i, col = "green", cex = 0.3, lwd = 3)
    ints = randHDR(function(x) dbeta(x, a0, b0), function(n) rbeta(n, a0, b0), n = n)
    a_ply(t(array(ints, c(2, length(ints)/2))), 1, 
          function(xs) {lines(x = xs, y = c(1, 1), col = "black", lwd = 3)})
    points(x = c(a0/(a0+b0)), y = 1, col = "green", cex = 0.3, lwd = 3)
}

# Sample some flips
p = 0.65
xs = runif(200) < p

plotPriorBeta(xs, 1, 1)