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 ~ 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 ~ 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 ~ bi) +
ylab("Log density of beta")
ggsave(file="/tmp/beta3.pdf", plot=p)
```