Bernoulli law

Applies to all situations where the data generation process can only result in two mutually exclusive outcomes (e.g., a coin toss). At each trial, if we assume that \(\Pr(\text{heads}) = \theta\), then \(\Pr(\text{tails}) = 1 - \theta\).

Bernoulli process

If we have a series of independent and identically distributed throws \(\{Y_i\}\) (i.e., each throw has a Bernoulli probability distribution with probability \(\theta\)), the set of throws can be described by a Binomial distribution.

Bernoulli process

Knowing that the trials are independent of each other, the probability of obtaining this sequence is \((1 - \theta) \times (1 - \theta) \times (1 - \theta) \times \theta \times \theta\), that is: \(\theta^{2} (1 - \theta)^{3}\).

Binomial coefficient

The binomial coefficient allows computing the number of possible combinations resulting in \(y\) successes for \(n\) throws in the following way (read “\(y\) among \(n\)” or “number of combinations of \(y\) among \(n\)”)¹:

\[ \left(\begin{array}{l} n \\ y \end{array}\right) = C_{y}^{n} = \frac{n !}{y !(n - y) !} \]

# computing the total number of possible configurations in R
choose(n = 3, k = 1)

[1] 3

Binomial distribution

\[ p(y \given \theta) = \Pr(Y = y \given \theta) = \left(\begin{array}{l} n \\ y \end{array}\right) \theta^{y}(1 - \theta)^{n - y} \] The binomial distribution allows us to calculate the probability of obtaining \(y\) successes in \(n\) trials, for a given \(\theta\). Example of the binomial distribution for an unbiased coin (\(\theta = 0.5\)), indicating the probability of obtaining \(n\) heads in 10 throws (in R: dbinom(x = 0:10, size = 10, prob = 0.5)).

Sampling binary data

library(tidyverse)
set.seed(666) # for reproducibility

sample(x = c(0, 1), size = 500, prob = c(0.4, 0.6), replace = TRUE) %>% # theta = 0.6
        data.frame() %>%
        mutate(x = seq_along(.), y = cummean(.) ) %>%
        ggplot(aes(x = x, y = y) ) +
        geom_line(lwd = 1) +
        geom_hline(yintercept = 0.6, lty = 3) +
        labs(x = "Number of trials", y = "Proportion of heads") +
        ylim(0, 1)

Defining the model (likelihood)

Likelihood function:

• We consider \(y\) to be the number of successes.
• We consider the number of observations \(n\) to be a constant.
• We consider \(\theta\) to be the parameter of our model (i.e., the probability of success).

Likelihood versus probability

A coin with a bias \(\theta\) is tossed again (where \(\theta\) represents the probability of getting heads). This coin is tossed twice and we obtain a Heads and a Tails.

Likelihood versus probability

# Graphical representation of the likelihood function for y = 1 and n = 2

y <- 1 # number of heads
n <- 2 # number of trials

data.frame(theta = seq(from = 0, to = 1, length.out = 1e3) ) %>%
  mutate(likelihood = dbinom(x = y, size = n, prob = theta) ) %>%
  ggplot(aes(x = theta, y = likelihood) ) +
  geom_area(color = "orangered", fill = "orangered", alpha = 0.5) +
  labs(x = expression(paste(theta, " - Pr(head)") ), y = "Likelihood")

Likelihood versus probability

If we calculate the area under the curve of this function, we get:

\[\int_{0}^{1} 2 \theta(1 - \theta) \mathrm{d} \theta = \frac{1}{3}\]

f <- function(theta) {2 * theta * (1 - theta) }
integrate(f = f, lower = 0, upper = 1)

0.3333333 with absolute error < 3.7e-15

Likelihood versus probability

Note that the likelihood of \(\theta\) for a particular data item is equal to the probability of the data for this value of \(\theta\). However, the distribution of these likelihoods (in columns) is not a probability distribution. In a usual Bayesian analysis, the data are considered fixed and the value of \(\theta\) is considered a random variable.

Defining the model (prior)

How can we define a prior for \(\theta\)?

Beta distribution

\[ \begin{align} \color{steelblue}{p(\theta \given a, b)} \ &\color{steelblue}{= \mathrm{Beta}(\theta \given a, b)} \\ & \color{steelblue}{= \theta^{a - 1}(1 - \theta)^{b - 1} / B(a, b)} \\ & \color{steelblue}{\propto \theta^{a - 1}(1 - \theta)^{b - 1}} \end{align} \]

where \(a\) and \(b\) are two parameters such that \(a \geq 0\), \(b \geq 0\), and \(B(a, b)\) is a normalisation constant.

Interpretation of parameters

