Beta_distribution.Rmd

---
title: "The beta distribution"
output: html_notebook
---


## Preliminaries

Load necessary libraries.

```{r}
# Load packages
library(tidyverse)
```

We define the `bayes_plot_continuous` function that we can use to plot priors, likelihoods, and posteriors.
 
```{r}
# bayes_plot_continuous takes a prior function 
# and a likelihood function on a range of values, 
# calculates the posterior, and then plots all three functions.
bayes_plot_continuous <- function(prior, likelihood, from, to) {
    # Calculate the area under the likelihood function (by integration).
    likelihood_area <- integrate(likelihood, from, to)
    
    # Scale the likelihood function by its area.
    likelihood_scaled <- function(theta) {
        likelihood(theta)/likelihood_area$value
    }
    
    # Calculate the numerator of the posterior function.
    posterior_numer <- function(theta) {
        prior(theta) * likelihood(theta)
    }
    
    # Calculate the denominator of the posterior function (by integration).
    posterior_denom <- integrate(posterior_numer, from, to)
    
    # The posterior is just the ratio.
    posterior <- function(theta) {
        posterior_numer(theta)/posterior_denom$value
        }
    
    # Plot the posterior function.
    posterior_plot <- ggplot(NULL, aes(x = x, color = col, 
                                       linetype = col)) +
        stat_function(data = tibble(x = c(from, to), 
                                        col = factor(1)),
                      fun = prior) +
        stat_function(data = tibble(x = c(from, to), 
                                        col = factor(2)),
                      fun = posterior) +
        stat_function(data = tibble(x = c(from, to), 
                                        col = factor(3)),
                      fun = likelihood_scaled) +
        theme_bw() +
        theme(panel.grid = element_blank()) +
        labs(title = "Prior, Posterior, and Scaled Likelihood",
             x = expression(theta),
             y = NULL) +
        scale_colour_manual(name = "Function",
                            values = c("blue", "black", "red"),
                            labels = c("Prior",
                                       "Posterior",
                                       "Scaled Likelihood")) +
        scale_linetype_manual(name = "Function",
                              values = c("dotted", "solid", "dashed"),
                              labels = c("Prior",
                                         "Posterior",
                                         "Scaled Likelihood"))
    posterior_plot
}
```


## The beta distribution

It is useful to have probability density functions that are defined on the unit interval (i.e. from 0 to 1). These can serve, for example, as models for the possible values of a proportion. Although in elementary statistics classes we often use normal models for this purpose, if you think about it, normal models really don't make a lot of sense. For one, normal models are functions with domain from negative infinity to infinity, and therefore, even narrow ones with means between 0 and 1 will place some probability in their tails (even if quite small) on values less than 0 and greater than 1. This obviously makes no sense for proportions. Furthermore, normal models are always symmetric, so they can never model skewed distributions.

The family of models known as beta distributions will fix both these problems.

Beta distributions, denoted $beta(\theta \mid a, b)$, depend on two positive real numbers that we will call $a$ and $b$. These are often called "shape" parameters, and their interpretation will be explained later. The probability density function for the beta distribution is given by the following formula, defined for values of $\theta$ between 0 and 1:

$$
p(\theta) = beta(\theta \mid a, b) = \frac{\theta^{(a - 1)} \left(1 - \theta\right)^{(b - 1)}}{B(a, b)}.
$$

Don't worry about the denominator $B(a, b)$. It's just a constant that ensures that the total area under the curve is 1 so that we really do have a probability density function. (It's similar to the term $\frac{1}{\sigma \sqrt{2 \pi}}$ in the formula for the normal distribution.) All the action is in the numerator. To emphasize that, we will often write
$$
p(\theta) \propto \theta^{(a - 1)} \left(1 - \theta\right)^{(b - 1)},
$$
using the "proportional to" symbol and ignore the normalizing constant.^[It's minorly annoying that the exponents are $a - 1$ and $b - 1$ instead of just $a$ and $b$. There are places where that convention simplifies the math, but we won't need to worry about it here.]

Here are some plots for various values of $a$ and $b$. Pay attention to how the values of $a$ and $b$ affect the shape of the function.

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta priors",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 1),
                  aes(color = "(1, 1)")) +
    stat_function(fun = dbeta, args = list(shape1 = 5, shape2 = 5),
                  aes(color = "(5, 5)")) +
    stat_function(fun = dbeta, args = list(shape1 = 10, shape2 = 10),
                  aes(color = "(10, 10)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(1, 1)" = "red",
                                  "(5, 5)" = "blue",
                                  "(10, 10)" = "green"),
                       breaks =  c("(1, 1)", "(5, 5)", "(10, 10)"))
