Bayesian inference with mathematical models

Estimated time: 85 minutes

Slides

Introduction

This session introduces the fundamental concepts of Bayesian inference through simulation with a simple epidemiological model. We will build intuition by manually calculating priors, likelihoods, and posteriors using an SIR (susceptible-infectious-recovered) epidemic model, before turning to algorithms that explore the posterior in the next session.

By the end of the course you will have met several different ways of doing this, and the closing session puts them side by side so you can work out which one your own problem needs. Ask that of each method as you meet it rather than saving it for the end. Each one can do something the others cannot, and gives something up to do it.

Objectives

The aim of this first session is to set you up with a framework for model fitting using simulation-based intuition. We will use Distributions.jl for statistical operations and build understanding through hands-on simulation and manual calculations, before moving on to automated inference.

In this session, you will

1. understand priors, likelihood, and posterior through simulation
2. manually try parameter combinations to see how they shape the posterior
3. explore the posterior systematically through grid search
4. see why we will need algorithmic approaches in later sessions

Setup

Source file

The source file of this session is located at sessions/introduction.qmd.

Load libraries

using Random ## for random numbers
using Distributions ## for probability distributions
using DifferentialEquations ## for differential equations
using DataFrames ## for data frames
using CSV ## for CSV file reading
using Plots ## for plots
using StatsPlots ## for statistical plots
using DrWatson

NoteComing from R?

using X in Julia is equivalent to library(X) in R — it loads a package and makes its exports available.

# R
library(deSolve)
library(ggplot2)

Initialisation

We set a random seed for reproducibility. Setting this ensures that you should get exactly the same results on your computer as we do.

Random.seed!(1234)

TaskLocalRNG()

A deterministic SIR model in Julia

We’ll work with a simple deterministic SIR model. We parameterise it with the basic reproductive number $R_0$ and the infectious period $D_\mathrm{inf}$, from which we derive the transmission rate $\beta = R_0/D_\mathrm{inf}$ and the recovery rate $\nu = 1/D_\mathrm{inf}$ that appear in the equations:

$$\begin{cases} \begin{aligned} \frac{\mathrm{d}S}{\mathrm{d}t} &= - \beta S \frac{I}{N}\\ \frac{\mathrm{d}I}{\mathrm{d}t} &= \beta S \frac{I}{N} - \nu I\\ \frac{\mathrm{d}R}{\mathrm{d}t} &= \nu I\\ \end{aligned} \end{cases}$$

To use this model in Julia we first define a function that represents the ODEs:

# Define ODE function
function sir_ode!(du, u, p, t)
    S, I, R = u
    R_0, D_inf = p
    N = S + I + R

    β = R_0 / D_inf
    ν = 1 / D_inf

    du[1] = -β * S * I / N  # dS/dt
    du[2] = β * S * I / N - ν * I  # dI/dt
    du[3] = ν * I  # dR/dt
end

TipWhy is there an exclamation mark in sir_ode!()?

In Julia, the exclamation mark indicates that the function can modify its arguments. In this case, the function modifies du in the last 3 lines to indicate the current change in the system as given by the ODEs.

NoteComing from R?

In R’s deSolve, ordinary differential equation (ODE) functions return list(c(dS, dI, dR)) and take arguments in the order (t, state, parameters). In Julia’s DifferentialEquations.jl, the function modifies du in place (note the ! convention) and the argument order is (du, u, p, t).

# R (deSolve)
sir_ode <- function(t, state, parameters) {
  with(as.list(c(state, parameters)), {
    N <- S + I + R
    beta <- R_0 / D_inf
    nu <- 1 / D_inf
    dS <- -beta * S * I / N
    dI <- beta * S * I / N - nu * I
    dR <- nu * I
    list(c(dS, dI, dR))
  })
}

In order to simulate this model, i.e. numerically integrate the differential equations, we define it as an ODEproblem (which tells Julia that the function we defined above represents the dynamics of the system) which we then solve:

## define initial state
init_state = [999, 1, 0]
## define start and end times
t0, tend = (0, 30)
## define parameters: R_0 = 2.5, D_inf = 2
sir_params = [2.5, 2]
# Define ODE problem
prob = ODEProblem(sir_ode!, init_state, (t0, tend), sir_params)
 # Solve ODE
sol = solve(prob, saveat = 1)

TipMore details on these commands

Here is a closer look at what each of these commands does:

init_state = [999, 1, 0]

This creates an 1-dimensional Array (i.e., a vector) with three elements that we will later use to represent the initial state of our forward simulation (999 susceptibles, 1 infective and 0 recovered). The output tells us that the elements of the vector are of type Int64, i.e. integer-valued numbers.

t0, tend = (0, 30)

This assigns values to the variables t0 (0) and tend (30). It is just a shorter version of writing

t0 = 0
tend = 30

prob = ODEProblem(sir_ode!, init_state, (t0, tend), sir_params)

This defines our ODE system. The first argument to ODEproblem() is a function that contains the equations of the ODE, in our case the sir_ode!() function we defined earlier. The second argument is the initial state. The third argument defines that start and end time of the simulation, here given as a tuple indicated by round parentheses, (t0, tend). A tuple is similar to an array (defined with brackets as we did above in defining init_state) but more flexible in allowing different data types, as in e.g. lists in R. Here we could also have used a vector [t0, tend]. Either will be correctly interpreted by ODEProblem(). The fourth argument gives the parameters of the model sir_params as a vector.

sol = solve(prob, saveat = 1)

This solves the ODE problem given by prob. The saveat argument tells the solver to save state information at time steps of size 1. If this is not given, the time steps at which information is saved will be chosen by the ODE solver - this is often fine but as we will later want to compare this to observations at integer time steps we here save information at these distinct steps.

NoteComing from R?

In R, deSolve::ode() defines and solves the ODE in one call. Julia splits this into two steps: ODEProblem() defines the problem, then solve() runs it.

# R (deSolve) — one call does everything
sol <- ode(
  y = c(S = 999, I = 1, R = 0),
  times = seq(0, 30, by = 1),
  func = sir_ode,
  parms = c(R_0 = 2.5, D_inf = 2)
)

The sol object now contains information on the system dynamics (contained in sol.u) at different times (contained in sol.t).

We can plot this using

plot(sol, labels=["S" "I" "R"])

Part 1: Understanding priors

What is a prior?

A prior represents what we believe about parameters before seeing any data. It encodes our background knowledge, expertise, or assumptions about plausible parameter values.

For our SIR model, let us set the priors to:

• $R_0 \sim \text{Uniform}(1, 20)$: We believe $R_0$ could be anywhere from 1 to 20, with all values equally likely
• $D_\mathrm{inf} \sim \text{Uniform}(1, 14)$: We believe infectious periods from 1 to 14 days are equally plausible

Warning

Uniform distributions are a common first choice, but the hard boundary at the edges can cause numerical problems for many samplers. We use them here for simplicity; we’ll introduce better-behaved priors below.

We can explore what these priors mean through simulation:

# Define our priors as Distribution objects
R_0_prior = Uniform(1.0, 20.0)
D_inf_prior = Uniform(1.0, 14.0)

# Sample from priors many times
n_samples = 1000
R_0_samples = rand(R_0_prior, n_samples)
D_inf_samples = rand(D_inf_prior, n_samples)

println(
    "R_0 prior samples - mean: ", round(mean(R_0_samples), digits = 2),
    ", range: ", round(minimum(R_0_samples), digits = 2), " to ",
    round(maximum(R_0_samples), digits = 2)
)
println(
    "D_inf prior samples - mean: ", round(mean(D_inf_samples), digits = 2),
    ", range: ", round(minimum(D_inf_samples), digits = 2), " to ",
    round(maximum(D_inf_samples), digits = 2)
)

R_0 prior samples - mean: 10.94, range: 1.0 to 20.0
D_inf prior samples - mean: 7.56, range: 1.0 to 13.98

TipCode explainer

• Uniform(a, b) - creates a uniform distribution object from a to b. This is from Distributions.jl.
• rand(dist, n) - draws n random samples from the distribution dist. You can also use rand(dist) for a single sample.

NoteComing from R?

R uses function families for each distribution: runif/dunif/punif/qunif for random draws, density, cumulative distribution function and quantiles respectively. In Distributions.jl, you create a single distribution object and use generic functions:

Operation	R	Julia
Create	(implicit)	dist = Uniform(1, 20)
Random draw	runif(10, 1, 20)	rand(dist, 10)
Density	dunif(x, 1, 20)	pdf(dist, x)
Log density	dunif(x, 1, 20, log=TRUE)	logpdf(dist, x)
CDF	punif(x, 1, 20)	cdf(dist, x)
Quantile	qunif(p, 1, 20)	quantile(dist, p)

The same pattern works for all distributions: Normal, Poisson, Gamma, etc.

# Visualise our priors
p1 = histogram(
    R_0_samples, title = "R_0 Prior", xlabel = "R_0", ylabel = "Frequency",
    alpha = 0.7, label = "Uniform(1, 20)"
)
p2 = histogram(
    D_inf_samples, title = "D_inf Prior", xlabel = "D_inf (days)",
    ylabel = "Frequency", alpha = 0.7, label = "Uniform(1, 14)"
)
plot(p1, p2, layout = (1,2))

Prior predictive simulation

Priors encode our beliefs about parameters, but we usually also have expectations about the kinds of outcomes the model should produce. A prior predictive check connects the two: we draw parameter values from the priors, simulate the model, and inspect the resulting trajectories. If those trajectories look implausible, that tells us to reconsider either the priors themselves or how the model uses the parameters to produce outcomes.

Let’s do this for our SIR model by simulating epidemics using random draws from our priors:

# Prepare plot where we'll add the individual trajectories
p = plot(
    title = "Prior Predictive Trajectories",
    xlabel = "Time", ylabel = "Infected"
)

ntraj = 20

for i in 1:ntraj
    R_0 = R_0_samples[i]
    D_inf = D_inf_samples[i]

    prob = ODEProblem(
        sir_ode!, init_state, (t0, tend), [R_0, D_inf]
    )
    sol = solve(prob)
    plot!(p, sol, alpha = 0.6, label = "", color = :grey, idxs = (0, 2))
end

plot!(p)

NoteHow did we tell Julia to plot just the number of infected?

The plot() commands now contain an argument idxs, which allows us to choose what to plot. This uses the plot recipe for differential equations. Index 0 corresponds to the time variable, and any indices >0 to the variables of the differential equations. In our case, we specify 0, 2, that is we plot time (on the x axis) against the second variable, I (on the y axis).

Note the use of plot! (with !) which adds to an existing plot, following the same convention as push! - the ! indicates the function modifies something in place rather than creating a new object.

NoteComing from R?

plot!() in Julia adds to an existing plot, similar to lines() or points() in base R, or adding layers with + geom_*() in ggplot2.

# Base R
plot(x, y)       # new plot    → Julia: plot(x, y)
lines(x2, y2)    # add to it   → Julia: plot!(x2, y2)

# ggplot2
ggplot(df, aes(x, y)) + geom_line() + geom_point()

TipExercise: Evaluate the prior predictive trajectories

Imagine these trajectories are for a flu-like respiratory infection in a closed population of 1000 people. Are there any that don’t look like how you would expect such an epidemic to unfold?

• Do any grow too fast or too slow?
• Are the peak sizes plausible for a population of 1000?
• Do any last too long or too short?

TipWhat to look for

The flat priors allow $R_0$ up to 20 and an infectious period up to 14 days, and both extremes produce epidemics nobody has seen.

At the top end, a high $R_0$ infects essentially the whole population within a few days, with a peak approaching 1000 — that would mean almost everyone infectious simultaneously. Real respiratory epidemics peak at a small percentage of the population at any one time.

