Practical Tutorial: Posteriors, Likelihoods, and Priors

Exploring how interactions between priors and likelihoods influence the posterior.

To get experience designing priors, we will explore how informative priors interact with the likelihood to influence the posterior.

Note

In this activity, we will revisit our simulated data and fitted model from the previous activity: Model of the Mean. All of the necessary code will be included here, but consider completing that activity first if you haven’t already.

Introduction

Bayesian posteriors are influenced by the data (i.e. the likelihood) as well as the priors. If the priors are uninformative, then the data signal prevails and the inferences will be similar (if not equivalent) to a frequentist approach using maximum likelihood. If the priors are informative, then they compete with the data to influence the posterior distribution. Informative priors are particularly influential when the data are weak (e.g. small sample size, high noise-to-signal ratio).

So, how do we design priors to avoid making our analysis too subjective? In most cases, we try to design priors that define a reasonable range of values for each parameter while remaining relatively flat withing the primary parameter space so that the data signal prevails. If we are working in a system where the signal from data is weak and we have objective prior information that we can formalise as informative priors, then this can objectively improve our inferences.

It takes some experience to get a feel for the sensitivity of posteriors to informative priors in different data scenarios. This excercise allows you to play with various scenarios using simulated data with a very simple model.

Model of the Mean

We’ll start with the model of the mean that we used previously. This model specifies a Gaussian likelihood (i.e. normal distribution) for the data with minimally informative priors for the mean \(\mu\) and standard deviation \(\sigma\).

\[ \begin{aligned} y_i &\sim normal(\mu, \sigma) \\ \\ \mu &\sim normal(0, 100) \\ \sigma &\sim uniform(0, 1000) \end{aligned} \]

We will re-run the model first without any changes and then we’ll come back to modify the priors.

Stan Model

data {
  int<lower=0> n; // sample size
  vector[n] y; // data
}
parameters {
  real mu; // mean
  real<lower=0> sigma; // standard deviation
}
model {
  // likelihood
  y ~ normal(mu, sigma);
  
  // priors
  mu ~ normal(0, 100); 
  sigma ~ uniform(0, 1000);
}

Important

Save this code into a file called model_of_the_mean.stan in your R working directory for this exercise.

Analysis Script

Once you have model_of_the_mean.stan saved into your R working directory, run the following script to simulate new data and refit the model using MCMC.

First, let’s just setup our R environment and working directory.

# cleanup
rm(list = ls())
gc()
cat("\014")
try(dev.off(), silent = T)

# check working directory
getwd()

# [optional] update/fix working directory
# setwd(file.path('/my/working/directory/model_of_the_mean'))

# directories
wd <- getwd()

Then, we’ll simulate data and run the model.

# load libraries
library(cmdstanr)

#---- SIMULATE DATA ----#

# set seed for random number generators
seed <- round(runif(1, 1, 1e6))
set.seed(seed)

# simulated population parameters
mu <- 5 # true population mean
sigma <- 3 # true population standard deviation

# simulate observations
n <- 1e3 # sample size
y <- rnorm(n = n, mean = mu, sd = sigma)

# Stan model data
md <- list(
  n = n,
  y = y,
  mu_true = mu,
  sigma_true = sigma,
  seed = seed
)

# save data to disk
saveRDS(object = md, file = file.path(wd, "md.rds"))


#---- MCMC ----#

# compile the stan model
mod <- cmdstan_model(file.path(wd, "model_of_the_mean.stan"))

# function to generate initial values
init_generator <- function(md = md, chain_id = 1) {
  result <- list()

  result[["mu"]] <- runif(1, -1e3, 1e3)
  result[["sigma"]] <- runif(1, 0, 1e3)

  return(result)
}

# MCMC configuration
chains <- 4
warmup <- 1e3
samples <- 2e3
inits <- lapply(1:chains, function(id) init_generator(md = md, chain_id = id))

# run MCMC
fit <- mod$sample(
  data = md,
  parallel_chains = chains,
  init = inits,
  iter_sampling = samples,
  iter_warmup = warmup,
  save_warmup = TRUE,
  seed = md$seed
)

# save fitted model to disk
fit$save_object(file = file.path(wd, "fit.rds"))

Visualise Posterior, Likelihood, and Priors

Now, let’s visualise the distributions of the priors, likelihood, and posterior for our base model.

A few helper functions are provided to help with this. Understanding the helper functions is not essential for this exercise. They are provided for your convenience and will be used in subsequent code chunks.

# log-likelihood function
# note: the density function--e.g. dnorm()-- should match the distribution of the likelihood from the model being assessed

Lfun <- function(mu, sigma, dat = md$y) {
  sum(dnorm(x = dat, mean = mu, sd = sigma, log = TRUE))
}

# rescale log-likelihoods to min and max for plotting
rescale_log_likelihoods <- function(log_likelihoods, new_min, new_max) {
  # log-sum-exp trick (thanks, ChatGPT 4o)
  max_log_likelihoods <- max(log_likelihoods)
  likelihoods <- exp(log_likelihoods - max_log_likelihoods)

  # rescale (thanks again, ChatGPT 4o)
  old_min <- min(likelihoods, na.rm = TRUE)
  old_max <- max(likelihoods, na.rm = TRUE)
  scaled_likelihoods <- ((likelihoods - old_min) / (old_max - old_min)) *
    (new_max - new_min) +
    new_min

  return(scaled_likelihoods)
}

# plot likelihood, prior, and posterior
myplot <- function(
  fit,
  md,
  par = "mu",
  xlim = c(4, 6),
  prior_mean = 0,
  prior_sd = 100
) {
  # fit = fit
  # md = md
  # par = "mu"
  # xlim = c(4, 6)
  # prior_mean = 0
  # prior_sd = 100

  # posterior estimate for mu
  # "from" and "to" define the range of values you want to plot
  density_mu <- density(
    fit$draws(par),
    from = xlim[1],
    to = xlim[2]
  )

  # take a look at the result from the density function
  # y = probability densities from our posterior estimate for mu
  # x = values of mu
  density_mu
  # plot(density_mu)

  # extract the probability density values (i.e. the values plotted on the y-axis)
  posterior_mu <- density_mu$y

  # probability density of the prior defined in the model
  #** Note: This should be updated to match the prior for mu that was defined in the model being assessed *
  prior_mu <- dnorm(
    x = density_mu$x,
    mean = prior_mean,
    sd = prior_sd
  )
  # plot(x = density_mu$x, y = prior_mu, type = 'l')

  # likelihood estimate: i.e. the likelihood of the data given mu, i.e. L(x|mu)
  log_likelihood_mu <- sapply(
    density_mu$x,
    function(m) {
      Lfun(
        mu = m,
        sigma = md$sigma_true
      )
    }
  )

  # rescale the log-likelihood values to match scale of posterior probabilities
  # Note: This is necessary just to make a nice plot with posteriors, priors, and likelihoods together

  likelihood_mu <- rescale_log_likelihoods(
    log_likelihood_mu,
    new_min = min(posterior_mu),
    new_max = max(posterior_mu)
  )

  ## make plot

  # posterior
  plot(
    x = density_mu$x,
    y = posterior_mu,
    type = "l",
    lwd = 0,
    col = "grey",
    xlab = "Parameter value",
    ylab = "Probability density or scaled likelihood"
  )

  polygon(
    x = c(density_mu$x, rev(density_mu$x)),
    y = c(posterior_mu, rep(0, length(posterior_mu))),
    col = "grey",
    border = NA
  )

  # prior
  lines(x = density_mu$x, y = prior_mu, lty = 3, lwd = 2)

  # likelihood
  lines(x = density_mu$x, y = likelihood_mu, lty = 2, lwd = 2)

  # true value
  abline(v = md$mu_true, col = "red")

  # legend
  legend(
    "topright",
    legend = c("Posterior", "Prior", "Likelihood", "True value"),
    fill = c("grey", NA, NA, NA),
    border = c(NA, NA, NA, NA),
    lty = c(NA, 3, 2, 1),
    lwd = c(NA, 2, 2, 2),
    col = c(NA, "black", "black", "red")
  )
}

