Add Bayesian posterior sampling to Analysis1d

docs/docs/tutorials/bayesian.ipynb

@@ -0,0 +1,306 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "eac0b8bb",
+   "metadata": {},
+   "source": [
+    "# Bayesian analysis\n",
+    "\n",
+    "Fitting with `fit()` finds the single set of parameter values that best matches the data, and reports an uncertainty derived from the curvature of $\\chi^2$ at that point. That uncertainty is only trustworthy when the parameters are uncorrelated and their uncertainties are close to Gaussian, which in QENS is often not the case.\n",
+    "\n",
+    "A **Bayesian** analysis answers a different question: instead of one best point, it maps out the whole *posterior distribution* over the parameters. From that you can read off credible intervals that stay honest when parameters are correlated or their distributions are skewed, and you can see the correlations directly.\n",
+    "\n",
+    "EasyDynamics does this with the DREAM sampler from [BUMPS](https://bumps.readthedocs.io/), through `bayesian.sample()`."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "02cb7aec",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pooch\n",
+    "\n",
+    "import easydynamics as edyn\n",
+    "import easydynamics.sample_model as sm\n",
+    "from easydynamics.analysis.analysis1d import Analysis1d\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0499fea7",
+   "metadata": {},
+   "source": [
+    "## Load the data\n",
+    "\n",
+    "We use the same artificial vanadium measurement as the [Analysis 1D](analysis1d.ipynb) tutorial, and analyse a single Q slice."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fb407621",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vanadium_experiment = edyn.Experiment('Vanadium')\n",
+    "\n",
+    "file_path = pooch.retrieve(\n",
+    "    url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/master/docs/docs/tutorials/data/vanadium_data_example.h5',\n",
+    "    known_hash='16cc1b327c303feeb88fb9dda5390dc4880b62396b1793f98c6fef0b27c7b873',\n",
+    ")\n",
+    "\n",
+    "vanadium_experiment.load_hdf5(filename=file_path)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fcdfd395",
+   "metadata": {},
+   "source": [
+    "## Build the model and fit it\n",
+    "\n",
+    "As in [Tutorial 1](tutorial1_brownian.ipynb), a vanadium measurement is modelled with the Gaussian as the *sample*: what is being measured is the resolution function itself, so there is nothing to convolve it with.\n",
+    "\n",
+    "Sampling does not require a fit first, but it benefits from one: DREAM starts its chains in a small ball around the parameters' current values, so beginning from fitted values means less burn-in is needed before the chains reach the interesting region."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "de3297cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "vanadium_components = sm.ComponentCollection()\n",
+    "vanadium_components.append_component(sm.Gaussian(width=0.1, area=1, name='Res. Gauss'))\n",
+    "\n",
+    "instrument_model = sm.InstrumentModel(\n",
+    "    background_model=sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001])),\n",
+    ")\n",
+    "\n",
+    "analysis = Analysis1d(\n",
+    "    display_name='Vanadium Analysis',\n",
+    "    experiment=vanadium_experiment,\n",
+    "    sample_model=sm.SampleModel(components=vanadium_components),\n",
+    "    instrument_model=instrument_model,\n",
+    "    Q_index=5,\n",
+    ")\n",
+    "\n",
+    "fit_result = analysis.fit()\n",
+    "print(f'reduced chi-squared = {fit_result.reduced_chi2:.4f}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b0a709f7",
+   "metadata": {},
+   "source": [
+    "## Bounds are the prior\n",
+    "\n",
+    "In DREAM, each parameter's `min` and `max` define a uniform prior, so **every free parameter must have finite bounds** before sampling. Most parameters start with at least one infinite bound, so `bayesian.sample()` would refuse to run.\n",
+    "\n",
+    "`bayesian.suggest_bounds()` proposes bounds from the fitted values and uncertainties. It is advisory: it changes nothing until you call `.apply()`, and it only ever fills in an *infinite* bound, so physical limits you have already set (an area that cannot go below zero, say) are left alone."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a004cd83",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "suggestions = analysis.bayesian.suggest_bounds()\n",
+    "print(suggestions)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ab263f41",
+   "metadata": {},
+   "source": [
+    "The defaults are deliberately generous — 10 standard deviations plus 20% of the value. Because the bounds are a uniform prior, being too *narrow* is the dangerous mistake: it truncates the posterior and makes the uncertainty look smaller than it is. The 20% term is there for parameters whose fitted uncertainty comes back as zero. All three settings (`n_sigma`, `relative_pad`, `absolute_floor`) can be adjusted, and you can always set `min` and `max` by hand.\n",
+    "\n",
+    "It is worth reading the table before applying it. A suggestion many orders of magnitude larger than the parameter itself is a useful warning sign: it means the fit returned a huge uncertainty, which usually happens because two parameters are **degenerate** — the data determines only some combination of them, so one can grow while the other shrinks with no effect on the fit. That is a problem to fix in the model, not with the sampler."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1a92a899",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "changed = suggestions.apply()\n",
+    "print(f'Applied bounds to: {[parameter.name for parameter in changed]}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8cf71e53",
+   "metadata": {},
+   "source": [
+    "## Sample the posterior\n",
+    "\n",
+    "`bayesian.sample()` runs the chains. The three numbers that matter are:\n",
+    "\n",
+    "- `samples` — how many draws to collect in total. More is better, at linear cost.\n",
+    "- `burn` — generations discarded at the start, while the chains are still travelling towards the bulk of the posterior.\n",
+    "- `thin` — keep only every n-th generation, which reduces the correlation between neighbouring draws.\n",
+    "\n",
+    "Sampling never moves your parameters: their values are restored afterwards, so the model is left exactly as the fit left it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "861fa264",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "results = analysis.bayesian.sample(samples=4000, burn=300, thin=2)\n",
+    "\n",
+    "print(f'Collected {results.draws.shape[0]} draws for {results.draws.shape[1]} parameters.')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f2ce5232",
+   "metadata": {},
+   "source": [
+    "## Did the chains converge?\n",
+    "\n",
+    "Always look at the traces before trusting the numbers. A converged chain looks like a \"hairy caterpillar\": noisy, but flat and stationary. A visible drift or slow wander means the chain has not settled and needs a longer burn-in or more samples."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b2fe638c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.plot_trace()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2b3ea7b4",
+   "metadata": {},
+   "source": [
+    "## Summarize the posterior\n",
+    "\n",
+    "`bayesian.summary()` reports the median and the 68% credible interval of each parameter, under the parameter's own name and unit. The interval is asymmetric in general, which is precisely the information a single symmetric error bar throws away."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8470375e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "60b4a448",
+   "metadata": {},
+   "source": [
+    "## Correlations between parameters\n",
+    "\n",
+    "The corner plot is the part least available from a least-squares fit. The diagonal shows each parameter's own distribution; each off-diagonal panel shows a pair. A round blob means the two are independent, while a tilted, narrow ridge means they are correlated and the data constrains only a combination of them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c96509cc",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.plot_corner()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8f819c75",
+   "metadata": {},
+   "source": [
+    "## Does the model actually describe the data?\n",
+    "\n",
+    "The posterior predictive plot re-evaluates the model for a sample of posterior draws and shades the region they cover. If the data wanders outside the band in a systematic way, the model is missing a feature, and no amount of parameter tuning will fix it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "cce95199",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.plot_posterior_predictive(n_draws=100)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5947da0f",
+   "metadata": {},
+   "source": [
+    "## Continuing and storing a chain\n",
+    "\n",
+    "If the traces suggest the chain needs to run longer, `extend_sampling()` continues the existing chain rather than starting over, so nothing already computed is thrown away."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "974e486b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "extended = analysis.bayesian.extend(additional_samples=1000, thin=2)\n",
+    "print(f'Chain now holds {extended.draws.shape[0]} draws.')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "69793d31",
+   "metadata": {},
+   "source": [
+    "Chains are expensive, so they can be saved and reloaded with `analysis.bayesian.save(path)` and `analysis.bayesian.load(path)`. A reloaded chain can be summarized, plotted, or extended further, exactly like a fresh one."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a3448aee",
+   "metadata": {},
+   "source": [
+    "## Things to watch out for\n",
+    "\n",
+    "**Data without uncertainties.** If your data carries no variances, the weights fall back to 1, which means the sampler assumes a noise level of 1 in whatever units the intensity happens to be. Least-squares does not care, since that scale cancels out of the best-fit position, but a posterior *does*: its width scales directly with the assumed noise, so the credible intervals will be wrong by whatever factor the true noise differs from 1. Bayesian analysis is not a way to avoid needing uncertainties on your data.\n",
+    "\n",
+    "**Sampling only some parameters.** `bayesian.sample(parameters=[...])` restricts the chain to a subset, which is faster because the number of chains scales with the number of parameters. Be careful with the result: the other parameters are held *fixed*, which is not the same as averaging over them. The intervals you get are conditional on those fixed values, and will be too narrow whenever the parameters are correlated.\n",
+    "\n",
+    "**Degenerate parameters.** As seen above, they are a modelling problem rather than a sampling one. `bayesian.suggest_bounds()` returning absurd values, or a warning that the posterior has piled up against its bounds, are both signs to go back and look at the model."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "default",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

docs/docs/tutorials/index.md

@@ -62,3 +62,6 @@ tutorials.
 - [Analysis](analysis.ipynb) - Learn how to fit a model to your data.
 - [Analysis 1D](analysis1d.ipynb) - Learn how to fit a model to your
   data at a particular Q.
+- [Bayesian analysis](bayesian.ipynb) - Learn how to map out the full
+  posterior distribution of your parameters, including their
+  correlations, instead of a single best-fit point.

docs/mkdocs.yml

@@ -217,6 +217,7 @@ nav:
           - Experiment: tutorials/experiment.ipynb
           - Analysis: tutorials/analysis.ipynb
           - Analysis 1D: tutorials/analysis1d.ipynb
+          - Bayesian analysis: tutorials/bayesian.ipynb
   - API Reference:
       - API Reference: api-reference/index.md
       - analysis: api-reference/analysis.md

pixi.toml

@@ -100,7 +100,13 @@ user = { features = ['py-max', 'user'] }
 unit-tests = 'python -m pytest tests/unit/ --color=yes -v'
 functional-tests = 'python -m pytest tests/functional/ --color=yes -v'
 integration-tests = 'python -m pytest tests/integration/ --color=yes -n auto -v'
-notebook-tests = 'python -m pytest --nbmake docs/docs/tutorials/ --nbmake-timeout=1200 --color=yes -n auto -v'
+# Warm the pooch cache first. Several notebooks fetch the same file, and running them with
+# '-n auto' has the workers race: one writes the file while another opens it, which fails on
+# Windows. Fetching up front leaves the parallel run with nothing to do but read.
+prefetch-tutorial-data = 'python tools/prefetch_tutorial_data.py'
+notebook-tests = { cmd = 'python -m pytest --nbmake docs/docs/tutorials/ --nbmake-timeout=1200 --color=yes -n auto -v', depends-on = [
+  'prefetch-tutorial-data',
+] }
 
 test = { depends-on = ['unit-tests'] }
 

pyproject.toml

@@ -23,17 +23,18 @@ classifiers = [
 ]
 requires-python = '>=3.12'
 dependencies = [
-  'easyscience', # The base library of the EasyScience framework
-  'pooch',       # Data downloader
-  'darkdetect',  # Detecting dark mode (system-level)
-  'plopp',       # Plotting library
-  'jupyterlab',  # Jupyter notebooks
-  'pixi-kernel', # Pixi Jupyter kernel
-  'ipykernel',   # Jupyter kernel (required for running notebooks)
-  'ipywidgets',  # Widgets (needed for interactive matplotlib backends)
-  'ipympl',      # Matplotlib Jupyter widget backend (%matplotlib widget)
-  'IPython',     # Interactive Python shell
-  'sympy',       # Symbolic mathematics (used for expression components)
+  'easyscience>=2.5.1', # The base library of the EasyScience framework. 2.5.1 adds fitting.Sampler
+  'matplotlib',         # Plotting (posterior trace, corner, and predictive plots)
+  'pooch',              # Data downloader
+  'darkdetect',         # Detecting dark mode (system-level)
+  'plopp',              # Plotting library
+  'jupyterlab',         # Jupyter notebooks
+  'pixi-kernel',        # Pixi Jupyter kernel
+  'ipykernel',          # Jupyter kernel (required for running notebooks)
+  'ipywidgets',         # Widgets (needed for interactive matplotlib backends)
+  'ipympl',             # Matplotlib Jupyter widget backend (%matplotlib widget)
+  'IPython',            # Interactive Python shell
+  'sympy',              # Symbolic mathematics (used for expression components)
 ]
 
 [project.optional-dependencies]