At the bottom end, $R_0$ near 1 produces no epidemic at all: a handful of cases that fade out. Those trajectories are not wrong, but they are not what you are trying to model.

A long infectious period drags the epidemic over months rather than weeks. For a flu-like infection you would expect $R_0$ around 1.2-2 and an infectious period of two to five days, so the flat prior spends most of its probability mass on epidemics that could not happen.

Choosing different priors

What if we had different beliefs? Here are more informative (narrower) priors:

# More informative priors based on epidemiological knowledge.
# R_0 is centred on 2.5, where most respiratory infections sit, and truncated to
# keep it positive. Note that 3% of this prior still sits below 1, i.e. on
# parameter values that produce no epidemic at all.
R_0_informed = truncated(Normal(2.5, 0.8), lower = 0.0)
# The infectious period has mean 2 days, with 95% of the mass between 0.24 and
# 5.6 days and about two thirds of it between 1 and 4.
D_inf_informed = Gamma(2, 1)

# Sample and visualise
R_0_informed_samples = rand(R_0_informed, 1000)
D_inf_informed_samples = rand(D_inf_informed, 1000)

p1 = histogram(
    R_0_samples, xlabel = "R_0", alpha = 0.7, label = "Uniform(1, 20)"
)
histogram!(p1, R_0_informed_samples, alpha = 0.7, label = "truncated Normal(2.5, 0.8)")

p2 = histogram(
    D_inf_samples, xlabel = "D_inf", alpha = 0.7, label = "Uniform(1, 14)"
)
histogram!(p2, D_inf_informed_samples, alpha = 0.7, label = "Gamma(2, 1)")

plot(p1, p2, layout = (1,2))

NoteComing from R?

Watch out for the Gamma parameterisation: Julia’s Gamma(shape, scale) matches R’s scale argument, but many R users default to the rate parameterisation where rate = 1/scale.

# R — these are equivalent:
dgamma(x, shape = 2, scale = 1)
dgamma(x, shape = 2, rate = 1)   # same here since scale = 1

# Julia equivalent:
# Gamma(2, 1)  — always uses scale, never rate

Be careful when scale ≠ 1: Gamma(2, 3) in Julia corresponds to dgamma(x, 2, scale = 3) or dgamma(x, 2, rate = 1/3) in R.

# Compare epidemic predictions
p = plot(
    title = "Prior Predictive: Flat vs Informed",
    xlabel = "Time", ylabel = "Infected"
)

# Flat priors (gray)
for i in 1:10
    R_0 = R_0_samples[i]
    D_inf = D_inf_samples[i]

    prob = ODEProblem(
        sir_ode!, init_state, (t0, tend), [R_0, D_inf]
    )
    sol = solve(prob)
    plot!(p, sol, alpha = 0.6, label = "", color = :gray, idxs = (0, 2))
end

# Informed priors (blue)
for i in 1:10
    R_0 = R_0_informed_samples[i]
    D_inf = D_inf_informed_samples[i]

    prob = ODEProblem(
        sir_ode!, init_state, (t0, tend), [R_0, D_inf]
    )
    sol = solve(prob)
    plot!(p, sol, alpha = 0.6, label = "", color = :blue, idxs = (0, 2))
end

plot!(p)

The informed priors are not just narrower — they are also better behaved, avoiding the hard boundaries of the uniforms. Truncating the Normal at zero keeps $R_0$ positive, and the Gamma is naturally supported on the positive reals.

Priors determine what kinds of epidemics our model considers plausible before seeing data.

TipExercise: Experiment with priors

Try modifying the priors and re-running the prior predictive simulation above. For example:

1. What happens if you use R_0_prior = Uniform(1.0, 5.0) instead of Uniform(1.0, 20.0)?
2. What if you use LogNormal(log(2.5), 0.3) for $R_0$? (Hint: this centres the prior around 2.5 but on a logarithmic scale)
3. Can you find a prior that produces only “realistic-looking” epidemic curves?

There is no single right answer — the point is to build intuition for how prior choices affect the range of predictions.

TipWhat each change does

Uniform(1.0, 5.0) removes the explosive trajectories and leaves a range that is broadly flu-like. Most of what you see will now look like a plausible epidemic.

LogNormal(log(2.5), 0.3) puts roughly 95% of its mass between about 1.4 and 4.5, centred on 2.5. Almost every trajectory then looks realistic, because the prior has done the work of ruling out the rest.

On the third question: you can always find a prior that produces only realistic curves, by making it narrow enough. The catch is that a prior narrow enough to exclude every implausible epidemic also leaves the data little room to change your mind. If the truth were $R_0 = 6$, a LogNormal(log(2.5), 0.3) prior would fight the data hard. That trade-off between ruling out nonsense and staying open to being wrong is the whole difficulty of prior choice, and it is why we check priors by simulating from them rather than by staring at the parameters.

NotePrinciples for choosing priors

1. Use domain knowledge: What do published studies say about plausible parameter ranges? For infectious diseases, systematic reviews often provide estimates.
2. Avoid impossible values: If a parameter must be positive, use a distribution with positive support (Gamma, LogNormal) rather than truncating a Normal.
3. Check with prior predictive simulation: Do samples from your priors generate plausible epidemic trajectories? If not, your priors may be too vague or misspecified.
4. Weakly informative is often better than flat: Completely flat priors can cause computational problems and may place substantial probability on implausible values.
5. Document your choices: Record why you chose each prior - this helps with reproducibility and peer review.

Evaluating prior densities

So far we’ve been sampling from our priors to understand what they represent. But for Bayesian inference, we also need to evaluate how probable specific parameter values are under our priors.

We define some specific parameter combinations to explore throughout our analysis:

# Parameter combinations we'll explore throughout our analysis
test_params = [
    (R_0 = 1.1, D_inf = 2.0),
    (R_0 = 1.5, D_inf = 2.0),
    (R_0 = 2.5, D_inf = 2.0),
    (R_0 = 4.0, D_inf = 2.0),
    (R_0 = 2.5, D_inf = 1.0),
    (R_0 = 2.5, D_inf = 4.0)
]

# Calculate prior densities for both flat and informed priors
log_priors_flat = []
log_priors_informed = []

println("Flat vs Informed prior log-densities:")
for theta in test_params
    # Flat priors
    flat_prior = logpdf(Uniform(1, 20), theta.R_0) + logpdf(Uniform(1, 14), theta.D_inf)
    push!(log_priors_flat, flat_prior)

    # Informed priors
    informed_prior = logpdf(R_0_informed, theta.R_0) + logpdf(D_inf_informed, theta.D_inf)
    push!(log_priors_informed, informed_prior)

    println(
        "R_0 = $(theta.R_0), D_inf = $(theta.D_inf): ",
        "Flat = $(round(flat_prior, digits = 1)), ",
        "Informed = $(round(informed_prior, digits = 1))"
    )
end

Flat vs Informed prior log-densities:
R_0 = 1.1, D_inf = 2.0: Flat = -5.5, Informed = -3.5
R_0 = 1.5, D_inf = 2.0: Flat = -5.5, Informed = -2.8
R_0 = 2.5, D_inf = 2.0: Flat = -5.5, Informed = -2.0
R_0 = 4.0, D_inf = 2.0: Flat = -5.5, Informed = -3.8
R_0 = 2.5, D_inf = 1.0: Flat = -5.5, Informed = -1.7
R_0 = 2.5, D_inf = 4.0: Flat = -5.5, Informed = -3.3

TipCode explainer

• logpdf(dist, x) - returns the log probability density of x under distribution dist. We use log probabilities because they’re numerically more stable (multiplying tiny probabilities can underflow).
• Normal(μ, σ) - creates a normal distribution with mean μ and standard deviation σ.
• Gamma(α, θ) - creates a gamma distribution with shape α and scale θ. Useful for positive parameters like durations.
• push!(array, value) - appends value to the end of array. The ! indicates it modifies the array in place, which you can think of as shorthand for array = push(array, value).

NoteComing from R?

logpdf(dist, x) in Julia is equivalent to calling a density function with log = TRUE in R. For example:

Julia	R
logpdf(Normal(2.5, 0.8), x)	dnorm(x, 2.5, 0.8, log = TRUE)
logpdf(Gamma(2, 1), x)	dgamma(x, 2, scale = 1, log = TRUE)
logpdf(Poisson(10), x)	dpois(x, 10, log = TRUE)

The key difference is that in Julia you first create the distribution object, then pass it to logpdf. In R, you pass the parameters directly to the density function each time.

Part 2: Understanding the likelihood

What is a likelihood?

The likelihood quantifies how probable our observed data is, given specific parameter values. It asks: “If the parameters were X, how likely would we be to see the data we actually observed?”

The prior predictive trajectories we looked at earlier were purely a function of the model and its specification. The likelihood, on the other hand, is about comparing the model to data. We start by loading some observed epidemic data (which we simulated from the model for you, but in principle this could be any data) representing the number of infectious individuals at different time points:

# Load epidemic data from CSV
epi1 = CSV.read(datadir("epi1_synthetic.csv"), DataFrame)

first(epi1, 5)

5×2 DataFrame

Row	time	obs
	Int64	Int64
1	1	1
2	2	3
3	3	4
4	4	11
5	5	19

TipWhat is datadir()?

datadir() comes from the DrWatson.jl package we loaded in the setup. It returns the path to the data/ folder in the project root, so datadir("epi1_synthetic.csv") expands to something like /path/to/mfiidd/data/epi1_synthetic.csv. This keeps file paths portable across different machines.

We first visualise this data set.

# Visualise the observed data
scatter(
    epi1.time, epi1.obs, xlabel = "Time", ylabel = "Observed Cases",
    title = "Epidemic Data", label = "Observed", markersize = 4, color = :red
)

This data was generated from an SIR model with $R_0 ≈ 2.5$ and $D_\mathrm{inf} = 2$ days. If we didn’t know this, could we figure it out by exploring the likelihood?

Exploring likelihood through simulation

Before we compute the likelihood, we need to decide how observations relate to the model’s output. Here we assume each observation is drawn from a Poisson distribution whose mean is the number of infectious individuals predicted by the model at that time point — a common choice for count data. Assuming the observations are conditionally independent given the model, the total log-likelihood is then the sum of log-likelihoods at each time point. We revisit both assumptions in the observation-model section below.

TipExercise: Predict the best fit

Before running the code below, consider the parameter combinations we are about to test. Based on your understanding of SIR dynamics and the data plot:

• Which value of $R_0$ do you expect to give the best fit? Why?
• How would changing $D_\mathrm{inf}$ affect the epidemic curve?

Then run it and check the log-likelihood values. Did the results match your intuition? Is the visual best fit also the one with the highest log-likelihood?

With that in mind, we explore the likelihood by trying the parameter combinations from before:

# Calculate and store likelihood values for our test parameters
log_likelihoods = []

# Generate plot
p = plot(title = "Model Fit Comparison", xlabel = "Time", ylabel = "Cases")

# For each parameter set, calculate log-likelihood and plot
for theta in test_params

    R_0 = theta.R_0
    D_inf = theta.D_inf

    # Simulate from the model
    prob = ODEProblem(
        sir_ode!, init_state, (t0, tend), [R_0, D_inf]
    )
    sol = solve(prob, saveat = epi1.time)

    # I compartment is second compartment. max(..., 0.0) guards against tiny
    # negative values that the ODE solver can return as the epidemic dies out
    I = max.(sol[2, :], 0.0)
    # assuming Poisson log-likelihood at each time point
    obs_dists = [Poisson(It) for It in I]
    # calculate total loglikelihood as sum of log likelihoods at each time point
    total_loglik = sum(logpdf.(obs_dists, epi1.obs))
    # store in vector
    push!(log_likelihoods, total_loglik)

    println(
        "R_0 = $(R_0), D_inf = $(D_inf) → ",
        "log-likelihood = $(round(total_loglik, digits = 1))"
    )
    plot!(
        p, sol, idxs = (0, 2),
        label = "R_0 = $(theta.R_0), D_inf = $(theta.D_inf)", linewidth = 2
    )
end

scatter!(
    p, epi1.time, epi1.obs, label = "Observed", markersize = 4, color = :black,
    markershape = :circle
)
plot!(p)

R_0 = 1.1, D_inf = 2.0 → log-likelihood = -5763.6
R_0 = 1.5, D_inf = 2.0 → log-likelihood = -2853.9
R_0 = 2.5, D_inf = 2.0 → log-likelihood = -160.9
R_0 = 4.0, D_inf = 2.0 → log-likelihood = -2456.2
R_0 = 2.5, D_inf = 1.0 → log-likelihood = -5495.7
R_0 = 2.5, D_inf = 4.0 → log-likelihood = -2821.5

TipCode explainer

• max.(sol[2, :], 0.0) - extracts the second row (the I compartment) from the ODE solution matrix, where : means “all columns” (all time points), and clamps away the tiny negative values a solver can return as the epidemic dies out. The clamp trades a crash for a quieter failure: Poisson(0.0) is valid, and scoring a day with observed cases against it returns -Inf, which zeroes that parameter set rather than raising an error.
• [Poisson(It) for It in I] - an array comprehension that creates a Poisson distribution for each expected value It in the vector I. If you know R, this is like lapply(I, function(It) ...) or purrr::map(I, ~ ...) — it applies an operation to each element and collects the results.
• Poisson(λ) - creates a Poisson distribution with rate parameter λ. Used here as our observation model.
• logpdf.(obs_dists, epi1.obs) - the . broadcasts the logpdf function element-wise, calculating log-likelihood for each observation against its corresponding Poisson distribution.

NoteComing from R?

R vectorises most functions implicitly; Julia requires an explicit . for element-wise operations (broadcasting).

# R — automatically vectorises over obs and lambda
dpois(obs, lambda = expected, log = TRUE)

# Julia — need the dot for element-wise operation
logpdf.(Poisson.(expected), obs)

Forgetting the . on a vector operation is a common gotcha for R users.

TipSolution

$R_0 = 2.5$ with $D_\mathrm{inf} = 2.0$ fits best, which is unsurprising once you know these data were simulated from parameters close to that.

$R_0 = 1.1$ barely produces an epidemic, so it cannot reach the observed peak at all. $R_0 = 4.0$ produces one that is too fast and too large, peaking well before the data do.

Changing $D_\mathrm{inf}$ at fixed $R_0$ changes the timescale rather than the final size, because $\beta = R_0 / D_\mathrm{inf}$ and $\nu = 1 / D_\mathrm{inf}$ move together. A shorter infectious period gives a faster, sharper epidemic; a longer one gives a slower, broader epidemic reaching a similar total.

On the last question: here the visual and numerical best fits agree, but they do not have to. The Poisson log-likelihood weights each day by how improbable the observation is under the prediction, so missing by 20 cases on a day expecting 15 loses far more log-likelihood than missing by 20 on a day expecting 300. A curve that looks worse to the eye can score better if it fits the high-count days closely.

The observation model

Our likelihood makes specific assumptions about how observations relate to the underlying epidemic:

1. Poisson noise: Each observation is drawn from a Poisson distribution centred on the model’s predicted prevalence. This implies the variance equals the mean - when we expect 100 cases, we expect variation of about ±10.

2. Conditional independence: Given the deterministic trajectory, observations at different time points are independent. The observation error at time $t$ doesn’t affect the error at time $t+1$.

The total log-likelihood is the sum across all time points: if the model predicts well at most times but poorly at one, that one bad fit drags down the total.

TipExercise: When do these assumptions break down?

Think about when these assumptions might be violated:

• When might the variance not equal the mean?
• When might observations not be independent?
• What other aspects of the observation process might you want to model?

TipSolution

Variance exceeds the mean whenever cases arrive in clumps: superspreading events, household clusters, or an outbreak investigation that finds a dozen cases at once. This is overdispersion, and it is the normal state of surveillance data rather than the exception. The observation models session replaces Poisson with a negative binomial for exactly this reason.

Independence fails for several everyday reasons. Weekend effects depress Saturday and Sunday reports and inflate Monday. Reporting delays smear one day’s infections across the following week, so consecutive observations share cases. Retrospective case-finding adds a block of old cases to whichever day the review happened.

Other things worth modelling: under-reporting, so that only a fraction of infections are ever observed; the delay between infection and report; and ascertainment that changes through an epidemic as testing capacity and case definitions shift.

A helper function

We’ll be calculating likelihoods repeatedly. Rather than copy-pasting the ODE simulation and Poisson likelihood code from above each time, let’s extract it into a reusable function:

function simulate_and_loglik(R_0, D_inf, data; init_state = [999, 1, 0], tspan = (0, 30))
    # Simulate ODE
    prob = ODEProblem(sir_ode!, init_state, tspan, [R_0, D_inf])
    sol = solve(prob, saveat = data.time)
    # max(..., 0.0) guards against tiny negative values from the ODE solver
    I = max.(sol[2, :], 0.0)

    # Calculate Poisson log-likelihood
    loglik = sum(logpdf.(Poisson.(I), data.obs))

    return (I = I, loglik = loglik)
end

This function simulates the SIR model and returns both the trajectory and the log-likelihood. We’ll use it throughout the rest of this session.

TipCode explainer: keyword arguments and named tuples

The semicolon ; in the function signature separates positional arguments (R_0, D_inf, data) from keyword arguments (init_state, tspan). Keyword arguments have defaults and must be passed by name: simulate_and_loglik(2.5, 2.0, epi1, init_state = [499, 1, 0]).

The function returns a named tuple (I = I, loglik = loglik), which you access with dot syntax: result.I or result.loglik. Named tuples are lightweight and immutable — think of them as a simple way to return multiple values with meaningful names.

NoteComing from R?

Keyword arguments with defaults work like R function arguments:

# R
simulate_and_loglik <- function(R_0, D_inf, data,
                                init_state = c(999, 1, 0),
                                tspan = c(0, 30)) {
  # ...
  list(I = I, loglik = loglik)  # named list
}
result$I  # access with $

Julia’s named tuples are similar to R’s named lists but lighter-weight. The main difference is the ; to separate keyword from positional arguments — R doesn’t make this distinction.

ImportantTake a break here

You’ve completed Parts 1 and 2: priors and likelihoods. Part 3 combines them into the posterior and covers grid search.

If you’re working through this in one sitting, stop here for 10-15 minutes. The concepts in Part 3 build directly on what you’ve learned, and a short break will help consolidate your understanding.

Part 3: Understanding the posterior

What is a posterior?

The posterior combines our prior beliefs with the data evidence using Bayes’ theorem:

$$p(\theta|\text{data}) \propto p(\text{data}|\theta) \times p(\theta)$$

• $p(\theta|\text{data})$: Posterior - what we believe about parameters after seeing data
• $p(\text{data}|\theta)$: Likelihood - how well parameters explain the data
• $p(\theta)$: Prior - what we believed before seeing data

In log space: $\log p(\theta|\text{data}) = \log p(\text{data}|\theta) + \log p(\theta) + \text{constant}$

We can combine our stored prior and likelihood values to calculate the posterior:

# Posterior = Prior + Likelihood (using our stored values)
println("Posterior breakdown using flat priors:")
for (i, theta) in enumerate(test_params)
    log_post_flat = log_priors_flat[i] + log_likelihoods[i]

    println("R_0=$(theta.R_0), D_inf=$(theta.D_inf)")
    println("  Prior: $(round(log_priors_flat[i], digits=1)), Likelihood: $(round(log_likelihoods[i], digits=1)), Posterior: $(round(log_post_flat, digits=1))")
end

Posterior breakdown using flat priors:
R_0=1.1, D_inf=2.0
  Prior: -5.5, Likelihood: -5763.6, Posterior: -5769.1
R_0=1.5, D_inf=2.0
  Prior: -5.5, Likelihood: -2853.9, Posterior: -2859.4
R_0=2.5, D_inf=2.0
  Prior: -5.5, Likelihood: -160.9, Posterior: -166.4
R_0=4.0, D_inf=2.0
  Prior: -5.5, Likelihood: -2456.2, Posterior: -2461.7
R_0=2.5, D_inf=1.0
  Prior: -5.5, Likelihood: -5495.7, Posterior: -5501.2
R_0=2.5, D_inf=4.0
  Prior: -5.5, Likelihood: -2821.5, Posterior: -2827.1

How priors influence the posterior

The following compares flat vs informed priors for all our test parameters:

# Compare flat vs informative priors using our stored values
println("Flat vs Informed posterior comparison:")
for (i, theta) in enumerate(test_params)
    # Calculate posteriors using stored values
    log_post_flat = log_priors_flat[i] + log_likelihoods[i]
    log_post_informed = log_priors_informed[i] + log_likelihoods[i]

    println("R_0=$(theta.R_0), D_inf=$(theta.D_inf)")
    println("  Flat posterior: $(round(log_post_flat, digits=1)), Informed posterior: $(round(log_post_informed, digits=1))")
    println("  Difference: $(round(log_post_informed - log_post_flat, digits=1))")
    println()
end

Flat vs Informed posterior comparison:
R_0=1.1, D_inf=2.0
  Flat posterior: -5769.1, Informed posterior: -5767.1
  Difference: 2.0

R_0=1.5, D_inf=2.0
  Flat posterior: -2859.4, Informed posterior: -2856.7
  Difference: 2.7

R_0=2.5, D_inf=2.0
  Flat posterior: -166.4, Informed posterior: -162.9
  Difference: 3.5

R_0=4.0, D_inf=2.0
  Flat posterior: -2461.7, Informed posterior: -2459.9
  Difference: 1.7

R_0=2.5, D_inf=1.0
  Flat posterior: -5501.2, Informed posterior: -5497.4
  Difference: 3.8

R_0=2.5, D_inf=4.0
  Flat posterior: -2827.1, Informed posterior: -2824.9
  Difference: 2.2

When we have strong prior beliefs that align with the data, the posterior gets a boost. When priors conflict with data, they pull the posterior away from the likelihood.

Systematic posterior exploration

Now that we understand the components of Bayesian inference, we can systematically explore the posterior by trying many parameter combinations and calculating their log-posterior values manually.

Grid search to find the best parameters

We calculate log-posterior values for a grid of parameter combinations:

# Define grid of parameter values to try
R_0_values = 1.5:0.2:4.0
D_inf_values = 1.0:0.2:3.0

# Store results
best_log_posterior = -Inf
best_params = nothing
results = []

for R_0 in R_0_values
    for D_inf in D_inf_values
        # Prior
        log_prior = logpdf(Uniform(1, 20), R_0) + logpdf(Uniform(1, 14), D_inf)

        # Likelihood (using our helper)
        result = simulate_and_loglik(R_0, D_inf, epi1)

        # Posterior = prior × likelihood
        log_post = log_prior + result.loglik
        push!(results, (R_0 = R_0, D_inf = D_inf, log_posterior = log_post))

        if log_post > best_log_posterior
            best_log_posterior = log_post
            best_params = (R_0 = R_0, D_inf = D_inf)
        end
    end
end

println("Best parameters: R_0 = $(best_params.R_0), D_inf = $(best_params.D_inf)")
println("Best log-posterior: $(round(best_log_posterior, digits = 1))")

Best parameters: R_0 = 2.3, D_inf = 2.0
Best log-posterior: -98.5

TipCode explainer

• 1.5:0.2:4.0 - creates a range starting at 1.5 in steps of 0.2, stopping before it passes 4.0. That gives 13 values for R_0, the last of which is 3.9.
• The nested for loops try every combination of R_0 and D_inf values, which is 13 × 11 = 143 combinations.
• For each combination, we calculate the log-posterior (prior + likelihood) and track the best one found so far.

NoteComing from R?

Julia’s start:step:stop syntax is equivalent to seq(start, stop, by = step) in R. Note that R’s : operator (e.g. 1:10) only increments by 1, whereas Julia’s : supports any step size.

# R
seq(1.5, 4.0, by = 0.2)

# Julia
1.5:0.2:4.0

Is this the result you expected?

Apart from finding the best parameters we can also use our results to visualise the whole posterior surface.

# Visualise the posterior surface
R_0_unique = unique(r.R_0 for r in results)
D_inf_unique = unique(r.D_inf for r in results)

# The loop above iterates R_0 outer, D_inf inner. Column-major reshape into
# (n_D_inf, n_R_0) then gives a matrix whose rows index D_inf and columns
# index R_0 — the orientation heatmap expects for (x = R_0, y = D_inf).
post_matrix = reshape(
    [r.log_posterior for r in results],
    length(D_inf_unique), length(R_0_unique)
)

heatmap(
    R_0_unique, D_inf_unique, post_matrix,
    xlabel = "R_0", ylabel = "D_inf (days)",
    title = "Log-Posterior Surface", color = :viridis
)

TipExercise: Read the posterior surface

Look at the heatmap.

• Whereabouts on the grid is the posterior mass, and what shape is it?
• Most of the grid is a single colour. What does that tell you?
• The grid is 143 combinations for two parameters. What happens to that number for a model with four?

TipWhat to look for

The peak sits at $R_0 = 2.3$, $D_\mathrm{inf} = 2.0$, with a log-posterior of $-98.5$. The rest of the grid runs down to $-10530$, so the surface spans over ten thousand log units. Only seven of the 143 combinations come within 100 of the peak and only one comes within 10, which is why almost the whole heatmap is one flat colour: the colour scale is being spent on a range in which the data have already ruled everything out.

The high region is a diagonal band rather than a blob. The best $D_\mathrm{inf}$ rises with $R_0$, from 1.0 days at $R_0 = 1.5$ to 3.0 days at $R_0 = 3.9$, because the data constrain how fast the epidemic grows and a larger $R_0$ can be partly offset by a longer infectious period. The two parameters are correlated, and the grid shows that correlation directly.

Note also that the best grid point is not the $R_0 \approx 2.5$ the data were simulated from. The grid steps in 0.2, and one realisation of Poisson noise is enough to move the peak by a step; the best-fitting parameters are not the true ones, which is the whole reason for wanting a posterior rather than a point.

Four parameters on the same grid would be roughly $143^2 \approx 20{,}000$ combinations, each an ODE solve, and most of them spent on the flat region that the colour scale already shows is worthless. That is the argument for the sampling methods in the rest of the course: they spend their effort where the posterior mass is.

Posterior predictive checking

Finding the best parameters is only half the job - we also need to check whether the model actually fits the data well. We do this by simulating data from the fitted model and comparing it to our observations.

# Get the trajectory for the best parameters
best_result = simulate_and_loglik(best_params.R_0, best_params.D_inf, epi1)

# Plot: deterministic trajectory + observation noise
p = plot(title = "Best Fit: Observation Noise Only", xlabel = "Time", ylabel = "Cases")

# Multiple realisations with Poisson observation noise (same parameters each time)
for i in 1:10
    pred_data = rand.(Poisson.(best_result.I))
    plot!(p, epi1.time, pred_data, alpha = 0.3, color = :blue, label = "")
end

# Deterministic trajectory
plot!(p, epi1.time, best_result.I, linewidth = 2, color = :red,
      label = "Model (R_0=$(best_params.R_0), D_inf=$(best_params.D_inf))")
scatter!(p, epi1.time, epi1.obs, label = "Observed", markersize = 4, color = :black)

The blue lines show variation from observation noise only - the underlying trajectory is identical each time.

Impact of parameter uncertainty

But we don’t know the true parameters - we have a distribution of plausible values. What happens when we also vary the parameters?

First, a check on the grid itself. The posterior probability of a grid point is proportional to $\exp$ of its log-posterior, normalised over the grid, so we can ask how many points carry the bulk of the mass.

# Subtracting the maximum first leaves the normalised weights unchanged and stops
# exp of a large negative number underflowing to zero
log_posts = [r.log_posterior for r in results]
post_weights = exp.(log_posts .- maximum(log_posts))
post_weights ./= sum(post_weights)

n_95 = count(cumsum(sort(post_weights, rev = true)) .< 0.95) + 1
println("$(n_95) of $(length(post_weights)) grid points carry 95% of the posterior mass")

1 of 143 grid points carry 95% of the posterior mass

One point out of 143 holds essentially all of it. The grid steps in units of 0.2 while the posterior is only about 0.03 wide, so the grid is roughly six posterior standard deviations coarse and cannot resolve the distribution at all. Everything the grid can tell us is “the answer is near 2.3”. The surface above does not show that starkly, because it plots the log posterior: the best cell beats the next one by about 14 log units, a factor of roughly a million, which the colour scale flattens into a single shade of yellow. The concentration only becomes visible once those values are exponentiated and normalised, which is what the cell above does.

To see the parameter uncertainty we have to refine the grid around the mode.

