Tuesday, August 25, 2015

Explore the Beta Distribution

Here is the R code:

plotbeta = function(a, b) {
    x = seq(0, 1, by=.002)
    dat = NULL
    for(ai in a) {
        for(bi in b) {
            y = dbeta(x, ai, bi, log=T)
            dat_ab = data.frame(x, y, ai, bi)
            dat = rbind(dat, dat_ab)
        }
    }
    dat
}

a = c(.01, .1, .9)
b = a
dat = plotbeta(a, b)
p = ggplot(dat, aes(x, y)) +
    geom_line() +
    facet_grid(ai &#126 bi) +
    ylab("Log density of beta")
ggsave(file="/tmp/beta1.pdf", plot=p)


a = c(1.10023, 5, 50)
b = a
dat = plotbeta(a, b)
p = ggplot(dat, aes(x, y)) +
    geom_line() +
    facet_grid(ai &#126 bi) +
    ylab("Log density of beta")
ggsave(file="/tmp/beta2.pdf", plot=p)


a = c(1.10023, 5, 50)
b = c(.01, .1, .9)
dat = plotbeta(a, b)
p = ggplot(dat, aes(x, y)) +
    geom_line() +
    facet_grid(ai &#126 bi) +
    ylab("Log density of beta")
ggsave(file="/tmp/beta3.pdf", plot=p)

0 comments: