Bayesian Statistics for Data Scientists – Python Guide 2026

Bayesian statistics offers a fundamentally different approach to inference: instead of asking “what is the probability of the data given a hypothesis?”, it asks “what is the probability of the hypothesis given the data?” This framework is more intuitive, handles small samples better, and naturally incorporates prior knowledge. This guide teaches Bayesian thinking and implementation in Python with PyMC.

Bayes’ Theorem

The foundation is Bayes’ theorem: P(hypothesis | data) = P(data | hypothesis) × P(hypothesis) / P(data). The prior P(hypothesis) encodes what you believe before seeing data. The likelihood P(data | hypothesis) is how probable the data is under that hypothesis. The posterior P(hypothesis | data) is your updated belief after seeing the data. This update process is the core of Bayesian inference — you start with a prior belief and revise it as evidence arrives.

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

# Example: estimating a coin's bias
# Prior: Beta(2, 2) — slightly informed, centred at 0.5
# Likelihood: Binomial — we observe 14 heads in 20 flips

heads, n = 14, 20
alpha_prior, beta_prior = 2, 2

# Beta-Binomial conjugate: posterior is also Beta
alpha_post = alpha_prior + heads
beta_post  = beta_prior + (n - heads)

theta = np.linspace(0, 1, 200)
prior     = stats.beta.pdf(theta, alpha_prior, beta_prior)
posterior = stats.beta.pdf(theta, alpha_post,  beta_post)
likelihood = stats.binom.pmf(heads, n, theta)

plt.figure(figsize=(10, 5))
plt.plot(theta, prior,     label=f'Prior Beta({alpha_prior},{beta_prior})', lw=2)
plt.plot(theta, likelihood/likelihood.max()*posterior.max(),
                           label='Likelihood (scaled)', lw=2, ls='--')
plt.plot(theta, posterior, label=f'Posterior Beta({alpha_post},{beta_post})', lw=2)
plt.axvline(heads/n, color='red', ls=':', label=f'MLE = {heads/n:.2f}')
plt.xlabel('Coin bias θ')
plt.legend()
plt.title('Bayesian Coin Bias Estimation')
plt.show()

posterior_mean = alpha_post / (alpha_post + beta_post)
ci_low, ci_high = stats.beta.ppf([0.025, 0.975], alpha_post, beta_post)
print(f'Posterior mean: {posterior_mean:.3f}')
print(f'95% Credible Interval: [{ci_low:.3f}, {ci_high:.3f}]')

PyMC for Probabilistic Programming

pip install pymc arviz

import pymc as pm
import arviz as az
import numpy as np

np.random.seed(42)

# Simulated data: true mu=10, sigma=2
true_mu    = 10
true_sigma = 2
data       = np.random.normal(true_mu, true_sigma, size=50)

with pm.Model() as model:
    # Priors
    mu    = pm.Normal('mu',    mu=0, sigma=20)
    sigma = pm.HalfNormal('sigma', sigma=10)

    # Likelihood
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=data)

    # MCMC sampling
    trace = pm.sample(2000, tune=1000, chains=4,
                      target_accept=0.9, random_seed=42)

# Summarise results
print(az.summary(trace, var_names=['mu', 'sigma']))
az.plot_posterior(trace, var_names=['mu', 'sigma'])
az.plot_trace(trace)
plt.tight_layout()
plt.show()

Bayesian A/B Testing

import pymc as pm
import numpy as np
import arviz as az

# Conversion data
control_visitors, control_conversions   = 1000, 120  # 12%
treatment_visitors, treatment_conversions = 1000, 145  # 14.5%

with pm.Model() as ab_model:
    # Priors — weakly informative Beta(1,1) = uniform
    p_control   = pm.Beta('p_control',   alpha=1, beta=1)
    p_treatment = pm.Beta('p_treatment', alpha=1, beta=1)

    # Likelihoods
    obs_control   = pm.Binomial('obs_control',   n=control_visitors,
                                 p=p_control, observed=control_conversions)
    obs_treatment = pm.Binomial('obs_treatment', n=treatment_visitors,
                                 p=p_treatment, observed=treatment_conversions)

    # Derived quantity: lift
    lift = pm.Deterministic('lift', p_treatment - p_control)
    rr   = pm.Deterministic('relative_lift', lift / p_control)

    trace = pm.sample(3000, tune=1000, chains=2, random_seed=42)

# Probability that treatment beats control
prob_better = (trace.posterior['lift'] > 0).mean().item()
expected_lift = trace.posterior['lift'].mean().item()
rel_lift = trace.posterior['relative_lift'].mean().item()

print(f'P(treatment > control): {prob_better:.1%}')
print(f'Expected absolute lift: {expected_lift:.4f}')
print(f'Expected relative lift: {rel_lift:.1%}')
print(az.summary(trace, var_names=['p_control', 'p_treatment', 'lift']))

Hierarchical (Multi-Level) Models

import pymc as pm
import numpy as np

# Sales data across 8 stores — some have very little data
stores    = 8
n_sales   = np.array([5, 200, 15, 80, 3, 150, 25, 60])
conv_rate = np.array([0.12, 0.15, 0.10, 0.14, 0.08, 0.16, 0.11, 0.13])
conversions = (n_sales * conv_rate).astype(int)

store_idx = np.arange(stores)

with pm.Model() as hierarchical_model:
    # Hyperpriors — global distribution of store rates
    mu_global    = pm.Beta('mu_global', alpha=2, beta=10)
    kappa_global = pm.HalfNormal('kappa_global', sigma=20)

    alpha = pm.Deterministic('alpha', mu_global * kappa_global)
    beta  = pm.Deterministic('beta',  (1 - mu_global) * kappa_global)

    # Store-level rates drawn from global distribution
    p_stores = pm.Beta('p_stores', alpha=alpha, beta=beta, shape=stores)

    # Likelihood
    obs = pm.Binomial('obs', n=n_sales, p=p_stores, observed=conversions)

    trace = pm.sample(2000, tune=1000, chains=2, random_seed=42)

import arviz as az
print(az.summary(trace, var_names=['p_stores', 'mu_global']))
# Stores with little data are "shrunk" toward the global mean —
# this is partial pooling and avoids overfitting small samples

Bayesian vs Frequentist

Frequentist statistics gives you p-values (“probability of data this extreme if null is true”) and confidence intervals (“if we repeated this experiment 100 times, 95 of the intervals would contain the true value”). Bayesian statistics gives you credible intervals (“there is 95% probability the parameter is in this range”) and direct probability statements about hypotheses. For business A/B testing, Bayesian is more intuitive: “There is an 87% probability that treatment B is better” is more actionable than “p=0.04”. For scientific publishing where replicability is paramount, frequentist conventions still dominate.

Conclusion

Bayesian statistics is especially valuable when sample sizes are small, when you need to incorporate prior domain knowledge, or when you need honest uncertainty quantification rather than binary significance thresholds. PyMC makes Bayesian modelling accessible — you describe your model in plain Python and let MCMC do the hard work of computing posteriors. Start with the coin-bias example to build intuition, then apply Bayesian A/B testing to your next experiment for more nuanced, decision-relevant results.