docs: chapter 1 tutorial 4/6 review feedback

.gitignore

@@ -4,6 +4,10 @@ root.log
 __pycache__/
 *.pyc
 **/images/
+# Exception: tutorial 4 embeds these six committed figures (bad/okay/good fit and
+# normalized residual maps) via raw GitHub URLs; without this they are silently
+# skipped by `git add` and the tutorial renders six broken images.
+!scripts/chapter_1_introduction/images/
 
 output/
 dataset/

llms-full.txt

@@ -23,7 +23,7 @@ AUTO-GENERATED by PyAutoHands — do not edit by hand; regenerate with generate.
 - [Tutorial 6: Gradients](scripts/chapter_1_introduction/tutorial_6_gradients.py): In tutorial 3, describing how a maximum likelihood estimator works, we said it "evaluates the likelihood at nearby points to estimate the gradient, determining the direction in which to move up in parameter space". That sentence slipped by quickly, but it hides one of the most important ideas in model-fitting: the search does not have to be blind, because the likelihood surface has a slope, and the slope tells us which way is up.
   - Contents: Data, Model, Analysis, Gradients, Finite Differencing, JAX and Autodiff, Maximum Likelihood Estimation (MLE), Markov Chain Monte Carlo (MCMC), Nested Sampling, Errors From Curvature, Wrap Up
 - [Tutorial 7: The Details](scripts/chapter_1_introduction/tutorial_7_the_details.py): In tutorial 6 we opened up the non-linear search and looked at how it moves. We saw that some searches only ever ask the likelihood function "what is the value here?", whilst gradient based searches also ask "and which way is uphill?", and that **PyAutoFit** can answer the second question exactly by using JAX to differentiate the likelihood function.
-  - Contents: Data, Model, Analysis, Parameterization, Assertions, Plateaus, Comparing Searches, Clipping, NaN Diagnostics, Unit Cube Vs Physical, Summary
+  - Contents: Data, Model, Analysis, Parameterization, Assertions, Plateaus, Comparing Searches, Clipping, NaN Diagnostics, Hamiltonian Diagnostics, Unit Cube Vs Physical, Summary
 - [Tutorial 8: Scientific Workflow](scripts/chapter_1_introduction/tutorial_8_scientific_workflow.py): You can now compose a model, fit it to data and interpret its results. A scientific study often repeats these steps for many datasets, competing models and different non-linear searches. The next challenge is keeping those fits organized so that you can understand and compare them.
 - [Tutorial Optional: Bayesian Formalism](scripts/chapter_1_introduction/tutorial_optional_bayesian_formalism.py): Every tutorial in this chapter has been an exercise in Bayesian inference, and not one of them wrote down a single equation of probability theory. That was deliberate. The **HowToFit** lectures are built on the conviction that you can learn to fit models to data properly, and interpret the results correctly, without first sitting through a course on Bayesian statistics.
   - Contents: Bayes Theorem, Data, The Model, The Likelihood, The Prior, The Posterior, Maximum Likelihood Vs Maximum A Posteriori, Marginalization, The Evidence, Gradients, Wrap Up

notebooks/chapter_1_introduction/tutorial_4_why_modeling_is_hard.ipynb

@@ -263,8 +263,8 @@
     "To define the Analysis class for this model-fit, we need to ensure that the `log_likelihood_function` can handle an \n",
     "instance containing multiple 1D profiles. Below is an expanded explanation and the corresponding class definition:\n",
     "\n",
-    "The log_likelihood_function will now assume that the instance it receives consists of multiple Gaussian profiles. \n",
-    "For each Gaussian in the instance, it will compute the model_data and then sum these to create the overall `model_data` \n",
+    "The `log_likelihood_function` will now assume that the instance it receives consists of multiple Gaussian profiles. \n",
+    "For each Gaussian in the instance, it will compute the `model_data` and then sum these to create the overall `model_data` \n",
     "that is compared to the observed data."
    ]
   },
@@ -573,7 +573,7 @@
     "plt.plot(range(data.shape[0]), model_data, color=\"r\")\n",
     "for model_data_1d_individual in model_data_list:\n",
     "    plt.plot(range(data.shape[0]), model_data_1d_individual, \"--\")\n",
-    "plt.title(f\"Fit (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Fit (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(\"Profile normalization\")\n",
     "plt.show()\n",
@@ -613,7 +613,7 @@
     "    capsize=2,\n",
     "    linestyle=\"\",\n",
     ")\n",
-    "plt.title(f\"Residuals (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Residuals (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(\"Residuals\")\n",
     "plt.show()\n",
@@ -645,7 +645,7 @@
     "residual_map = data - model_data\n",
     "normalized_residual_map = residual_map / noise_map\n",
     "plt.plot(xvalues, normalized_residual_map, color=\"k\")\n",