Uninformative Prior

Now, we are ready to plot the likelihood and prior alongside the posterior distribution. We will focus on the mean parameter \(\mu\) from our model.

myplot(
  fit = fit,
  md = md,
  par = "mu",
  xlim = c(3, 7), # lower and upper x-axis values for the plot
  prior_mean = 0, # this should match the mean of the prior on mu in model_of_the_mean.stan
  prior_sd = 1000 # match the sd of the prior on mu in model_of_the_mean.stan
)

Remember that we set an uninformative prior on the \(\mu\) parameter in our model (i.e. \(\mu \sim normal(0, 100)\)), so it appears as a flat line at the bottom of this plot. Notice that our Bayesian posterior is essentially equivalent to the likelihood in this case because the prior is uninformative.

Note

Likelihoods are not probabilities (i.e. they do not range from zero to one and they do not sum to one across all possible outcomes). We have rescaled our likelihood to have the same height as our posterior just so that we can plot them together.

Informative Prior Aligned With True Parameter Value

Now, let’s test a strongly informative prior that is well aligned with the data signal. Here we will use the prior \(\mu \sim normal(5, 0.1)\) (remember, 5 is our true population mean that we defined to generate our simulated data). We are using the exact same simulated dataset and initial values from our original model.

We can see the the posterior now has less uncertainty (i.e. less variance) than the likelihood. Why is that?

Is our posterior well aligned with our likelihood? What are some reasonable scenarios when the posterior would differ from the likelihood under these conditions?

Informative Prior Not Aligned With True Parameter Value

Now, let’s test a strongly informative prior that is not well aligned with the true parameter value. In other words, our prior is strong and wrong (a bad combo)! Here we will use the prior \(\mu \sim normal(4, 0.1)\) (remember, 5 is our true population mean that we defined to generate our simulated data). We are using the exact same simulated dataset and initial values from our original model.

Our posterior now falls somewhere between the prior and the likelihood. Why is that?

What do you think would happen if the prior was slightly less informative (e.g. slightly more variance in the prior)? Can you test that?

Going Further

Now we have a process to explore how priors interact with likelihoods to influence the posterior distribution in various scenarios. Take some time to investigate on your own by modifying the priors and the data to see how results may change. A few questions are provided below to guide you into some interesting topics to explore, but feel free to follow your own intuition and answer your own questions.

The full analysis script is provided below for your convenience.

What effect does a very strong prior have when it is centered very far from the true parameter value? For example, try the prior \(\mu \sim normal(1000, 0.01)\). Note that this prior should be flat within the parameter space occupied by the data (e.g. around 5).

How does a moderately informative prior influence the posterior when the data signal is very weak? Try re-running the data simulation with a very small sample size (e.g. n=3) with a prior like \(\mu \sim normal(4, 1)\).

What would you expect from an informative uniform prior? Try \(\mu \sim uniform(4, 5)\).

Are there other scenarios that think would be interesting to explore? Give it a try!

Full Analysis Script

# cleanup
rm(list = ls())
gc()
cat("\014")
try(dev.off(), silent = T)

# check working directory
getwd()

# [optional] update/fix working directory
# setwd(file.path('/my/working/directory/model_of_the_mean'))

# directories
wd <- getwd()

# load libraries
library(cmdstanr)
library(posterior)
library(bayesplot)

#---- SIMULATE DATA ----#

# set seed for random number generators
seed <- round(runif(1, 1, 1e6))
set.seed(seed)

# simulated population parameters
mu <- 5 # true population mean
sigma <- 3 # true population standard deviation

# SAMPLE SIZE
n <- 1e3

# simulate observations
y <- rnorm(n = n, mean = mu, sd = sigma)

# Stan model data
md <- list(
  n = n,
  y = y,
  mu_true = mu,
  sigma_true = sigma,
  seed = seed
)

#---- MCMC ----#

# compile the stan model
mod <- cmdstan_model(file.path(wd, "model_of_the_mean.stan"))

# function to generate initial values
init_generator <- function(md = md, chain_id = 1) {
  result <- list()

  result[["mu"]] <- runif(1, -1e3, 1e3)
  result[["sigma"]] <- runif(1, 0, 1e3)

  return(result)
}

