Skip to content

Instantly share code, notes, and snippets.

@tel
Created September 12, 2011 23:47
Show Gist options
  • Save tel/1212793 to your computer and use it in GitHub Desktop.
Save tel/1212793 to your computer and use it in GitHub Desktop.
HPD Plotting for Beta-Binomial model
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment