The Core of Bayesian Estimation

An Intuitive and Mathematical Deep Dive into Bayes’ Theorem and MCMC
Author

Xiaoge Zhang, PhD

Published

September 17, 2026

1 A Simple Thought Experiment

Suppose we draw 30 numbers from an unknown Normal distribution. Our goal is to “guess” the true mean of this distribution.

  • The Frequentist approach: Uses methods like Maximum Likelihood Estimation (MLE), Method of Moments (MOM), or Generalized Method of Moments (GMM) to find the parameter value based purely on the observed 30 numbers. If we draw another 30 numbers, the estimate changes based on the new sample.
  • The Bayesian approach: Also considers the likelihood of the data. However, instead of finding the single parameter value that maximizes the likelihood (like MLE), it combines the likelihood with a Prior to capture how the probability changes across all possible parameter values. Thus, the Bayesian approach yields a full distribution (the Posterior) of the parameter, rather than a single point estimate.

Almost all parameter estimations in statistics are variations of this simple game, including those used in HTA or HEOR. In Health Technology Assessment (HTA), this Bayesian philosophy has become the cornerstone of advanced evidence synthesis (like Network Meta-Analysis). It allows us to naturally propagate uncertainty through complex decision models, rather than relying on the frequentist concept of repeating experiments to infinity.

In this chapter, we will:

  1. Explain the fundamental equation of Bayesian estimation.
  2. Use a simple Normal distribution example to show the concept.
  3. Visualize how the Prior and Likelihood merge into the Posterior.
  4. Explain why we need Markov Chain Monte Carlo (MCMC) for complex models like NMA.

2 The Fundamental Equation

2.1 Starting From Conditional Probability

For any two events \(A\) and \(B\) with \(P(B) > 0\), conditional probability is defined as

\[P(A \mid B) = \frac{P(A \cap B)}{P(B)}\]

Once we know \(B\) has happened, we restrict attention to the outcomes in which \(B\) is true, and ask what share of them also have \(A\).

The same definition works in the other direction, \(P(B \mid A) = P(A \cap B) / P(A)\). Both describe the same joint event \(A \cap B\), so

\[P(A \cap B) = P(A \mid B)\,P(B) = P(B \mid A)\,P(A)\]

Dividing through by \(P(B)\) gives Bayes’ Theorem in its general form:

\[P(A \mid B) = \frac{P(B \mid A)\,P(A)}{P(B)}\]

2.2 From Events to Parameters and Data

Bayesian estimation applies the same identity with two substitutions:

  • \(A\) becomes the unknown parameter \(\theta\) (the hypothesis).
  • \(B\) becomes the observed data \(y\) (the evidence).

The direction we can compute is \(P(y \mid \theta)\): given a value of the parameter, how probable are the data? The direction we want is \(P(\theta \mid y)\): given the data, what should we believe about the parameter? Bayes’ Theorem connects the two:

\[P(\theta | y) = \frac{P(y | \theta) P(\theta)}{P(y)}\]

When \(\theta\) is continuous, probabilities become densities and the denominator becomes an integral over every possible value of \(\theta\):

\[P(y) = \int P(y \mid \theta)\,P(\theta)\,d\theta\]

Where:

  • \(P(\theta | y)\) is the Posterior: What we believe about the parameter \(\theta\) after seeing the data \(y\).
  • \(P(y | \theta)\) is the Likelihood: How likely the data \(y\) is, given a specific parameter \(\theta\).
  • \(P(\theta)\) is the Prior: What we believed about the parameter \(\theta\) before seeing the data.
  • \(P(y)\) is the Marginal Likelihood (or Evidence): A normalizing constant that ensures the posterior integrates to 1.

Since \(P(y)\) is a constant with respect to \(\theta\), we often write the relationship as a proportionality:

\[\text{Posterior} \propto \text{Likelihood} \times \text{Prior}\]

3 Minimal Case: Guessing the Mean (Without Analytical Solutions)

Let’s continue with our thought experiment. Suppose we draw 30 numbers from an unknown distribution. We want to find the distribution of the true mean \(\theta\).

3.1 The Data (Likelihood)

To keep this example as simple as possible, suppose we already know the true standard deviation of the distribution (let’s assume it is 1), and we are only trying to estimate the unknown mean \(\theta\).

The Likelihood function is the product of the Probability Density Functions (PDF) of observing each of the 30 data points \(y_i\):

\[L(\theta | y) = \prod_{i=1}^{30} f(y_i | \theta) = \prod_{i=1}^{30} \frac{1}{\sqrt{2\pi}} e^{-\frac{(y_i - \theta)^2}{2}}\]

This function allows us to explore different values of \(\theta\). We are looking for values of \(\theta\) that make the left-hand side (the likelihood) larger. Essentially, we are assuming that the data points we actually observed are more likely to occur (having a higher probability) under the true parameter value.

If we take the Frequentist approach, our task ends here. Since we have the actual values of \(y_i\), this becomes a simple mathematical problem of finding the maximum of a function. The value of \(\theta\) that maximizes this likelihood is what we traditionally call Maximum Likelihood Estimation (MLE).

A Bayesian, on the other hand, does not simply solve this equation for a single maximum point. Instead, they use the same core idea but ask a broader question: “What does the entire distribution of \(\theta\) look like to make the likelihood of our data as large as possible?” Thus, the Bayesian approach seeks the distribution of the parameter, rather than a single peak.

3.2 The Prior (Our Guess)

Suppose through some other channels (such as previous literature or domain knowledge), we believe that these 30 numbers are drawn from a distribution with a mean around 0 (e.g., \(N(0, 1)\)), or at least we guess the shape of the distribution of \(\theta\) is like \(N(0, 1)\). We express this belief as our Prior distribution.

Recall that the Likelihood function \(L(\theta | y)\) above is simply a product of 30 PDF values. According to Bayes’ Theorem (\(\text{Posterior} \propto \text{Likelihood} \times \text{Prior}\)), we now take that product and multiply it by the PDF of our Prior (here it takes the form of \(N(0, 1)\)), \(f_{prior}(\theta)\).

\[P(\theta | y) \propto \prod_{i=1}^{30} f(y_i | \theta) \times f_{prior}(\theta)\]

This forms an even larger product of probability densities. Our ultimate goal is to find what kind of distribution of \(\theta\) makes this combined PDF (this new joint probability) as large as possible. This is how the Prior “pulls” the final estimate towards itself.

So, what kind of \(\theta\) will make this joint density larger? It either needs to make the left part (the Likelihood) larger, or make the Prior density larger, or ideally both. These are the values of \(\theta\) that will make our joint probability (the Posterior) larger.

If we have no prior thoughts or guesses about this \(\theta\), we can also replace this Prior with a constant (like 1), or a Uniform Distribution, or a Normal Distribution with a very large variance. In that case, the Prior provides no “pull”, and the Bayesian estimate reduces to the Maximum Likelihood Estimate (MLE).

3.3 The Posterior (The Final Distribution)

Now, how do we actually find this group of \(\theta\) values that make the joint density large?

In this simple example, we can calculate it directly. But for complex real-world models (like Network Meta-Analysis), finding the perfect mathematical formula is impossible. This is where Markov Chain Monte Carlo (MCMC) comes in.

MCMC is a simulation method. Instead of solving equations, it simply “tries” different values of \(\theta\) to see which ones result in a higher joint density. In a formal expression, it takes a random walk through different values of \(\theta\). Crucially, it is designed to obtain those \(\theta\) values that are in regions where the joint density (Likelihood \(\times\) Prior) is higher. By collecting thousands of samples of \(\theta\) from this walk, we “map out” the entire shape of the Posterior distribution.

By multiplying the Prior and the Likelihood, we get this final distribution. We don’t need a single “best guess” number; we want the entire shape of this new distribution to understand our uncertainty.

4 Visualizing the Update in R: A Simple MCMC Simulation

Let’s write a simple Markov Chain Monte Carlo (MCMC) algorithm (the Metropolis algorithm) to “find” the posterior distribution, just like we described above. Its target is the product of the Likelihood and the Prior, and the algorithm takes random walks to map it out.

With a single parameter, that same product is simply a function of \(\theta\), and we can plot it like any curve: compute its value at many closely spaced values of \(\theta\), join the points, and rescale so the area under it is 1. Putting the curve and the MCMC samples on one set of axes shows what the sampler has been mapping.

Code 
library(ggplot2)
library(dplyr)
library(knitr)

# 1. Generate the 30 observed numbers from a "true" distribution
set.seed(42) # for reproducibility
true_theta <- 1.5 # the "hidden" true parameter
y <- rnorm(30, mean = true_theta, sd = 1)

# 2. Define the Target (Log-Posterior)
log_target <- function(theta) {
    # sum of log-densities over the 30 points
    log_like <- sum(dnorm(y, mean = theta, sd = 1, log = TRUE))
    # log-prior: guess = 0, SD = 1
    log_prior <- dnorm(theta, mean = 0, sd = 1, log = TRUE)
    return(log_like + log_prior) # log of Likelihood x Prior
}

# 3. Metropolis Algorithm
n_iter <- 10000
samples <- numeric(n_iter)
current_theta <- 0 # starting point of the random walk

set.seed(42)
for (i in 1:n_iter) {
    # Propose a new theta (random walk)
    proposed_theta <- current_theta + rnorm(1, 0, 0.5)

    # Calculate acceptance ratio (on log scale)
    log_acc_ratio <- log_target(proposed_theta) - log_target(current_theta)

    # Accept or reject
1    if (log(runif(1)) < log_acc_ratio) {
        current_theta <- proposed_theta
    }
    samples[i] <- current_theta
}

# 4. Plot the same target directly as a curve
theta_grid <- seq(0.5, 2.5, length.out = 2000)
d_theta <- diff(theta_grid)[1]
log_post_grid <- sapply(theta_grid, log_target)
2h <- exp(log_post_grid - max(log_post_grid))
3post_dens <- h / sum(h * d_theta)

# 5. Plot the MCMC samples against the directly drawn curve
ggplot() +
    geom_histogram(
        data = data.frame(theta = samples), aes(x = theta, y = after_stat(density)),
        binwidth = 0.05, boundary = 0, fill = "lightblue", color = "white"
    ) +
    geom_line(data = data.frame(theta = theta_grid, density = post_dens), aes(x = theta, y = density), color = "blue", linewidth = 1) +
    coord_cartesian(xlim = c(0.8, 2.3)) +
    theme_minimal() +
    labs(
        title = "The same posterior, sampled and drawn",
        subtitle = "Bars: where the 10,000 MCMC draws fell. Line: Likelihood x Prior plotted directly",
        x = "True Mean (Theta)",
        y = "Density"
    )
1
Random threshold to allow “downhill” moves proportional to probability
2
Subtract the maximum before exponentiating, so the values do not underflow to zero
3
Rescale to area 1: the numerical counterpart of \(P(y)\)

Code 
# 6. Summarise both
probs <- c(0.025, 0.5, 0.975)
grid_mean <- sum(theta_grid * post_dens) * d_theta
grid_sd <- sqrt(sum((theta_grid - grid_mean)^2 * post_dens) * d_theta)
grid_q <- sapply(probs, function(p) theta_grid[which(cumsum(post_dens * d_theta) >= p)[1]])

summary_tbl <- rbind(
    "MCMC samples" = c(mean(samples), sd(samples), quantile(samples, probs)),
    "Plotted directly" = c(grid_mean, grid_sd, grid_q)
)
kable(summary_tbl, digits = 3, col.names = c("Mean", "SD", "2.5%", "Median", "97.5%"))
Mean SD 2.5% Median 97.5%
MCMC samples 1.520 0.176 1.167 1.519 1.854
Plotted directly 1.518 0.180 1.166 1.518 1.870

Each bar covers an interval of \(\theta\) 0.05 wide. Its height is the share of the 10,000 MCMC draws that fell in that interval, divided by 0.05 so that the bars have a total area of 1, like the area under the line. A tall bar marks a range of \(\theta\) where the random walk spent much of its time.