```

Note that $beta(\theta \mid 1, 1)$ is a uniform distribution. If $a = b$, then the beta distribution is symmetric with less spread as $a$ and $b$ increase.

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta priors",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 2, shape2 = 1),
                  aes(color = "(2, 1)")) +
    stat_function(fun = dbeta, args = list(shape1 = 5, shape2 = 1),
                  aes(color = "(5, 1)")) +
    stat_function(fun = dbeta, args = list(shape1 = 10, shape2 = 1),
                  aes(color = "(10, 1)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(2, 1)" = "red",
                                  "(5, 1)" = "blue",
                                  "(10, 1)" = "green"),
                       breaks =  c("(2, 1)", "(5, 1)", "(10, 1)"))
```

Beta distributions with large $a$ and $b = 1$ are severely left-skewed.

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta priors",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 2),
                  aes(color = "(1, 2)")) +
    stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 5),
                  aes(color = "(1, 5)")) +
    stat_function(fun = dbeta, args = list(shape1 = 1, shape2 = 10),
                  aes(color = "(1, 10)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(1, 2)" = "red",
                                  "(1, 5)" = "blue",
                                  "(1, 10)" = "green"),
                       breaks =  c("(1, 2)", "(1, 5)", "(1, 10)"))
```

These are the same as above, but exchanging the roles of $a$ and $b$ resulting in severely right-skewed distributions.

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta priors",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 2, shape2 = 8),
                  aes(color = "(2, 8)")) +
    stat_function(fun = dbeta, args = list(shape1 = 4, shape2 = 6),
                  aes(color = "(4, 6)")) +
    stat_function(fun = dbeta, args = list(shape1 = 6, shape2 = 4),
                  aes(color = "(6, 4)")) +
    stat_function(fun = dbeta, args = list(shape1 = 8, shape2 = 2),
                  aes(color = "(8, 2)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(2, 8)" = "red",
                                  "(4, 6)" = "blue",
                                  "(6, 4)" = "green",
                                  "(8, 2)" = "black"),
                       breaks =  c("(2, 8)", "(4, 6)", "(6, 4)", "(8, 2)"))
```

When $a + b$ is constant, the spread is about the same, but the skewness changes.

There is no reason the values of $a$ and $b$ must be integers. Here is a particular weird example of a bimodal distribution created from small values of $a$ and $b$.

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta priors",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 0.1, shape2 = 0.1),
                  aes(color = "(0.1, 0.1)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(0.1, 0.1)" = "red"),
                       breaks =  c("(0.1, 0.1)"))
```


## Binomial likelihood

We return to the example of binomial data from the previous set of notes: in 18 trials, we observe 12 successes. The likelihood function is expressed as follows:

$$
p(x = 12 \mid \theta) \propto \theta^{12} (1 - \theta)^{6}.
$$

The code below defines the binomial likelihood function for 12 successes in a sample of size 18.

```{r}
likelihood1 <- function(theta) {
    dbinom(x = 12, size = 18, prob = theta)
    }
```

Look at a plot of this likelihood function:

```{r}
likelihood1_plot <- ggplot(tibble(x = c(0,1)), aes(x = x)) +
    stat_function(fun = likelihood1) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Binomial Likelihood",
         x = expression(theta),
         y = NULL)
likelihood1_plot
```


## Combining a beta prior with a binomial likelihood

Recall that the posterior density function is proportional to the prior times the likelihood:

$$
p(\theta \mid x) \propto p(\theta) p(x \mid \theta).
$$

If we substitute the formula for the beta distribution for the prior and the formula for the binomial distribution for the likelihood, we get the following posterior:

$$
\begin{align*}
p(\theta \mid x)
    &\propto \theta^{(a - 1)}
        \left(1 - \theta\right)^{(b - 1)} \theta^{12} (1 - \theta)^{6} \\
    &= \theta^{\left( (a + 12) - 1\right)}
        \left(1 - \theta\right)^{\left( (b + 6) - 1 \right)}
\end{align*}
$$

Look carefully: the posterior distribution is also a beta distribution except now the shape parameters are $a + 12$ and $b + 6$.

When dealing with data generated from the binomial distribution, we say that the beta distribution is a *conjugate prior*. A conjugate prior is one that when combined with a specific likelihood function results in a posterior that comes from the same family. A beta prior, when combined with a binomial likelihood, results in a beta posterior. Keep in mind that the beta distribution will not be a conjugate prior for all likelihood functions. It only works in the presence of the binomial likelihood.


## Visualization

Using the `bayes_plot_continuous` function, we can see the effect of using various beta priors.

First, we'll use a uniform distribution, which can be represented as a $beta(\theta \mid 1, 1)$

```{r}
prior1 <- function(theta) {dbeta(theta, shape1 = 1, shape2 = 1)}
bayes_plot_continuous(prior1, likelihood1, from = 0, to = 1)
```

So it's no surprise that the posterior is just the scaled likelihood. The bigger surprise is that the posterior is simply $beta(\theta \mid 13, 7)$:

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    stat_function(fun = dbeta, args = list(shape1 = 13, shape2 = 7),
                  aes(color = "(13, 7)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(13, 7)" = "red"),
                       breaks =  c("(13, 7)"))
```

What if we start with a $beta(\theta \mid 5, 5)$ prior?

```{r}
prior2 <- function(theta) {dbeta(theta, shape1 = 5, shape2 = 5)}
bayes_plot_continuous(prior2, likelihood1, from = 0, to = 1)
```

The posterior is just $beta(\theta \mid 17, 11)$:

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    labs(title = "Beta prior",
         x = expression(theta),
         y = NULL) +
    stat_function(fun = dbeta, args = list(shape1 = 17, shape2 = 11),
                  aes(color = "(17, 11)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(17, 11)" = "red"),
                       breaks =  c("(17, 11)"))
