Binomial Likelihood

Let’s tackle the first piece of Bayes’ theorem: the likelihood, \(p(y|\theta)\). Here, we are essentially asking the question: “how likely is the data, given a particular value for \(\theta\)?” Here, we introduce the binomial likelihood function:

\[p(y|\theta)=\theta^y(1-\theta)^{n-y}\]

where \(y\) is the number of successes and \(n\) is the number of trials.

Let’s return to our gold merchant and see how we can express the likelihood in terms of the data the merchant observes. Imagine our merchant collects data by digging each of the 10 x 10 squares in the sample. The merchant finds that 14 out of 100 spaces in the sample contain gold:

Here’s how we express the data in terms of the binomial likelihood function:

\[p(y|\theta)=\theta^{14}(1-\theta)^{100-14}\]

If we compute \(p(y|\theta)\) across the grid of \(\theta\) values from 0 to 1, we get the following probability density:

Reproducing this in R is fairly simple - we could substitute the values into the binomial formula, or just use the built-in dbinom function. We then create a dataframe containing the likelihood for each theta and use ggplot2 from the tidyverse to draw the plot:

library(tidyverse)

# use seq(.01, 1, .01) to generate a sequence from .01 to 1 in .01 increments
likelihood <- dbinom(x = 14, size = 100, prob = seq(.01, 1, .01))

# creating a data frame with 'likelihood' for plotting
df <- data.frame(theta = seq(.01, 1, .01),
                 likelihood)

# create the plot
ggplot(df, aes(x = theta, y = likelihood)) +
  geom_line(size = 1) +
  labs(x = TeX("$\\theta$"),
       y = "Likelihood",
       title = TeX("$p(y|\\theta) = \\theta^{14}(1 - \\theta)^{100 - 14}$"),
       subtitle = "Dashed red line indicates maximum likelihood estimate.",
       caption = "Using nonnormalized likelihood.") +
  geom_vline(xintercept = .14, size = 1.25, color = "red", linetype = 2) +
  scale_x_continuous(breaks = seq(0, 1, by = .10)) +
  theme_minimal(base_size = 14) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.title = element_text(hjust = .5, size = 20))

Not surprisingly, the most likely value of \(\theta\) (the maximum likelihood estimate) is .14. In Bayesian statistics, however, we aren’t content with just the point estimate and instead we use the entire density to express uncertainty about \(\theta\).

Incorporating Prior Information

Imagine that our merchant hears rumours that each 100 square meters of land has a 20% probability of containing gold. We now want to incorporate this prior information into our model. \(p(\theta)\) is the prior probability of \(\theta\) and it can represent past data or even subjective knowledge about the likely value and uncertainty of \(\theta\).

A common approach when choosing priors is to identify a conjugate prior: a formula for expressing the prior that has a similar data structure to that of the likelihood. In this case, a conjugate prior for the binomial likelihood is the beta function:

\[\text{beta}(a,b)=\theta^{a-1}(1-\theta)^{b-1}\]

where \(a\) can represent successes and \(b\) can represent failures.

Notice the similarity between the formulas for the binomial and beta functions. They have identical data structures, which makes the beta a conjugate prior for the binomial likelihood. Let’s use a \(\text{beta}(2,8)\) as a prior, representing our knowledge of a rumored 20% probability of gold.

Similarly, you can reproduce this using the dbeta function in R:

# here we use x as a sequence from .01-1
b_likelihood <- dbeta(x = seq(.01, 1, .01), shape1 = 2, shape2 = 8)

df <- data.frame(theta = seq(.01, 1, .01),
                 b_likelihood)

df %>%
  ggplot(aes(x = theta, y = b_likelihood)) +
  geom_line(size = 1) +
  labs(x = TeX("$\\theta$"),
       y = NULL,
       title = TeX("beta(2,8)")) +
  scale_x_continuous(breaks = seq(0, 1, by = .10)) +
  theme_minimal(base_size = 14) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.title = element_text(hjust = .5, size = 20))

Why use a beta prior instead of another binomial density? The reason is that the binomial density is dependent on a sample size parameter n, whereas the beta density is not. It is more convenient to express our prior knowledge without n.

Getting the Posterior

The posterior probability, \(p(\theta|y)\) takes into account the data and the prior. It is often viewed as a compromise between the data likelihood and the prior probability. This is more easily understood graphically:

Here’s how to produce this figure:

theta <- seq(.01, 1, .01)
prior <- dbeta(x = theta, shape1 = 2, shape2 = 8)

#multiply the binomial likelihood by 100 to express it in 'probability terms' like the beta
likelihood <- dbinom(x = 14, size = 100, prob = theta) * 100
posterior_density <- likelihood * prior

df <- data.frame(theta = theta,
                 likelihood = likelihood,
                 prior = prior,
                 pd = posterior_density)

# use pivot_longer to create a variable out of likelihood vs. prior vs. posterior
data_long <- pivot_longer(df, cols = c(likelihood, prior, pd))

data_long %>%
  ggplot(aes(x = theta, y = value, color = name)) +
  geom_line(size = 1.5) +
  labs(x = TeX("$\\theta$"),
       y = NULL,
       color = NULL,
       title = "Posterior, Likelihood, and Prior") +
  scale_color_discrete(labels = c("Likelihood", "Posterior", "Prior")) +
  scale_x_continuous(breaks = seq(0, 1, by = .10)) +
  theme_minimal(base_size = 16) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        panel.grid.major.y = element_blank(),
        panel.grid.minor.y = element_blank(),
        plot.title = element_text(hjust = .5, size = 20))

The above chart shows the posterior, likelihood, and prior side by side. Notice how the prior (in blue) contains less certainty than the likelihood. The posterior density is simply the prior multiplied by the likelihood, thus it contains information from both sources. It’s a general rule that in instances when there is a lot of data and our prior is mostly uninformative that the likelihood will overwhelm the prior.

If we were to extend this further to make actual conclusions about our uncertainty of the proportion of gold in the hills, we would need to draw random samples from the posterior distribution and generate interval estimates based on these draws. For this post, however, we will stick to a graphical evaluation of \(p(\theta|y)\).