The bars and the line are one distribution reached two ways: the MCMC by a random walk, the line by computing the formula along the whole axis. Plotting works here only because there is one parameter: at 100 values per parameter, ten parameters would need \(100^{10}\) calculations. The MCMC computes \(\text{Likelihood} \times \text{Prior}\) only at the points it visits, so its cost does not explode as parameters are added. That is why it can handle Network Meta-Analysis models, which carry dozens of parameters (see the next section).

Note: With 30 samples, the data provides stronger evidence. If you reduce the sample size to 10, you will see the posterior peak drift more due to random sampling error, as the likelihood becomes less sharp and the prior at 0 exerts a relatively stronger pull.

4.1 With Half the Data

What happens with only 15 observations? To separate the effect of sample size from the luck of one particular sample, the comparison is run on two independent draws from the same true distribution (\(\theta = 1.5\)). Draw 1 is the 30 numbers used above; Draw 2 is a fresh set of 30. In each draw, the 15-observation case uses the first 15 of its 30 numbers. The prior and the algorithm are unchanged.

Code 
# The same Metropolis algorithm as above, written as a function of the data
run_metropolis <- function(y_obs) {
    log_target_obs <- function(theta) {
        sum(dnorm(y_obs, mean = theta, sd = 1, log = TRUE)) + dnorm(theta, mean = 0, sd = 1, log = TRUE)
    }
    draws <- numeric(n_iter)
    current_theta <- 0
    set.seed(42)
    for (i in 1:n_iter) {
        proposed_theta <- current_theta + rnorm(1, 0, 0.5)
        if (log(runif(1)) < log_target_obs(proposed_theta) - log_target_obs(current_theta)) {
            current_theta <- proposed_theta
        }
        draws[i] <- current_theta
    }
    draws
}

set.seed(43)
y2 <- rnorm(30, mean = true_theta, sd = 1) # a second, independent sample of 30

runs <- list(
    list(draw = "Draw 1", n = 30, y_obs = y),
    list(draw = "Draw 1", n = 15, y_obs = y[1:15]),
    list(draw = "Draw 2", n = 30, y_obs = y2),
    list(draw = "Draw 2", n = 15, y_obs = y2[1:15])
)

df_runs <- bind_rows(lapply(runs, function(r) {
    data.frame(draw = r$draw, n = paste(r$n, "observations"), theta = run_metropolis(r$y_obs))
}))
df_runs$n <- factor(df_runs$n, levels = c("30 observations", "15 observations"))

cri <- df_runs %>%
    group_by(n, draw) %>%
    summarise(
        Mean = mean(theta),
        SD = sd(theta),
        Lower = quantile(theta, 0.025),
        Upper = quantile(theta, 0.975),
        .groups = "drop"
    ) %>%
    mutate(Width = Upper - Lower)

ggplot(df_runs, aes(x = theta)) +
    geom_histogram(aes(y = after_stat(density)), binwidth = 0.05, boundary = 0, fill = "lightblue", color = "white") +
    geom_vline(xintercept = true_theta, color = "grey40") +
    geom_vline(data = cri, aes(xintercept = Lower), color = "blue", linetype = "dashed") +
    geom_vline(data = cri, aes(xintercept = Upper), color = "blue", linetype = "dashed") +
    geom_text(data = cri, aes(x = Lower, y = Inf, label = sprintf("%.2f", Lower)), color = "blue", size = 4, hjust = 1.1, vjust = 1.5) +
    geom_text(data = cri, aes(x = Upper, y = Inf, label = sprintf("%.2f", Upper)), color = "blue", size = 4, hjust = -0.1, vjust = 1.5) +
    facet_wrap(~ n + draw, ncol = 1, labeller = label_wrap_gen(multi_line = FALSE)) +
    scale_x_continuous(breaks = seq(0.5, 2.5, by = 0.5)) +
    coord_cartesian(xlim = c(0.5, 2.7)) +
    theme_minimal() +
    labs(
        title = "Half the data: a wider interval and a less reliable centre",
        subtitle = "Dashed lines: 95% credible interval. Grey line: true value (1.5)",
        x = "True Mean (Theta)",
        y = "Density"
    )