• The absence of prior knowledge can be expressed by setting \(a = b = 1\) (orange distribution).
• A prior in favour of an absence of bias can be expressed by setting \(a = b > 1\) (green distribution).
• A bias in favour of Heads can be expressed by setting \(a > b\) (blue distribution).
• A bias in favour of Tails can be expressed by setting \(a < b\) (purple distribution).

Interpretation of parameters

The level of certainty increases with the sum \(\kappa = a + b\).

• No idea where the coin comes from: \(a = b = 1\) -> flat prior.
• While waiting for the experiment to begin, the coin was tossed 10 times and we observed 5 “Heads”: \(a = b = 5\) -> weakly informative prior.
• The coin comes from the Bank of France: \(a = b = 50\) -> strongly informative prior.

Interpretation of parameters

Suppose we have an estimate of the most likely value of the parameter \(\theta\). We can re-parametrise the Beta distribution as a function of the mode \(\omega\) and the level of certainty \(\kappa\):

\[ \begin{align} a &= \omega(\kappa - 2) + 1 \\ b &= (1 - \omega)(\kappa - 2) + 1 &&\mbox{for } \kappa > 2 \end{align} \]

Conjugate prior

Formally, if \(\mathcal{F}\) is a class of sampling distributions \(p(y \given \theta)\), and \(\mathcal{P}\) is a class of prior distributions for \(\theta\), then the class \(\mathcal{P}\) is conjugate to \(\mathcal{F}\) if:

\[ p(\theta \given y) \in \mathcal{P} \text{ for all } p(\cdot \given \theta) \in \mathcal{F} \text{ and } p(\cdot) \in \mathcal{P} \]

(p.35, Gelman et al., 2013). In other words, a prior is called a conjugate prior if, when converted to a posterior by being multiplied by the likelihood, it keeps the same form. In our case, the Beta prior is a conjugate prior for the Binomial likelihood, because the posterior is a Beta distribution as well.

Analytical derivation of the posterior distribution

Assume a prior defined by: \(\ \color{steelblue}{p(\theta \given a, b) = \mathrm{Beta}(a, b) = \frac{\theta^{a - 1}(1 - \theta)^{b - 1}}{B(a, b)} \propto \theta^{a - 1}(1 - \theta)^{b - 1}}\)

An example to help you digest

We observe \(y = 7\) correct answers out of \(n = 10\) questions. We choose a prior \(\mathrm{Beta}(1, 1)\), that is, a uniform prior on \([0, 1]\). This prior is equivalent to a prior knowledge of 0 successes and 0 failures (i.e., a flat prior).

An example to help you digest

Influence of prior on posterior distribution

Case where \(n < a + b, (n = 10, a = 4, b = 16)\).

Influence of prior on posterior distribution

Case where \(n = a + b, (n = 20, a = 4, b = 16)\).

Influence of prior on posterior distribution

Case where \(n > a + b, (n = 40, a = 4, b = 16)\).

Take-home message

The posterior distribution is always a compromise between the prior distribution and the likelihood function (Kruschke, 2015).

Take-home message

The more data we have, the less influence the prior has in estimating the posterior distribution (and vice versa). Warning: When the prior assigns a probability of 0 to certain values of \(\theta\), the model is unable to learn (these values are then considered “impossible”).

Marginal likelihood

\[ \text{Posterior} = \frac{\text{Likelihood} \times \text{Prior}}{\text{Marginal likelihood}} \propto \text{Likelihood} \times \text{Prior} \]

\[ p(\theta \given \text{data}) = \frac{p(\text{data} \given \theta) \times \ p(\theta)}{p(\text{data})} \propto p(\text{data} \given \theta) \times p(\theta) \]

Marginal likelihood

New issue: \(p(\text{data})\) is obtained by calculating the sum (for discrete discrete variables) or the integral (for continuous variables) of the joint density \(p(\text{data}, \theta)\) over all possible values of of \(\theta\). This can become tricky when the model includes many continuous parameters…

Marginal likelihood

There are three ways to get around this problem:

• Analytical solution \(\longrightarrow\) Using a conjugate prior (e.g., the Beta-Binomial model).

• Discretised solution \(\longrightarrow\) Computing the posterior on a finite set of points (grid method).

• Approximated solution \(\longrightarrow\) The parameter space is “cleverly” sampled (e.g., MCMC methods, cf. Course n°03).

Discrete distributions

Continuous distributions

Problem: This solution is very restrictive. Ideally, the model (likelihood + prior) model should be defined on the basis of the interpretation of the parameters of these distributions, and not to facilitate the calculations…

Posterior distribution, grid method

Posterior distribution, grid method

