Add Bayesian posterior sampling to Analysis and ParameterAnalysis

docs/docs/tutorials/bayesian.ipynb

@@ -27,7 +27,8 @@
     "import easydynamics.sample_model as sm\n",
     "from easydynamics.analysis.analysis1d import Analysis1d\n",
     "\n",
-    "%matplotlib inline"
+    "# Make the plots interactive; the Q sliders need the widget backend\n",
+    "%matplotlib widget"
    ]
   },
   {
@@ -152,7 +153,9 @@
     "- `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."
+    "Sampling never moves your parameters: their values are restored afterwards, so the model is left exactly as the fit left it.\n",
+    "\n",
+    "For a long run, `progress=True` shows a single self-updating line with the percentage of generations completed, closed with `Sampling: done`. The percentage is based on the backend's own estimate of the run length, which can be too high, so a finished run may close the line before reaching 100%."
    ]
   },
   {
@@ -162,7 +165,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "results = analysis.bayesian.sample(samples=4000, burn=300, thin=2)\n",
+    "results = analysis.bayesian.sample(samples=4000, burn=300, thin=2, progress=True)\n",
     "\n",
     "print(f'Collected {results.draws.shape[0]} draws for {results.draws.shape[1]} parameters.')"
    ]
@@ -227,14 +230,56 @@
     "analysis.bayesian.plot_corner()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "e52ab7b4",
+   "metadata": {},
+   "source": [
+    "### One parameter at a time\n",
+    "\n",
+    "`plot_marginal()` pulls a single parameter's posterior out of the chain: a histogram of its draws, with the median and the 16/84 percentiles — the same numbers `summary()` reports — marked on it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2cdf2451",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.plot_marginal('Res. Gauss width')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "683943ef",
+   "metadata": {},
+   "source": [
+    "### The correlation matrix at a glance\n",
+    "\n",
+    "Where the corner plot shows every pairwise distribution, `plot_correlations()` reduces each panel to a single number — the Pearson correlation between the two parameters — and colour-codes the grid. It is the quickest way to spot which parameters the data cannot tell apart."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d577ef22",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "analysis.bayesian.plot_correlations()"
+   ]
+  },
   {
    "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."
+    "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.\n",
+    "\n",
+    "The band defaults to the 68% credible interval; `credible_interval=95.0` widens it to 95%."
    ]
   },
   {
@@ -276,6 +321,106 @@
     "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": "ef0998e4",
+   "metadata": {},
+   "source": [
+    "## Several Q values at once\n",
+    "\n",
+    "Everything so far used `Analysis1d`, a single Q slice. A full `Analysis` can sample too, either way round:\n",
+    "\n",
+    "- `fit_method='independent'` gives each Q its own chain. Cheaper, and the Q values cannot influence one another.\n",
+    "- `fit_method='simultaneous'` runs a single chain over every Q at once, which is what you need when parameters are shared across Q. It costs considerably more, because DREAM runs a number of chains proportional to the parameter count and a simultaneous run has every Q's parameters in play together.\n",
+    "\n",
+    "Sampling is much slower than fitting, so it is worth trying a few Q values before committing to all of them. Passing `Q_index` samples just that one."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b9d3b565",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Fresh models, so this analysis is independent of the single-Q one above rather than\n",
+    "# sharing its already-sampled components.\n",
+    "all_q_components = sm.ComponentCollection()\n",
+    "all_q_components.append_component(sm.Gaussian(width=0.1, area=1, name='Res. Gauss'))\n",
+    "\n",
+    "full_analysis = edyn.Analysis(\n",
+    "    display_name='Vanadium, all Q',\n",
+    "    experiment=vanadium_experiment,\n",
+    "    sample_model=sm.SampleModel(components=all_q_components),\n",
+    "    instrument_model=sm.InstrumentModel(\n",
+    "        background_model=sm.BackgroundModel(components=sm.Polynomial(coefficients=[0.001])),\n",
+    "    ),\n",
+    ")\n",
+    "full_analysis.fit(fit_method='independent')\n",
+    "\n",
+    "for Q_index in (4, 8, 12):\n",
+    "    full_analysis.analysis_list[Q_index].bayesian.suggest_bounds().apply()\n",
+    "    full_analysis.bayesian.sample(\n",
+    "        fit_method='independent', Q_index=Q_index, samples=3000, burn=200, thin=2\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "740fa625",
+   "metadata": {},
+   "source": [
+    "`bayesian.summary()` gathers the per-Q chains into one table, labelled by Q index. Each row is a marginal distribution, and a marginal is well defined within its own chain, so collecting them says nothing that was not sampled."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "511ef922",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "full_analysis.bayesian.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "12347890",
+   "metadata": {},
+   "source": [
+    "Corner plots are the one thing that cannot be gathered up. The chains were run separately, so no draw pairs a parameter at one Q with a parameter at another, and a combined figure would show correlations that came from how the sampling was run rather than from the data.\n",
+    "\n",
+    "So `plot_corner()` steps through them instead. The slider offers only the Q values that were actually sampled — 4, 8 and 12 here — and `plot_corner(Q_index=8)` goes straight to one of them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "38d0b23c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "full_analysis.bayesian.plot_corner()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2c16b5e4",
+   "metadata": {},
+   "source": [
+    "The other plots work the same way over independent chains: `plot_posterior_predictive()`, `plot_trace()`, `plot_marginal()` and `plot_correlations()` all show a Q slider in a notebook — the predictive plot through the same slider machinery as `plot_data_and_model()` — take `Q_index=` to go straight to one Q, and outside a notebook name the sampled Q indices instead."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "86b57835",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "full_analysis.bayesian.plot_posterior_predictive(n_draws=100)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "a3448aee",
@@ -287,7 +432,13 @@
     "\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."
+    "**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.\n",
+    "\n",
+    "**Several Q values at once.** An `Analysis` can sample either way. `fit_method='independent'` gives each Q its own chain, which is cheaper; `fit_method='simultaneous'` runs one chain over every Q, which is what you need when parameters are shared across Q. `bayesian.summary()` gathers the per-Q chains into one table either way.\n",
+    "\n",
+    "Corner plots are the exception. Independent chains share no draws, so nothing pairs a parameter at one Q with a parameter at another, and combining them would show correlations that came from how the sampling was run rather than from the data. `analysis.bayesian.plot_corner()` therefore shows one Q at a time: pass `Q_index`, or leave it out in a notebook to get a slider across the sampled Q values.\n",
+    "\n",
+    "**Reproducibility.** Two identical `sample()` calls will not give identical chains: the DREAM backend draws from global random state and exposes no seed. Judge results by whether the summary is stable when the chain is extended, not by exact repetition."
    ]
   }
  ],

docs/docs/tutorials/tutorial1_brownian.ipynb

@@ -531,6 +531,71 @@
     "parameter_analysis.get_all_parameters()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "163c27bb",
+   "metadata": {},
+   "source": [
+    "### How certain are the diffusion parameters?\n",
+    "\n",
+    "The uncertainties printed above come from the curvature of $\\chi^2$ at the best fit. That is a good estimate when the parameters are uncorrelated and their uncertainties are roughly Gaussian, but $D$ and the scale are fitted to the same curve and need not be either. A Bayesian analysis maps the full posterior instead, so we can check.\n",
+    "\n",
+    "The bounds act as the prior, so every free parameter needs finite ones first. `bayesian.suggest_bounds()` proposes them from the fit and is advisory until `.apply()` is called."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "de797763",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": [],
+   "source": [
+    "suggestions = parameter_analysis.bayesian.suggest_bounds()\n",
+    "print(suggestions)\n",
+    "suggestions.apply()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "604cd4e9",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": [],
+   "source": [
+    "parameter_analysis.bayesian.sample(samples=4000, burn=200, thin=2)\n",
+    "parameter_analysis.bayesian.summary()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "2e0cdab0",
+   "metadata": {},
+   "source": [
+    "The corner plot shows how the two parameters trade off against each other. A tilted, narrow ridge means the data pins down a combination of $D$ and the scale more tightly than either one separately."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "id": "6208979b",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": [],
+   "source": [
+    "parameter_analysis.bayesian.plot_corner()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "154d6137",
+   "metadata": {},
+   "source": [
+    "Notice that these credible intervals are **much narrower** than the uncertainties printed further up, and that difference is worth understanding rather than trusting.\n",
+    "\n",
+    "The least-squares fit of the widths has a reduced $\\chi^2$ of about 150: the Brownian model does not describe the fitted widths to within their error bars. `lmfit` responds by inflating its reported uncertainties by the square root of that, roughly a factor of 12, on the assumption that a poor fit means the input uncertainties were understated. The sampler makes no such adjustment — it takes the stated uncertainties at face value — so its intervals come out around twelve times tighter.\n",
+    "\n",
+    "Neither is simply right. The gap is a signal that the two-step model is not capturing the data, which is exactly what we address next by fitting the diffusion model to all the data at once."
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "fc2f8434",

src/easydynamics/analysis/__init__.py

@@ -8,12 +8,14 @@
 from easydynamics.analysis.posterior import ParameterPosterior
 from easydynamics.analysis.posterior import PosteriorSummary
 from easydynamics.analysis.posterior_labels import ParameterLabels
+from easydynamics.analysis.posterior_sampling import MultiQPosteriorSampler
 from easydynamics.analysis.posterior_sampling import PosteriorSampler
 
 __all__ = [
     'Analysis',
     'BoundsSuggestion',
     'BoundsSuggestions',
+    'MultiQPosteriorSampler',
     'ParameterAnalysis',
     'ParameterLabels',
     'ParameterPosterior',