# Refine around the mode the coarse grid found, rather than around a window
# written in by hand: if the data or the model change, the window follows.
R_0_fine = (best_params.R_0 - 0.3):0.01:(best_params.R_0 + 0.3)
D_inf_fine = (best_params.D_inf - 0.3):0.01:(best_params.D_inf + 0.3)

fine_results = []
for R_0 in R_0_fine, D_inf in D_inf_fine
    log_prior = logpdf(Uniform(1, 20), R_0) + logpdf(Uniform(1, 14), D_inf)
    result = simulate_and_loglik(R_0, D_inf, epi1)
    push!(
        fine_results,
        (R_0 = R_0, D_inf = D_inf, log_posterior = log_prior + result.loglik)
    )
end

fine_log_posts = [r.log_posterior for r in fine_results]
fine_weights = exp.(fine_log_posts .- maximum(fine_log_posts))
fine_weights ./= sum(fine_weights)

fine_n_95 = count(cumsum(sort(fine_weights, rev = true)) .< 0.95) + 1
println("$(fine_n_95) of $(length(fine_weights)) refined grid points carry 95% of the mass")

# Draw parameter values from the posterior, so that each grid point is drawn as
# often as the posterior says it should be
draws = fine_results[rand(Categorical(fine_weights), 15)]

# Plot trajectories from the drawn parameter values
p_uncertainty = plot(
    title = "Parameter + Observation Uncertainty", xlabel = "Time", ylabel = "Cases"
)

for r in draws
    result = simulate_and_loglik(r.R_0, r.D_inf, epi1)
    pred_data = rand.(Poisson.(result.I))  # Add observation noise too
    plot!(p_uncertainty, epi1.time, pred_data, alpha = 0.3, color = :blue, label = "")
end

scatter!(
    p_uncertainty, epi1.time, epi1.obs, label = "Observed", markersize = 4,
    color = :black
)

88 of 3721 refined grid points carry 95% of the mass

TipCode explainer

• exp.(log_posts .- maximum(log_posts)) - converts log-posterior values back to a probability scale.
• Categorical(weights) - a distribution over the indices 1:length(weights), where index i has probability weights[i]. Drawing from it and indexing into fine_results is how you sample from a posterior that is defined on a grid.
• for R_0 in R_0_fine, D_inf in D_inf_fine - a single loop over both ranges, equivalent to the nested loops used for the coarse grid.

TipWhat’s the difference?

In the first plot we fixed the parameters at their best values and varied only the observation noise. In the second we draw the parameters from the posterior as well, so the spread accounts for our uncertainty about the parameters on top of the noise in the observations. The difference is small here, because 30 daily observations pin this two-parameter model down tightly. It grows with weaker data and with models that have more parameters to be uncertain about.

TipExercise: Posterior predictive checking

Does the real data look consistent with our simulated datasets? This is posterior predictive checking - a key model validation technique that we’ll explore in depth in the model checking session.

TipWhat to look for

Ask whether the observed series sits inside the bulk of the simulated ones or consistently along one edge. Systematically high or low means the model is missing something structural rather than merely being uncertain.

Then check the features you care about: is the peak the right height, and does it arrive at the right time? A model can track the shape and still get the timing wrong by several days, which matters a great deal if you are forecasting.

Finally, look at the spread. If the real data wander outside the envelope more often than the envelope’s nominal coverage allows, the observation model is too tight. With a Poisson likelihood that is the usual finding, and it is the diagnostic that motivates the negative binomial later in the course.

Why we need algorithmic approaches

We have seen how being able to simulate a model, in combination with prior densities and observation probabilities (expressed in the likelihood) enables us to fit this model to data and develop an understanding of parametric uncertainty. But we have done so using grid search, which has serious limitations:

# The curse of dimensionality
n_params = [2, 3, 4, 5, 6, 10]
n_values_per_param = 15

println("Grid search scaling problems:")
for n in n_params
    total_evals = n_values_per_param^n
    println("$(n) parameters: $(total_evals) evaluations")
end

Grid search scaling problems:
2 parameters: 225 evaluations
3 parameters: 3375 evaluations
4 parameters: 50625 evaluations
5 parameters: 759375 evaluations
6 parameters: 11390625 evaluations
10 parameters: 576650390625 evaluations

This shows that the computational cost of grid search, i.e. the number of simulations that we need to run, reaches very large numbers even at moderate levels of complexity.

We have already met the problem in miniature. The first grid, 143 points across a plausible range, could not resolve a posterior two parameters wide, and refining it took 3721 evaluations to describe a distribution that a sampler summarises in a few thousand draws whatever the dimension.

In addition, this kind of sampling is wasteful, because most parameter combinations will be far from the ones that lead to a good fit to the data. Almost all of those 3721 refined evaluations landed on parameter values with negligible posterior weight.

Lastly, it is difficult to establish whether the grid is fine enough (or important parameter regions might be in the “gaps”), and whether the limits of the grid are wide enough. Note that we only knew where to refine because the coarse grid had already told us roughly where the mode was.

Learning points

TipLearning points

• Priors express beliefs before seeing data - explore them through simulation (prior predictive checking)
• Likelihood quantifies how well parameters explain observed data, given an observation model
• Posterior combines prior and likelihood via Bayes’ theorem: $p(\theta|\text{data}) \propto p(\text{data}|\theta) \times p(\theta)$
• Observation models encode assumptions (e.g., Poisson noise, conditional independence) - think critically about when they’re violated
• Grid search can explore posteriors but doesn’t scale beyond a few parameters
• Posterior predictive checking validates whether fitted models generate realistic data

References

• Distributions.jl documentation - Probability distributions in Julia
• DifferentialEquations.jl - ODE solving in Julia

See the further reading page for course-wide references on Bayesian inference, infectious disease modelling, and Julia.

Next session

In MCMC sampling we will see how Markov chain Monte Carlo algorithms solve these problems by exploring the posterior automatically. We start with simple algorithms you can code yourself, then see how Turing.jl automates the same ideas with better algorithms.