# MCMC configuration
chains <- 4
warmup <- 1e3
samples <- 2e3
inits <- lapply(1:chains, function(id) init_generator(md = md, chain_id = id))

# run MCMC
fit <- mod$sample(
  data = md,
  parallel_chains = chains,
  init = inits,
  iter_sampling = samples,
  iter_warmup = warmup,
  save_warmup = TRUE,
  seed = md$seed
)

#---- HELPER FUNCTIONS FOR PLOTTING ----#

# log-likelihood function
Lfun <- function(mu, sigma, dat = md$y) {
  sum(dnorm(x = dat, mean = mu, sd = sigma, log = TRUE))
}

# rescale log-likelihoods to min and max for plotting
rescale_log_likelihoods <- function(log_likelihoods, new_min, new_max) {
  # log-sum-exp trick (thanks, ChatGPT 4o)
  max_log_likelihoods <- max(log_likelihoods)
  likelihoods <- exp(log_likelihoods - max_log_likelihoods)

  # rescale (thanks again, ChatGPT 4o)
  old_min <- min(likelihoods, na.rm = TRUE)
  old_max <- max(likelihoods, na.rm = TRUE)
  scaled_likelihoods <- ((likelihoods - old_min) / (old_max - old_min)) *
    (new_max - new_min) +
    new_min

  return(scaled_likelihoods)
}

# plot likelihood, prior, and posterior
myplot <- function(
  fit,
  md,
  par = "mu",
  xlim = c(4, 6),
  prior_mean = 0,
  prior_sd = 100
) {
  # posterior estimate for mu
  density_mu <- density(
    fit$draws(par),
    from = xlim[1],
    to = xlim[2]
  )

  # extract the probability density values (i.e. the values plotted on the y-axis)
  posterior_mu <- density_mu$y

  # probability density of the prior defined in the model
  prior_mu <- dnorm(
    x = density_mu$x,
    mean = prior_mean,
    sd = prior_sd
  )

  # likelihood estimate: i.e. the likelihood of the data given mu, i.e. L(x|mu)
  log_likelihood_mu <- sapply(
    density_mu$x,
    function(m) {
      Lfun(
        mu = m,
        sigma = md$sigma_true
      )
    }
  )

  likelihood_mu <- rescale_log_likelihoods(
    log_likelihood_mu,
    new_min = min(posterior_mu),
    new_max = max(posterior_mu)
  )

  # plot
  plot(
    x = density_mu$x,
    y = posterior_mu,
    type = "l",
    lwd = 0,
    col = "grey",
    xlab = "Parameter value",
    ylab = "Probability density or scaled likelihood"
  )
  polygon(
    x = c(density_mu$x, rev(density_mu$x)),
    y = c(posterior_mu, rep(0, length(posterior_mu))),
    col = "grey",
    border = NA
  )
  lines(x = density_mu$x, y = prior_mu, lty = 3, lwd = 2)
  lines(x = density_mu$x, y = likelihood_mu, lty = 2, lwd = 2)
  abline(v = md$mu_true, col = "red")
  legend(
    "topright",
    legend = c("Posterior", "Prior", "Likelihood", "True value"),
    fill = c("grey", NA, NA, NA),
    border = c(NA, NA, NA, NA),
    lty = c(NA, 3, 2, 1),
    lwd = c(NA, 2, 2, 2),
    col = c(NA, "black", "black", "red")
  )
}

#---- MAKE THE PLOT ----#

### USER OPTIONS ###

# lower and upper x-axis values for the plot
xlim <- c(3, 7)

# mean of the prior on mu
# note: this should match the prior from your current model (e.g. model_of_the_mean.stan)
prior_mean <- 0

# sd of the prior on mu
# note: this should match the prior from your current model (e.g. model_of_the_mean.stan)
prior_sd <- 1000

#####################

# make the plot
myplot(
  fit = fit,
  md = md,
  par = "mu",
  xlim = xlim,
  prior_mean = prior_mean,
  prior_sd = prior_sd
)