Bayesian Statistics for Data Science: A Practical Introduction (2026)

Most data scientists learn statistics through the frequentist lens — p-values, confidence intervals, hypothesis tests. But there’s another entire framework for statistical reasoning: Bayesian statistics. Rather than asking “how often would this result occur by chance?”, Bayesian statistics asks “given what I’ve observed, what should I believe?” This shift in framing often produces more intuitive answers and richer insights, especially when you have prior knowledge or limited data.

Bayes’ Theorem: The Foundation

Bayes’ theorem relates prior beliefs to posterior beliefs after observing data. In its probability form: P(hypothesis | data) = P(data | hypothesis) × P(hypothesis) / P(data). In plain English: your updated belief equals your likelihood of seeing this data if the hypothesis is true, multiplied by your prior belief, normalised by the probability of seeing this data at all.

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# Classic example: Is this coin fair?
# Prior: believe the coin is fair (bias = 0.5) but uncertain
# We observe 7 heads out of 10 flips

# Grid approximation — simple Bayesian inference
bias_grid = np.linspace(0, 1, 1000)  # possible values of P(heads)

# Prior: uniform (all biases equally likely before seeing data)
prior = np.ones(len(bias_grid))
prior /= prior.sum()

# Likelihood: P(7 heads in 10 flips | bias) — binomial distribution
n_flips, n_heads = 10, 7
likelihood = stats.binom.pmf(n_heads, n_flips, bias_grid)

# Posterior = likelihood * prior (then normalise)
posterior = likelihood * prior
posterior /= posterior.sum()

# Summary statistics
posterior_mean = np.sum(bias_grid * posterior)
# 95% credible interval
cumulative = np.cumsum(posterior)
ci_lower = bias_grid[np.searchsorted(cumulative, 0.025)]
ci_upper = bias_grid[np.searchsorted(cumulative, 0.975)]

print(f"Posterior mean:        {posterior_mean:.3f}")
print(f"95% Credible Interval: [{ci_lower:.3f}, {ci_upper:.3f}]")
print("Interpretation: There is a 95% probability the true bias")
print(f"lies between {ci_lower:.2f} and {ci_upper:.2f}")

plt.figure(figsize=(10, 4))
plt.plot(bias_grid, prior, 'b--', label='Prior', alpha=0.6)
plt.plot(bias_grid, likelihood / likelihood.sum(), 'g-', label='Likelihood (scaled)')
plt.plot(bias_grid, posterior, 'r-', label='Posterior', linewidth=2)
plt.axvline(posterior_mean, color='red', linestyle=':', label=f'Posterior mean = {posterior_mean:.3f}')
plt.xlabel('Coin bias (P(heads))'); plt.ylabel('Probability density')
plt.title('Bayesian Coin Flip Analysis'); plt.legend()
plt.show()

Notice the key difference from frequentist output: the 95% credible interval has a direct probabilistic interpretation — “there is a 95% probability the parameter lies in this range.” A frequentist 95% confidence interval does NOT mean this (it means 95% of such intervals constructed this way would contain the true value).

Bayesian A/B Testing

Bayesian A/B testing is one of the most practically useful applications for product and growth data scientists. Unlike frequentist A/B tests (which require choosing sample size in advance and can’t be “peeked” at), Bayesian tests can be monitored continuously and give intuitive answers like “there’s an 87% probability that variant B is better than A”:

from scipy.stats import beta

def bayesian_ab_test(control_conversions, control_trials,
                     variant_conversions, variant_trials,
                     n_samples=100_000):
    # Bayesian A/B test for conversion rates.
    # Prior: Beta(1, 1) = uniform prior (no prior knowledge)
    # Posterior distributions
    # Beta(alpha, beta) is conjugate prior for binomial likelihood
    control_posterior = beta(
        a=1 + control_conversions,
        b=1 + (control_trials - control_conversions)
    )
    variant_posterior = beta(
        a=1 + variant_conversions,
        b=1 + (variant_trials - variant_conversions)
    )

    # Monte Carlo: sample from both posteriors
    control_samples = control_posterior.rvs(n_samples)
    variant_samples = variant_posterior.rvs(n_samples)

    # P(variant > control)
    prob_variant_better = (variant_samples > control_samples).mean()

    # Expected lift
    expected_lift = (variant_samples / control_samples - 1).mean() * 100

    # Credible interval on the difference
    diff_samples = variant_samples - control_samples
    ci = np.percentile(diff_samples, [2.5, 97.5])

    print(f"Control:  {control_conversions}/{control_trials} = {control_conversions/control_trials:.2%}")
    print(f"Variant:  {variant_conversions}/{variant_trials} = {variant_conversions/variant_trials:.2%}")
    print(f"P(variant > control): {prob_variant_better:.1%}")
    print(f"Expected lift: {expected_lift:.1f}%")
    print(f"95% CI on difference: [{ci[0]:.3f}, {ci[1]:.3f}]")

    return prob_variant_better

# Example: Button colour test
prob = bayesian_ab_test(
    control_conversions=450, control_trials=5000,
    variant_conversions=510, variant_trials=5000
)

Bayesian Inference with PyMC

For more complex models, PyMC provides a full probabilistic programming framework. You describe the model in terms of distributions, and PyMC uses Markov Chain Monte Carlo (MCMC) to sample from the posterior:

import pymc as pm
import numpy as np

# Bayesian linear regression
np.random.seed(42)
X = np.random.normal(0, 1, 100)
y = 2.5 * X + 0.8 + np.random.normal(0, 0.5, 100)  # true: slope=2.5, intercept=0.8

with pm.Model() as linear_model:
    # Priors
    slope     = pm.Normal('slope',     mu=0, sigma=10)
    intercept = pm.Normal('intercept', mu=0, sigma=10)
    sigma     = pm.HalfNormal('sigma', sigma=1)

    # Likelihood
    mu = slope * X + intercept
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

    # Sample from posterior
    trace = pm.sample(2000, tune=1000, return_inferencedata=True, progressbar=False)

print(pm.summary(trace, var_names=['slope', 'intercept', 'sigma']))

Frequently Asked Questions

When should I use Bayesian statistics instead of frequentist?
Bayesian methods shine when: you have genuine prior knowledge to incorporate, you need to update beliefs incrementally as data arrives, you want direct probability statements about parameters, or your sample size is small. Frequentist methods work well for large samples and when you want results that don’t depend on prior specification.

What’s a conjugate prior?
A conjugate prior is one where the prior and posterior belong to the same distribution family. For example, a Beta prior with a binomial likelihood gives a Beta posterior. This allows analytical (no sampling required) computation of the posterior — which is why the Beta distribution is used for A/B testing conversion rates.

Is Bayesian inference always computationally expensive?
Not always. Conjugate models have closed-form solutions. Grid approximation works for one or two parameters. MCMC is needed for complex multi-parameter models and can be slow but is highly parallelisable. Variational inference (used in PyMC and TensorFlow Probability) is faster but approximate.