```


## Interpretation of shape parameters

The prior distribution is supposed to encode previous information. Then the data comes in (through the likelihood) and changes the prior into the posterior. In the examples above, whatever $a$ and $b$ represented in the prior, the posterior simply added 12 to $a$ and 6 to $b$. Since the data was 12 successes and 6 failures, this suggests that $a$ has something to do with the number of successes and $b$ has something to do with the number of failures.

Consider again $beta(\theta \mid 1, 1)$. This is a uniform distribution which is supposed to reflect complete ignorance. How many prior successes and failures would we have observed if we had complete ignorance? None!

On the other hand, look above at the graphs of $beta(\theta \mid a, b)$ where $a = b$, say, for example, $beta(\theta \mid 5, 5)$. It becomes more peaked around 0.5, or 50%. What prior information about successes and failures would lead us to believe that the true proportion was likely around 50%? An equal number of successes and failures.

And what prior information would make us believe even more strongly that the true proportion was close to 50%? In other words, when would our prior be even more strongly peaked around 50%? A lot of successes and failures, in equal numbers. Look back at the graph of $beta(\theta \mid 10, 10)$.

Finally, what if you observed a disproportionate number of successes and relatively few failures. What would you believe about the true value of $\theta$? It would be closer to 1. Again, look back at the graphs with $a > b$ and see that they are all skewed with peaks closer to 1 and long tails of less probability approaching 0. The skewness reverses directions if $b > a$. Think about why that makes sense.

All of this together suggests that a $beta(\theta \mid a, b)$ can be thought of as the distribution that assumes you've already observed $a - 1$ successes and $b - 1$ failures. (The "minus one" part is necessary if you want to believe that $beta(\theta \mid 1, 1)$ represents no prior information.) And larger values of $a$ and $b$ correspond to larger sample sizes. The total prior sample size would be $a + b - 2$ (with $a - 1$ successes plus $b - 1$ failures).

Now what happens to $a$ and $b$ when you combine the prior with the data (in the form of the likelihood)? You get a beta with a new value of $a$ that is equal to the old value of $a$ plus the number of successes in the data, and a new value of $b$ that is equal to the old value of $b$ plus the number of failures in the data. In other words, your posterior reflects the successes you "imagined" to be part of the prior combined with an observed number of additional successes. Same goes for failures.


## Setting beta priors

Keep in mind that when formulating a prior in the real world, you likely do not have access to actual data in the form of prior successes and failures. Your job in formulating a realistic prior is to imagine a sample size that is commensurate with your idea of how "strong" your prior is, and to set the sizes of $a$ and $b$ relative to each other to be proportional to your prior idea of how much more/less you expect successes or failures.

Some people find it easier to set a beta prior based on location and spread. For example, perhaps I believe the prior truth to be 25% but I'm willing to acknowledge a range of values between 5% and 45%. One easy way to do this is to use the mode of the beta distribution.

Why the mode instead of the mean? Since the beta distribution is usually skewed, the mode is a better measure of center. After all, the mode is the peak of the distribution where most of the probability mass is.^[This is true for one-dimensional distributions, but paradoxically is not true in higher dimensions. In fact, using the mode as the "most probable value" is extremely problematic in higher dimensions. See https://mc-stan.org/users/documentation/case-studies/curse-dims.html for more information if you're curious.]

The mode $\omega$ of a beta distribution (which only makes sense if $a, b > 1$) is given by

$$
\omega = \frac{a - 1}{a + b - 2}.
$$

Now this one equation is not enough to solve uniquely for $a$ and $b$. One more equation we can use is the *concentration* $\kappa$:

$$
\kappa = a + b.
$$

Remember that $a + b$ is intuitively like a "prior sample size". (Technically, the equivalent prior sample size was $a + b - 2$, but it's close.) So we can pick $\kappa$ to represent the strength of our prior convictions.

If we use some algebra to solve the two equations above for $a$ and $b$, we get the following:

$$
\begin{align*}
a &= \omega \left( \kappa - 2 \right) + 1 \\
b &= \left( 1 - \omega \right) \left( \kappa - 2 \right) + 1
\end{align*}
$$

In practice, it's easier to set the mode and then experiment with values of $\kappa$ that give us the spread we desire in our prior.

We know we want $\omega = 0.25$. If we try out $\kappa = 10$ just for fun, we get:

```{r}
a1 <- 0.25 * (10 - 2) + 1
a1
```

```{r}
b1 <- (1 - 0.25) * (10 - 2) + 1
b1
```

Let's graph that prior and see if it gives us the spread we want:

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    stat_function(fun = dbeta, args = list(shape1 = 3, shape2 = 7),
                  aes(color = "(3, 7)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(3, 7)" = "red"),
                       breaks =  c("(3, 7)"))
```

That's too wide for what we proposed in our prior. Let's increase the value of $\kappa$ to 20:

```{r}
a2 <- 0.25 * (20 - 2) + 1
a2
```

```{r}
b2 <- (1 - 0.25) * (20 - 2) + 1
b2
```

```{r}
ggplot(tibble(x = c(0, 1)), aes(x)) +
    theme_bw() +
    theme(panel.grid = element_blank()) +
    stat_function(fun = dbeta, args = list(shape1 = 5.5, shape2 = 14.5),
                  aes(color = "(5.5, 14.5)")) +
    scale_color_manual(name = "(a, b)",
                       values = c("(5.5, 14.5)" = "red"),
                       breaks =  c("(5.5, 14.5)"))
```

That's quite a bit better.


## Do we need to use conjugate priors in Bayesian data analysis?

Before the age of fast computers, it was nearly impossible to compute posterior density functions. However, when using conjugate priors, the job became much easier because of the simple and easily understandable form of the posterior. (No calculus required!) Therefore, the use of conjugate priors (when they could be found) was practically the only way to do Bayesian data analysis historically.

Now that we have fast computers, this is much less of an issue. Modern MCMC samplers don't care one whit about conjugate priors; they do the same amount of work regardless of the form of the prior. Therefore, we will view the use of conjugate priors more as an interesting historical practice and not a computationally needful thing.