Posterior distribution, grid method

Posterior distribution, grid method

Posterior distribution, grid method

Posterior distribution, grid method

Problem with the number of parameters… Refining the grid increases the calculation time:

• 3 parameters with a \(10^3\) grid = \(10^9\) calculation points
• 10 parameters with a grid of \(10^3\) nodes = \(10^{30}\) calculation points

Another solution: Sampling the posterior distribution

To sample (cleverly) a posterior distribution, we can use different implementations of MCMC methods (e.g., Metropolis-Hastings, Hamiltonian Monte Carlo) which we will discuss in Course n°03. In the meantime, we will work with samples from the posterior distribution i) to get used to results from MCMC methods and ii) because it is simpler to compute summary statistics (e.g., mean or credible intervals) on samples rather than by computing integrals.

p_grid <- seq(from = 0, to = 1, length.out = 1000) # creates a grid
prior <- rep(1, 1000) # uniform prior
likelihood <- dbinom(x = 12, size = 20, prob = p_grid) # computes likelihood
posterior <- (likelihood * prior) / sum(likelihood * prior) # computes posterior
samples <- sample(x = p_grid, size = 1e3, prob = posterior, replace = TRUE) # sampling
hist(samples, main = "", xlab = expression(theta) ) # histogram

Posterior distribution, summary

Analytical solution for the Beta-Binomial model:

a <- b <- 1 # parameters of the Beta prior
n <- 9 # number of observations
y <- 6 # number of successes
p_grid <- seq(from = 0, to = 1, length.out = 1000)
posterior <- dbeta(p_grid, y + a, n - y + b) # plot(posterior)

p_grid <- seq(from = 0, to = 1, length.out = 1000)
prior <- rep(1, 1000) # uniform prior
likelihood <- dbinom(x = y, size = n, prob = p_grid)
posterior <- (likelihood * prior) / sum(likelihood * prior) # plot(posterior)

samples <- sample(x = p_grid, size = 1e4, prob = posterior, replace = TRUE) # hist(samples)

Posterior distribution, summary

Analytical solution

• The posterior distribution is described explicitly
• The model is strongly constrained (i.e., we have to pick conjugate priors)

Using samples to summarise the posterior distribution

Estimation of the central tendency: From a set of samples of samples from a posterior distribution, we can calculate the mean mean, mode, and median. For example, for a uniform prior, 10 coin tosses and 3 heads.

mode_posterior <- find_mode(samples) # in blue
mean_posterior <- mean(samples) # in orange
median_posterior <- median(samples) # in green

Using samples to summarise the posterior distribution

What is the probability that the bias of the coin \(\theta\) is greater than 0.5?

sum(samples > 0.5) / length(samples) # equivalent to mean(samples > 0.5)

[1] 0.112

sum(samples > 0.2 & samples < 0.4) / length(samples)

[1] 0.5482

Highest density interval (HDI)

Properties of the highest density interval:

• The HDI indicates the most likely values for the parameter (given the data and the priors)
• The narrower the HDI, the greater the degree of certainty
• The width of the HDI decreases as the number of measurements increases

Highest density interval (HDI)

Highest density interval (HDI)

library(imsb)

set.seed(666)
p_grid <- seq(from = 0, to = 1, length.out = 1e3)
pTheta <- dbeta(p_grid, 3, 10)
massVec <- pTheta / sum(pTheta)
samples <- sample(x = p_grid, size = 1e4, replace = TRUE, prob = pTheta)

posterior_plot(samples = samples, credmass = 0.89)

Region of practical equivalence (ROPE)

This procedure can be used to accept or reject a null value. The region of practical equivalence or region of practical equivalence (ROPE) defines an interval of values that are considered to be “equivalent” to the null value. The figure below summarises the possible decisions resulting from this procedure (Kruschke, 2018).

Region of practical equivalence (ROPE)

The value of the parameter (e.g., \(\theta = 0.5\)) is rejected if the HDI is entirely outside the ROPE. The parameter value (e.g., \(\theta = 0.5\)) is accepted if the HDI is entirely within the ROPE. If the HDI and the ROPE overlap, we cannot conclude…

posterior_plot(samples = samples, rope = c(0.49, 0.51) ) +
    labs(x = expression(theta) )

Model comparison

We tossed coin 200 times and got 115 “Heads”. Is the coin biased? We build two models and try to find out which one best accounts for the data.

Model comparaison

The Bayes factor is the ratio of the marginal likelihoods (of the two models).

