%matplotlib widget
Bayesian analysis¶
The previous tutorials demonstrate how to obtain maximum-likelihood estimates of the parameters of a model by minimising $\chi^2$. An alternative, and often more informative, view of the same problem is provided by a Bayesian analysis, in which the goal is to characterise the full posterior distribution over the parameters given the data.
Recall Bayes' theorem,
$$ p(\theta \mid d) \propto p(d \mid \theta)\, p(\theta), $$ (bayes)
where $\theta$ are the model parameters, $d$ is the observed data, $p(d \mid \theta)$ is the likelihood, and $p(\theta)$ is the prior. In easyscience, the min/max bounds of a Parameter are interpreted as a uniform prior, and a Gaussian likelihood is constructed from the data and supplied weights.
easyscience exposes a Bayesian Markov-chain Monte Carlo (MCMC) sampler through the Fitter.mcmc_sample method. Under the hood this uses BUMPS' DREAM sampler, so the underlying minimizer must be switched to BUMPS.
{note}
This tutorial focuses on Bayesian analysis with a simple QENS model for illustration. For dedicated QENS fitting with more sophisticated models, consider using [`EasyDynamics`](https://github.com/easyscience/easydynamics).
When should you use a Bayesian analysis?¶
A maximum-likelihood estimate (MLE) gives you a single best-fit value with a symmetric uncertainty, which is fast and often sufficient. A Bayesian analysis becomes valuable when:
- Your uncertainties are asymmetric — MLE error bars assume the parameter distribution is Gaussian, which is not always true. The posterior samples capture skew naturally.
- You have prior knowledge — if you know from physics that a parameter must lie in a certain range, encoding that as a prior ($p(\theta)$) is more principled than simply clamping bounds after the fit.
- You care about parameter correlations — the joint posterior (shown in the corner plot below) reveals trade-offs between parameters that a single covariance matrix can miss.
- You want to propagate uncertainty to predictions — with the posterior in hand, you can compute credible bands on any function of the parameters (see the posterior-predictive band section) without linearised error propagation.
The trade-off is computational cost: MCMC requires thousands of model evaluations, whereas an MLE fit may converge in dozens. For the simple 4-parameter model used here the difference is negligible, but for expensive models it is worth starting with MLE and only switching to MCMC when you need the richer output.
In this tutorial we re-use the QENS dataset and the Lorentzian-with-resolution model from the Fitting QENS tutorial, but instead of returning a single best-fit value with a symmetric error bar we will draw thousands of samples from the posterior.
Load the data¶
We load the same simulated QENS dataset that was used in the previous tutorial.
import matplotlib.pyplot as plt
import numpy as np
def fetch_data(name: str, known_hash: str) -> str:
"""
Fetch pre-prepared data from a remote source and return the path to the file.
"""
import pooch
return pooch.retrieve(
url=f'https://public.esss.dk/groups/scipp/dmsc-summer-school/2025/{name}',
known_hash=known_hash,
)
filename = fetch_data(
'4-reduction/energy_transfer_QENS_unknown_quasi_elastic_many_neutrons.dat',
known_hash='sha256:e49fa9a1d2ef5eeb714524903e8ffa9de6616e5a299a799cbb0b2f3ee8fd459a',
)
def load(filename: str):
"""Load three-column text data (x, y, error) and filter NaN values.
The file is expected to be a space-separated three-column text file
(commonly ``.dat`` or ``.txt``) with columns: x, y, error (1σ).
"""
x, y, error = np.loadtxt(filename, unpack=True)
selection = np.isfinite(y)
return x[selection], y[selection], error[selection]
omega, intensity_obs, intensity_error = load(filename)
# Restrict to the region of interest, as in the QENS tutorial.
selection = (omega > -0.06) & (omega < 0.06)
omega, intensity_obs, intensity_error = (
omega[selection],
intensity_obs[selection],
intensity_error[selection],
)
Downloading data from 'https://public.esss.dk/groups/scipp/dmsc-summer-school/2025/4-reduction/energy_transfer_QENS_unknown_quasi_elastic_many_neutrons.dat' to file '/home/runner/.cache/pooch/13c45bd8472df2cbc0ff550e2ec9ed15-energy_transfer_QENS_unknown_quasi_elastic_many_neutrons.dat'.
fig, ax = plt.subplots()
ax.errorbar(omega, intensity_obs, intensity_error, fmt='.')
ax.set(xlabel='$\\omega$/meV', ylabel='$I(\\omega)$')
plt.show()
Defining parameters with priors¶
Create four Parameter objects, for the area $A$, $\gamma$, $\omega_0$ and $\sigma$. The min and max arguments define a uniform prior on each parameter — the sampler will only consider values inside this range and will treat every value inside the range as equally plausible a priori.
| Parameter | Initial Value | Min | Max |
|---|---|---|---|
| $A$ (area) | 10 | 1 | 100 |
| $\gamma$ | 8.0 × 10-3 | 1.0 × 10-4 | 1.0 × 10-2 |
| $\omega_0$ | 1.0 × 10-3 | 0 | 2.0 × 10-3 |
| $\sigma$ | 1.0 × 10-3 | 1.0 × 10-5 | 1.0 × 10-1 |
from easyscience import Parameter
area = Parameter(name='area', value=10, fixed=False, min=1, max=100)
gamma = Parameter(name='gamma', value=8e-3, fixed=False, min=1e-4, max=1e-2)
omega_0 = Parameter(name='omega_0', value=1e-3, fixed=False, min=0, max=2e-3)
sigma = Parameter(name='sigma', value=1e-3, fixed=False, min=1e-5, max=1e-1)
The model¶
We re-use the convolution of a Lorentzian with a Gaussian resolution function from the QENS tutorial:
$$ I(\omega) = \frac{A\gamma}{\pi\big[(\omega - \omega_0)^2 + \gamma^2\big]} \;\ast\; \mathcal{N}(0, \sigma), $$ (model)
where $A$ is a scale factor (area), $\gamma$ is the Lorentzian half-width at half-maximum, $\omega_0$ is the centre offset and $\sigma$ is the width of the Gaussian resolution kernel.
You might want to look at other QENS models, as implemented in the EasyDynamics library.
{note}
The convolution below uses `np.convolve(..., 'same')` for simplicity. In `'same'` mode, `numpy` pads the signal edges with zeros, which can introduce small artefacts at the boundaries. Because our Gaussian kernel ($\sigma \approx 10^{-3}$) is much narrower than the data range ($\pm 0.06$ meV), these edge effects are negligible here. For production work with broader kernels, consider `scipy.signal.convolve` with an explicit boundary mode such as `'reflect'` or `'nearest'`.
from scipy.stats import norm
def lorentzian(x: np.ndarray) -> np.ndarray:
return area.value / np.pi * gamma.value / ((x - omega_0.value) ** 2 + gamma.value**2)
def intensity_model(x: np.ndarray) -> np.ndarray:
gauss = norm(0, sigma.value).pdf(x)
gauss /= gauss.sum()
return np.convolve(lorentzian(x), gauss, 'same')
Maximum-likelihood fit (for reference)¶
Before running MCMC, perform a quick maximum-likelihood fit. The result will give us a good starting point and lets us compare the central tendency of the posterior with the classical best-fit values.
from easyscience import Fitter
from easyscience import ObjBase
parameter_container = ObjBase(name='params', A=area, gamma=gamma, omega_0=omega_0, sigma=sigma)
mle_fitter = Fitter(parameter_container, intensity_model)
mle_result = mle_fitter.fit(x=omega, y=intensity_obs, weights=1 / intensity_error)
print(f'A = {area.value:.4g}')
print(f'gamma = {gamma.value:.4g}')
print(f'omega_0 = {omega_0.value:.4g}')
print(f'sigma = {sigma.value:.4g}')
ObjBase is deprecated and will be removed in a future version. Please migrate to ModelBase.
A = 5.227 gamma = 0.007161 omega_0 = 0.001142 sigma = 0.001087
Drawing posterior samples with DREAM¶
We now draw samples from the posterior distribution $p(\theta \mid d)$ using the BUMPS DREAM (DiffeRential Evolution Adaptive Metropolis) algorithm. DREAM is an ensemble MCMC method that runs multiple chains in parallel and automatically tunes the proposal distribution.
DREAM only works with the BUMPS minimizer. We reuse the mle_fitter created above, switch it to BUMPS, and call Fitter.mcmc_sample, which returns a dictionary with the following keys:
draws: a(n_samples, n_parameters)array of posterior samples: each row is one complete draw from the joint posterior (one value for every parameter simultaneously), and each column holds all sampled values for a single parameter;param_names: the unique names of the parameters, in the same column order asdraws;internal_bumps_object: the underlying BUMPSMCMCDrawobject, useful for advanced diagnostics;logp: the log-posterior of each retained sample.
The key sampling parameters are:
samples(4000): the total number of posterior draws to generate across all chains;burn(500): the number of initial burn-in iterations to discard — the sampler needs time to find the typical set of the posterior, and early samples are not representative;thin(2): the thinning interval — only every second sample is kept, which reduces autocorrelation between consecutive draws;
First, we switch to the BUMPS minimizer:
from easyscience import AvailableMinimizers
mle_fitter.switch_minimizer(AvailableMinimizers.Bumps)
result = mle_fitter.mcmc_sample(
x=omega,
y=intensity_obs,
weights=1 / intensity_error,
samples=10000,
burn=500,
thin=2,
)
print(f'Drew {result["draws"].shape[0]} samples for {result["draws"].shape[1]} parameters.')
print('parameters:', result['param_names'])
Drew 5000 samples for 4 parameters. parameters: ['Parameter_0', 'Parameter_1', 'Parameter_2', 'Parameter_3']
Convergence diagnostics¶
Before we trust the posterior samples, we should check that the MCMC chains have converged — that is, the sampler has found the typical set of the posterior and is no longer drifting. Two simple visual checks are:
- Trace plot — plot the sampled parameter values against the sample index. A well-converged chain looks like a "hairy caterpillar": it fluctuates around a stable mean with no long-term trends.
- Log-posterior plot — the log-posterior $\log p(\theta \mid d)$ should also stabilise after burn-in. If it is still climbing at the end of the run, the sampler has not yet converged.
The logp array returned by sample contains the log-posterior for every retained sample (i.e. after burn-in and thinning). This means a flat logp trace is exactly what we want to see.
draws = result['draws']
logp = result['logp']
if callable(logp):
_, logp = logp()
logp = logp.flatten()
name_to_col = {name: idx for idx, name in enumerate(result['param_names'])}
def column_for(parameter):
return draws[:, name_to_col[parameter.unique_name]]
fig, axes = plt.subplots(5, 1, figsize=(10, 12), sharex=True)
# Trace plots for each parameter
for ax, (label, par) in zip(
axes[:4],
(('area', area), ('gamma', gamma), ('omega_0', omega_0), ('sigma', sigma)),
):
ax.plot(column_for(par), lw=0.5)
ax.set_ylabel(label)
ax.set_xlim(0, len(draws) - 1)
# Log-posterior trace
axes[4].plot(logp, lw=0.5, color='C4')
axes[4].set_ylabel('log-posterior')
axes[4].set_xlabel('sample index')
fig.suptitle('MCMC trace plots — check for "hairy caterpillar" behaviour')
fig.tight_layout()
plt.show()
Posterior summaries¶
Summarise each marginal posterior by its median and the 16th/84th percentiles, which together give an asymmetric 68% credible interval. Compare these with the MLE values you obtained above.
Why might the Bayesian median differ from the MLE? The MLE finds the single point that maximises the likelihood, while the Bayesian median is the central value of the posterior — which also accounts for the prior $p(\theta)$. If a parameter's posterior is skewed (asymmetric), the median and the mode (which approximates the MLE) will not coincide. The table below shows both so you can spot any such differences.
Note that the columns of result['draws'] are ordered by result['param_names'] (which use the parameters' unique_name), so it is worth building a small helper to look up a column by friendly name.
summary_rows = []
for label, par in (
('area', area),
('gamma', gamma),
('omega_0', omega_0),
('sigma', sigma),
):
col = column_for(par)
lo, med, hi = np.percentile(col, [16, 50, 84])
summary_rows.append((label, med, med - lo, hi - med, par.value))
print(f'{"param":<8s} {"median":>12s} {"16th":>12s} {"84th":>12s} {"MLE":>12s}')
for label, med, low, high, mle in summary_rows:
print(f'{label:<8s} {med:12.4g} {low:12.4g} {high:12.4g} {mle:12.4g}')
param median 16th 84th MLE area 5.223 0.05534 0.05122 5.29 gamma 0.007143 0.000158 0.0001391 0.007199 omega_0 0.001132 0.0001028 9.926e-05 0.001125 sigma 0.001163 0.0005375 0.0005524 0.001176
Visualise the joint posterior¶
Marginal summaries hide correlations between parameters. A corner plot (a triangular grid of pairwise scatter plots and 1-D histograms) is the standard way to display them. Below we build one with plain matplotlib. To produce a publication-quality version, use plotting packages such as corner which can generate corner plots with a single function call.
labels = ['area', 'gamma', 'omega_0', 'sigma']
cols = np.column_stack([column_for(p) for p in (area, gamma, omega_0, sigma)])
n = len(labels)
fig, axes = plt.subplots(n, n, figsize=(8, 8))
for i in range(n):
for j in range(n):
ax = axes[i, j]
if j > i:
ax.set_visible(False)
continue
if i == j:
ax.hist(cols[:, i], bins=40, color='C0', histtype='stepfilled', alpha=0.7)
ax.set_yticks([])
else:
ax.hexbin(cols[:, j], cols[:, i], gridsize=30, cmap='Blues', mincnt=1)
if i == n - 1:
ax.set_xlabel(labels[j])
else:
ax.set_xticklabels([])
if j == 0:
ax.set_ylabel(labels[i])
else:
ax.set_yticklabels([])
fig.tight_layout()
plt.show()
Posterior-predictive band¶
With the full posterior in hand we can propagate uncertainty through the model without assuming Gaussianity. Unlike error propagation from an MLE fit — which assumes parameters are normally distributed and independent — the posterior draws naturally capture any skew, heavy tails, and correlations between parameters. Below we pick a few hundred random draws, evaluate the model at each, and plot the resulting 95% credible band alongside the data.
rng = np.random.default_rng(seed=0)
n_draws = 300
indices = rng.choice(draws.shape[0], size=n_draws, replace=False)
predictions = np.empty((n_draws, omega.size))
saved = {p.unique_name: p.value for p in (area, gamma, omega_0, sigma)}
try:
for k, idx in enumerate(indices):
area.value = draws[idx, name_to_col[area.unique_name]]
gamma.value = draws[idx, name_to_col[gamma.unique_name]]
omega_0.value = draws[idx, name_to_col[omega_0.unique_name]]
sigma.value = draws[idx, name_to_col[sigma.unique_name]]
predictions[k] = intensity_model(omega)
finally:
for p in (area, gamma, omega_0, sigma):
p.value = saved[p.unique_name]
lo = np.percentile(predictions, 2.5, axis=0)
hi = np.percentile(predictions, 97.5, axis=0)
mid = np.percentile(predictions, 50, axis=0)
fig, ax = plt.subplots()
ax.errorbar(omega, intensity_obs, intensity_error, fmt='.', label='data')
ax.fill_between(omega, lo, hi, color='C1', alpha=0.3, label='95% posterior band')
ax.plot(omega, mid, '-', color='C1', label='posterior median')
ax.set(xlabel='$\\omega$/meV', ylabel='$I(\\omega)$')
ax.legend()
plt.show()