Code 
kable(cri, digits = 3, col.names = c("Sample", "Draw", "Mean", "SD", "95% CrI lower", "95% CrI upper", "Width"))
Sample Draw Mean SD 95% CrI lower 95% CrI upper Width
30 observations Draw 1 1.520 0.176 1.167 1.854 0.687
30 observations Draw 2 1.325 0.174 0.981 1.655 0.673
15 observations Draw 1 1.857 0.255 1.374 2.359 0.985
15 observations Draw 2 1.218 0.251 0.742 1.711 0.968

Two things change when the sample is halved.

The interval widens, by about 40%, and by the same amount in both draws. The reason is how evidence accumulates. Each observation is a noisy clue about \(\theta\): on its own, a single number could easily sit a whole unit away from it. A candidate value of \(\theta\) has to fit every clue at once, and the more clues there are, the harder it is for a value far from the truth to fit them all, so fewer candidates survive. Formally, each observation adds a term \(-(y_i - \theta)^2/2\) to the log-likelihood: a downward parabola in \(\theta\) that bends by the same amount whatever the value of \(y_i\). Thirty of them add up to a parabola 30 times as steep, fifteen to one only half as steep, so the log-likelihood falls away more slowly as \(\theta\) moves from the peak and a wider range of \(\theta\) stays plausible, spreading the posterior out.

In this Normal example the arithmetic is exact. The posterior SD is \(1/\sqrt{n + 1}\), where the \(+1\) is the prior’s contribution, so halving the data widens the interval by \(\sqrt{31/16} \approx 1.39\). It depends only on how many observations there are, not on their values, which is why both draws widen alike. More generally, an interval narrows only with \(\sqrt{n}\): halving the width takes four times the data.

The centre becomes less reliable. With 30 observations the two draws are centred at about 1.5 and 1.3, 0.2 apart. With 15 they are centred at about 1.9 and 1.2, more than 0.6 apart. A smaller sample’s mean wanders further from the truth, and that is exactly the uncertainty the wider interval describes. Two draws are an illustration rather than a proof: the standard error of a sample mean is \(1/\sqrt{n}\), which is 0.18 for 30 observations and 0.26 for 15.

5 Moving Beyond Simple Examples: Why we need MCMC in NMA

In our simple example above, we could actually have calculated the posterior directly using math formulas (since the Normal distribution with a Normal prior has a known closed-form solution).

But in real-world Network Meta-Analysis (NMA), the likelihoods are much more complex (e.g., Logistic regression for binary outcomes, or Cox Proportional Hazards for survival data) and there are many shared parameters across multiple trials. The denominator \(P(y) = \int P(y|\theta)P(\theta) d\theta\) becomes a high-dimensional integral that is impossible to solve analytically.

This is why we use advanced MCMC software like JAGS or Stan. It is important to realize that they are doing exactly the same thing we demonstrated in our simple R code above. The underlying logic remains unchanged: evaluating the relative height of the numerator (\(\text{Likelihood} \times \text{Prior}\)) to decide which parameter values to keep or discard, thereby mapping the distribution.

The primary difference lies in their search and iteration methods (how they propose the next step):

  • JAGS relies heavily on Gibbs Sampling (and Metropolis-Hastings). It typically updates one parameter at a time. This is like taking random walks along the axes of the parameter space. While effective for simpler models, it can be slow and struggle to converge in complex, highly correlated NMA models.
  • Stan uses Hamiltonian Monte Carlo (HMC). Instead of taking blind random walks, it calculates the gradient of the log-posterior (the slope of the probability hill) and uses physical simulation concepts to glide efficiently through the parameter space.

Both tools are simply more sophisticated “walkers” designed to scale up to the dozens or hundreds of parameters required in professional HTA submissions.

6 Audit Point: The “Prior” in HTA Submissions

In NICE submissions, the choice of priors is heavily scrutinized:

  1. Vague/Non-informative Priors: Often used for treatment effects to “let the data speak”.
  2. Informative Priors: Used when data is extremely sparse (e.g., rare diseases). We might “borrow” a prior for the between-study heterogeneity variance (\(\tau^2\)) from a historical database of similar trials (e.g., the Turner or Rhodes priors).

A credible Bayesian analyst must always perform a Prior Sensitivity Analysis to prove that the conclusions are driven by the trial data, not by arbitrary choices of the prior.