\[\frac{p(\mathcal{M}_{0} \given \text{data})}{p(\mathcal{M}_{1} \given \text{data})} = \color{green}{\frac{p(\text{data} \given \mathcal{M}_{0})}{p(\text{data} \given \mathcal{M}_{1})}} \times \frac{p(\mathcal{M}_{0})}{p(\mathcal{M}_{1})}\]

Model checking

Two roles of the likelihood function:

• It is a function of \(\theta\) for calculating the posterior distribution: \(\mathcal{L}(\theta \given y, n)\).
• When \(\theta\) is known/fixed, it is a probability distribution: \(p(y \given \theta, n) \propto \theta^y(1 - \theta)^{(n - y)}\).

samples <- rbinom(n = 1e4, size = 10, prob = 0.6)

Model checking

In a Bayesian models, there are two sources of uncertainty when generating predictions:

• Uncertainty related to the sampling process
  -> We draw data from a Binomial pdf
• Uncertainty about the value of \(\theta\)
  -> Our knowledge of \(\theta\) is described by a (posterior) pdf

posterior <- rbeta(n = 1e4, shape1 = 16, shape2 = 10)
samples <- rbinom(n = 1e4, size = 10, prob = posterior)

Prior and posterior predictive checking

Posterior predictive checking

Practical work

An analyst who works in a factory that makes famous Swedish bread rolls read a book that raised a thorny question… Why does the toast always land on the butter side? Failing to come up with a plausible answer, she set out to verify this assertion. The first experiment was to drop a slice of buttered bread from the height of a table. The results of this experiment are available in the tartine1 dataset from the imsb package.

Retrieving the data

First task: Retrieving the data.

# importing the data
data <- open_data(tartine1)

# summary of the data
str(data)

'data.frame':   500 obs. of  2 variables:
 $ trial: int  1 2 3 4 5 6 7 8 9 10 ...
 $ side : int  1 1 0 1 0 0 1 1 1 0 ...

Questions

Proposed solution - Question 1

Since the toast only has two sides, the result is similar to a draw according to a binomial distribution with an unknown parameter \(\theta\). What is the posterior distribution of the parameter \(\theta\) given these data?

# number of trials
nbTrial <- length(data$trial)

# number of "successes" (i.e., when the toast lands on the butter side)
nbSuccess <- sum(data$side)

# size of the grid
grid_size <- 1e3

# generating the grid
p_grid <- seq(from = 0, to = 1, length.out = grid_size)

# uniform prior
prior <- rep(1, grid_size)

# computing the likelihod
likelihood <- dbinom(x = nbSuccess, size = nbTrial, prob = p_grid)

# computing the posterior
posterior <- likelihood * prior / sum(likelihood * prior)

Proposed solution - Question 2

Calculate the 95% HDI of the posterior distribution and give a graphical representation of the result.

samples <- sample(x = p_grid, prob = posterior, size = 1e3, replace = TRUE)
posterior_plot(samples = samples, credmass = 0.95) + labs(x = expression(theta) )

Proposed solution - Question 3

Can the null hypothesis that \(\theta = 0.5\) be rejected? No, because the HDI overlaps with the ROPE…

posterior_plot(
  samples = samples, credmass = 0.95,
  compval = 0.5, rope = c(0.49, 0.51)
  ) + labs(x = expression(theta) )

Proposed solution - Question 4

At this point, no conclusion can be drawn. The analyst decides to repeat a series of observations to refine her results.

data2 <- open_data(tartine2)
str(data2)

'data.frame':   100 obs. of  2 variables:
 $ trial: int  1 2 3 4 5 6 7 8 9 10 ...
 $ side : int  0 0 1 0 0 1 1 1 0 0 ...

nbTrial2 <- length(data2$trial) # number of trials
nbSucces2 <- sum(data2$side) # number of "successes"

Proposed solution - Question 4

We use the previous posterior as the prior for this new model.

mode1 <- find_mode(samples)
prior2 <- dbeta(p_grid, mode1 * (nbTrial - 2) + 1, (1 - mode1) * (nbTrial - 2) + 1)

data.frame(x = p_grid, y = prior2) %>%
  ggplot(aes(x = x, y = y) ) +
  geom_area(alpha = 0.8, fill = "steelblue") +
  geom_line(size = 0.8) +
  labs(x = expression(theta), y = "Probability density")

Proposed solution - Question 4

likelihood2 <- dbinom(x = nbSucces2, size = nbTrial2, prob = p_grid)
posterior2 <- likelihood2 * prior2 / sum(likelihood2 * prior2)
samples2 <- sample(p_grid, prob = posterior2, size = 1e4, replace = TRUE)

posterior_plot(
  samples = samples2, credmass = 0.95,
  compval = 0.5, rope = c(0.49, 0.51)
  ) + labs(x = expression(theta) )