-    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(r\"Normalized Residuals ($\\sigma$)\")\n",
     "plt.show()\n",
@@ -863,7 +863,7 @@
     "plt.plot(range(data.shape[0]), model_data, color=\"r\")\n",
     "for model_data_1d_individual in model_data_list:\n",
     "    plt.plot(range(data.shape[0]), model_data_1d_individual, \"--\")\n",
-    "plt.title(f\"Fit (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Fit (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(\"Profile normalization\")\n",
     "plt.show()\n",
@@ -873,7 +873,7 @@
     "residual_map = data - model_data\n",
     "normalized_residual_map = residual_map / noise_map\n",
     "plt.plot(xvalues, normalized_residual_map, color=\"k\")\n",
-    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(r\"Normalized Residuals ($\\sigma$)\")\n",
     "plt.show()\n",
@@ -1035,7 +1035,7 @@
     "plt.plot(range(data.shape[0]), model_data, color=\"r\")\n",
     "for model_data_1d_individual in model_data_list:\n",
     "    plt.plot(range(data.shape[0]), model_data_1d_individual, \"--\")\n",
-    "plt.title(f\"Fit (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Fit (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(\"Profile normalization\")\n",
     "plt.show()\n",
@@ -1045,7 +1045,7 @@
     "residual_map = data - model_data\n",
     "normalized_residual_map = residual_map / noise_map\n",
     "plt.plot(xvalues, normalized_residual_map, color=\"k\")\n",
-    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(r\"Normalized Residuals ($\\sigma$)\")\n",
     "plt.show()\n",
@@ -1215,7 +1215,7 @@
     "plt.plot(range(data.shape[0]), model_data, color=\"r\")\n",
     "for model_data_1d_individual in model_data_list:\n",
     "    plt.plot(range(data.shape[0]), model_data_1d_individual, \"--\")\n",
-    "plt.title(f\"Fit (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Fit (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(\"Profile normalization\")\n",
     "plt.show()\n",
@@ -1225,7 +1225,7 @@
     "residual_map = data - model_data\n",
     "normalized_residual_map = residual_map / noise_map\n",
     "plt.plot(xvalues, normalized_residual_map, color=\"k\")\n",
-    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood})\")\n",
+    "plt.title(f\"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})\")\n",
     "plt.xlabel(\"x values of profile\")\n",
     "plt.ylabel(r\"Normalized Residuals ($\\sigma$)\")\n",
     "plt.show()\n",

notebooks/chapter_1_introduction/tutorial_6_gradients.ipynb

@@ -803,7 +803,7 @@
     "advanced together in one compiled call using `vmap`, so twelve lanes cost far less than twelve times one lane.\n",
     "\n",
     "This search is given a `name` and `path_prefix`, so its results are written to the `output` folder as tutorial 5\n",
-    "did, because there is a file in there we are about to read."
+    "did, where you can inspect them once the fit completes."
    ]
   },
   {
@@ -877,40 +877,6 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Because this search wrote its results to hard disk, its output folder contains a file called `search.summary`: how\n",
-    "long the search took, how long one log likelihood evaluation took and, for gradient searches, diagnostics on how\n",
-    "often things went wrong. We read it back below."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "metadata": {},
-   "source": [
-    "search_summary_path = path.join(str(search.paths.output_path), \"search.summary\")\n",
-    "\n",
-    "if path.exists(search_summary_path):\n",
-    "    with open(search_summary_path) as f:\n",
-    "        print(f.read())"
-   ],
-   "outputs": [],
-   "execution_count": null
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "The block at the bottom, headed `Resampling Info`, contains two entries only a gradient search can report:\n",
-    "\n",
-    "- `Value-NaN Lane-Steps`: the number of times a lane stepped somewhere the log likelihood could not be computed at\n",
-    "  all, most often because it stepped outside the priors.\n",
-    "\n",
-    "- `Gradient-NaN Lane-Steps`: the number of times the log likelihood *was* computable but its gradient was not. This\n",
-    "  is the sneakier of the two, because such a lane does not crash or die, it simply stops moving while continuing to\n",
-    "  look perfectly healthy.\n",
-    "\n",
-    "Both counters are usually small and harmless, but they are the vocabulary you need to diagnose a gradient fit that\n",
-    "has gone quietly wrong. Tutorial 7 explains where they come from and what to do about them.\n",
-    "\n",
     "__Markov Chain Monte Carlo (MCMC)__\n",
     "\n",
     "In tutorial 3 we used `Emcee`, whose walkers propose a step, compute the likelihood there and accept or reject the\n",
@@ -951,18 +917,19 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now the gradient-aware alternative, `BlackJAXNUTS`, which is Hamiltonian Monte Carlo. The physical picture behind it\n",
-    "is genuinely helpful.\n",
-    "\n",
-    "Imagine the likelihood surface turned upside down, so its peak becomes a valley, and place a ball on the resulting\n",
-    "landscape. Give it a random flick and let it roll: it accelerates down slopes, coasts up the other side and travels\n",
-    "a long way while staying in regions the landscape favours. That trajectory is computed from the gradient at each\n",
-    "moment, which is what autodiff hands us for free. Where the ball stops becomes the next sample, and because it\n",
-    "travelled a long, informed distance rather than a small random hop, consecutive samples are far less similar.\n",
-    "\n",
-    "\"NUTS\" stands for the No U-Turn Sampler, which solves the awkward choice here: how long to let the ball roll. Roll\n",
-    "too briefly and you wasted the gradient; roll too long and the ball curves back on itself. NUTS stops the trajectory\n",
-    "when it starts doubling back. Being a gradient method, it needs the JAX analysis."
+    "Now the gradient-aware alternative, `BlackJAXNUTS`, which is Hamiltonian Monte Carlo. It is still MCMC, with walkers\n",
+    "moving through parameter space and proposals accepted or rejected, but the walker now knows which way to step,\n",
+    "because autodiff hands it the gradient at every point it visits.\n",
+    "\n",
+    "That changes how far one proposal can usefully travel. Instead of a single small random hop, the walker follows the\n",
+    "gradient along a trajectory of many small steps, staying in the regions the likelihood favours the whole way, and\n",
+    "ends up somewhere genuinely far from where it started. Consecutive samples are therefore far less correlated than\n",
+    "`Emcee`'s.\n",
+    "\n",
+    "The one thing left to choose is how long to follow that trajectory. Stop too early and the gradient was wasted;\n",
+    "carry on too long and the path curves back on itself and returns to where it began. \"NUTS\" stands for the No U-Turn\n",
+    "Sampler, which watches for that doubling back and ends the trajectory there. Being a gradient method, it needs the\n",
+    "JAX analysis."
    ]
   },
   {
@@ -993,40 +960,6 @@
    "outputs": [],
    "execution_count": null
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Hamiltonian sampling comes with its own diagnostics, stored in the `samples_info` dictionary. Three are worth\n",
-    "knowing:\n",
-    "\n",
-    "- `n_divergent`: the number of trajectories which \"diverged\", meaning the ball flew off to infinity instead of\n",
-    "  following the landscape. A handful is tolerable; many means the steps are too large and the samples cannot be\n",
-    "  trusted.\n",
-    "\n",
-    "- `ess_min`: the \"effective sample size\" of the worst constrained parameter. Consecutive samples are correlated, so\n",
-    "  300 samples are worth fewer than 300 independent draws, and this says how many they are worth.\n",
-    "\n",
-    "- `mean_acceptance`: the fraction of proposed trajectories accepted, which for NUTS should sit high, around the 0.8\n",
-    "  the warm up phase tunes towards. A low value means the sampler is struggling."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "metadata": {},
-   "source": [
-    "samples = result.samples\n",
-    "\n",
-    "print(\"Diagnostics of the Hamiltonian Monte Carlo fit:\\n\")\n",
-    "print(f\"Number of divergent trajectories = {samples.samples_info.get('n_divergent')}\")\n",
-    "print(f\"Minimum effective sample size    = {samples.samples_info.get('ess_min')}\")\n",
-    "print(\n",
-    "    f\"Mean acceptance rate             = {samples.samples_info.get('mean_acceptance')}\"\n",
-    ")"
-   ],
-   "outputs": [],
-   "execution_count": null
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1229,7 +1162,7 @@
    "source": [
     "__Wrap Up__\n",
     "\n",
-    "This tutorial took the one sentence tutorial 3 used to describe how an MLE search moves and unpacked it:\n",
+    "Tutorial 3 described how an MLE search moves in a single sentence. This tutorial unpacked that sentence:\n",
     "\n",
     "1. **Gradients**: the gradient of the log likelihood is a vector with one entry per free parameter, evaluated at a\n",
     "point, pointing in the direction the likelihood increases fastest.\n",
@@ -1244,8 +1177,8 @@
     "compiled call, making it far harder to trap in a local maximum.\n",
     "\n",
     "4. **MCMC**: `Emcee` proposes random steps and accepts or rejects them, never asking which way is up.\n",
-    "`BlackJAXNUTS` rolls a ball across the landscape using the gradient to shape its trajectory, producing far less\n",
-    "correlated samples for the same number of steps.\n",
+    "`BlackJAXNUTS` follows the gradient along a trajectory before taking each sample, producing far less correlated\n",
+    "samples for the same number of steps.\n",
     "\n",
     "5. **Nested sampling**: `Nautilus` replaces low likelihood live points with higher likelihood ones drawn from the\n",
     "priors, a procedure with no direction of travel for a gradient to inform. JAX speeds it up by evaluating batches of\n",