diff --git a/.gitignore b/.gitignore index 73edf3f..c1002ca 100644 --- a/.gitignore +++ b/.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/ diff --git a/llms-full.txt b/llms-full.txt index 343517f..9218edd 100644 --- a/llms-full.txt +++ b/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 diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard.md b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard.md index 7e19f28..da0d01e 100644 --- a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard.md +++ b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard.md @@ -56,6 +56,12 @@ import matplotlib.pyplot as plt import autofit as af ``` + .../PyAutoNerves/autonerves/workspace.py:31: UserWarning: Cannot verify the workspace at HowToFit/scripts/chapter_1_introduction is compatible with the installed library version (2026.8.17.1): no `version.minimum_library_version` or `version.workspace_version` key in config/general.yaml and no version.txt at the workspace root. + + If you cloned the workspace from `main` rather than a release tag, set `version.workspace_version_check: False` in config/general.yaml to silence this warning. The `main` branch updates more frequently than library releases, so version mismatches are expected and not actionable for `main`-branch users. + + You can also set the environment variable PYAUTO_SKIP_WORKSPACE_VERSION_CHECK=1 to disable temporarily. + warnings.warn(message) Working Directory has been set to `HowToFit` @@ -92,6 +98,19 @@ noise_map = af.util.numpy_array_from_json( ) ``` + .../PyAutoNerves/autonerves/workspace.py:31: UserWarning: The workspace at HowToFit records library version 2026.7.9.1, but the installed library is 2026.8.17.1 — more than 30 days newer. The workspace examples and configs may lag the installed API. Pull the latest workspace: + + cd HowToFit && git pull origin main + + To bypass this check, edit config/general.yaml: + + version: + workspace_version_check: False + + You can also set the environment variable PYAUTO_SKIP_WORKSPACE_VERSION_CHECK=1 to disable temporarily. + warnings.warn(message) + + Plotting the data reveals that the signal is more complex than a simple 1D Gaussian, as the wings to the left and right are more extended than what a single Gaussian profile can account for. @@ -188,8 +207,8 @@ __Analysis__ To define the Analysis class for this model-fit, we need to ensure that the `log_likelihood_function` can handle an instance containing multiple 1D profiles. Below is an expanded explanation and the corresponding class definition: -The log_likelihood_function will now assume that the instance it receives consists of multiple Gaussian profiles. -For each Gaussian in the instance, it will compute the model_data and then sum these to create the overall `model_data` +The `log_likelihood_function` will now assume that the instance it receives consists of multiple Gaussian profiles. +For each Gaussian in the instance, it will compute the `model_data` and then sum these to create the overall `model_data` that is compared to the observed data. @@ -276,7 +295,7 @@ class Analysis(af.Analysis): residual_map = self.data - model_data chi_squared_map = (residual_map / self.noise_map) ** 2.0 chi_squared = sum(chi_squared_map) - noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0)) + noise_normalization = np.sum(np.log(2 * np.pi * self.noise_map**2.0)) log_likelihood = -0.5 * (chi_squared + noise_normalization) return log_likelihood @@ -364,6 +383,27 @@ print(model.info) sigma UniformPrior [14], lower_limit = 0.0, upper_limit = 25.0 +The same model can also be visualized as a figure, making its structure easier to understand at a glance. + +The figure shows how the model is organized: which parameters belong to each component, and whether they are free, +fixed, shared, linked by an expression, solved during the fit, or not configured. `model.info` provides the +corresponding numerical details, including the prior assigned to each free parameter and the value of each fixed +parameter. + +Tutorial 3 searched 3 dimensions, whereas this model has 15 of them. Parameter space grows in volume so quickly with +that number that the same search which comfortably found the single `Gaussian` is about to struggle. + + +```python +af.ModelPlotter(model).figure() +``` + + + +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_17_0.png) + + + __Search__ We again use the nested sampling algorithm Dynesty to fit the model to the data. @@ -404,39 +444,18 @@ print("The search has finished run - you may now continue the notebook.") The non-linear search has begun running. This Jupyter notebook cell with progress once the search has completed - this could take a few minutes! - 2026-07-11 16:23:05,069 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. - - - 2026-07-11 16:23:05,079 - root - INFO - Output to hard-disk disabled, input a search name to enable. - - - 2026-07-11 16:23:05,080 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). - - - 2026-07-11 16:23:05,289 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores - - - 2026-07-11 16:23:05,324 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search - - + 2026-09-15 00:25:56,691 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. + 2026-09-15 00:25:56,692 - autofit.non_linear.search.abstract_search - INFO - On-the-fly updates of the maximum likelihood model are disabled. Set `updates: iterations_per_quick_update` in config/general.yaml to a finite number of iterations to enable them. + 2026-09-15 00:25:56,703 - root - INFO - Output to hard-disk disabled, input a search name to enable. + 2026-09-15 00:25:56,705 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). + 2026-09-15 00:25:56,956 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores + 2026-09-15 00:25:56,993 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search ~/venv/PyAuto/lib/python3.12/site-packages/dynesty/dynesty.py:194: UserWarning: Specifying slice option while using rwalk sampler does not make sense warnings.warn('Specifying slice option while using rwalk sampler' - - - 9039it [00:44, 201.70it/s, +50 | bound: 1188 | nc: 1 | ncall: 46185 | eff(%): 19.701 | loglstar: -inf < 181.115 < inf | logz: 2.376 +/- 1.690 | dlogz: 0.001 > 0.059] - - - - - 2026-07-11 16:23:51,837 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. - - - 2026-07-11 16:23:52,602 - root - INFO - Removing search internal folder. - - - 2026-07-11 16:23:52,655 - root - INFO - Search complete, returning result - - + 10227it [00:58, 176.00it/s, +50 | bound: 1347 | nc: 1 | ncall: 52117 | eff(%): 19.738 | loglstar: -inf < 117.335 < inf | logz: -85.584 +/- 1.726 | dlogz: 0.001 > 0.059] + 2026-09-15 00:26:57,999 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. + 2026-09-15 00:26:59,473 - root - INFO - Removing search internal folder. + 2026-09-15 00:26:59,545 - root - INFO - Search complete, returning result The search has finished run - you may now continue the notebook. @@ -450,8 +469,8 @@ of all 5 model components. print(result.info) ``` - Bayesian Evidence 2.37555245 - Maximum Log Likelihood 181.11496469 + Bayesian Evidence -85.58424959 + Maximum Log Likelihood 117.33513924 model Collection (N=15) gaussian_0 - gaussian_4 Gaussian (N=3) @@ -459,23 +478,23 @@ print(result.info) Maximum Log Likelihood Model: gaussian_0 - centre 49.879 + centre 50.001 ... [51 lines of output truncated] ... - centre 51.71 (41.37, 60.03) - normalization 0.00 (0.00, 0.00) - sigma 18.81 (16.40, 21.77) + centre 50.02 (50.00, 50.04) + normalization 74.01 (73.19, 74.79) + sigma 6.30 (6.26, 6.34) gaussian_2 - centre 50.00 (49.99, 50.00) - normalization 20.50 (20.32, 20.65) - sigma 1.01 (1.01, 1.02) + centre 7.43 (3.57, 15.65) + normalization 0.00 (0.00, 0.06) + sigma 1.66 (0.88, 2.43) gaussian_3 - centre 50.21 (50.11, 50.29) - normalization 125.71 (121.73, 129.59) - sigma 13.44 (13.28, 13.61) + centre 49.98 (49.94, 50.02) + normalization 204.21 (203.37, 204.96) + sigma 16.87 (16.80, 16.92) gaussian_4 - centre 49.96 (49.93, 50.00) - normalization 55.27 (54.46, 56.10) - sigma 5.58 (5.53, 5.62) + centre 82.16 (67.91, 97.41) + normalization 0.00 (0.00, 0.00) + sigma 20.89 (17.43, 23.49) instances @@ -517,7 +536,7 @@ plt.errorbar( plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -526,7 +545,7 @@ plt.close() -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_23_0.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png) @@ -554,7 +573,7 @@ plt.errorbar( capsize=2, linestyle="", ) -plt.title(f"Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Residuals") plt.show() @@ -564,7 +583,7 @@ plt.close() -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_0.png) @@ -583,23 +602,17 @@ to it being a noise fluctuation. residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") -plt.ylabel("Normalized Residuals ($\sigma$)") +plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() plt.clf() plt.close() ``` - <>:6: SyntaxWarning: invalid escape sequence '\s' - <>:6: SyntaxWarning: invalid escape sequence '\s' - /tmp/ipykernel_20726/582017931.py:6: SyntaxWarning: invalid escape sequence '\s' - plt.ylabel("Normalized Residuals ($\sigma$)") - - -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_1.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_29_0.png) @@ -745,6 +758,20 @@ print(model.info) sigma UniformPrior [37], lower_limit = 0.0, upper_limit = 25.0 +Tuning the priors changed no part of the model itself: it is still the same 5 `Gaussian`'s and the same 15 free +parameters. We have simply told the search a smaller region of parameter space to look in. + + +```python +af.ModelPlotter(model).figure() +``` + + + +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_0.png) + + + We now repeat the model-fit using these updated priors. First, you should note that the run time of the fit is significantly faster than the previous fit. This is because @@ -772,39 +799,18 @@ print("The search has finished run - you may now continue the notebook.") The non-linear search has begun running. This Jupyter notebook cell with progress once the search has completed - this could take a few minutes! - 2026-07-11 16:23:53,402 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. - - - 2026-07-11 16:23:53,412 - root - INFO - Output to hard-disk disabled, input a search name to enable. - - - 2026-07-11 16:23:53,412 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). - - - 2026-07-11 16:23:53,418 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores - - - 2026-07-11 16:23:53,453 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search - - + 2026-09-15 00:27:01,163 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. + 2026-09-15 00:27:01,165 - autofit.non_linear.search.abstract_search - INFO - On-the-fly updates of the maximum likelihood model are disabled. Set `updates: iterations_per_quick_update` in config/general.yaml to a finite number of iterations to enable them. + 2026-09-15 00:27:01,183 - root - INFO - Output to hard-disk disabled, input a search name to enable. + 2026-09-15 00:27:01,185 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). + 2026-09-15 00:27:01,197 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores + 2026-09-15 00:27:01,260 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search ~/venv/PyAuto/lib/python3.12/site-packages/dynesty/dynesty.py:194: UserWarning: Specifying slice option while using rwalk sampler does not make sense warnings.warn('Specifying slice option while using rwalk sampler' - - - 2731it [00:15, 178.52it/s, +50 | bound: 350 | nc: 1 | ncall: 14675 | eff(%): 19.015 | loglstar: -inf < 125.941 < inf | logz: 73.841 +/- 0.861 | dlogz: 0.001 > 0.059] - - - - - 2026-07-11 16:24:09,309 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. - - - 2026-07-11 16:24:09,531 - root - INFO - Removing search internal folder. - - - 2026-07-11 16:24:09,590 - root - INFO - Search complete, returning result - - + 4440it [00:45, 96.87it/s, +50 | bound: 560 | nc: 1 | ncall: 23317 | eff(%): 19.298 | loglstar: -inf < 176.353 < inf | logz: 90.086 +/- 1.157 | dlogz: 0.001 > 0.059] + 2026-09-15 00:27:48,582 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. + 2026-09-15 00:27:49,401 - root - INFO - Removing search internal folder. + 2026-09-15 00:27:49,493 - root - INFO - Search complete, returning result The search has finished run - you may now continue the notebook. @@ -827,7 +833,7 @@ plt.errorbar( plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -837,16 +843,16 @@ plt.close() residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") -plt.ylabel("Normalized Residuals ($\sigma$)") +plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() plt.clf() plt.close() ``` - Bayesian Evidence 73.84079536 - Maximum Log Likelihood 125.94072461 + Bayesian Evidence 90.08597221 + Maximum Log Likelihood 176.35326867 model Collection (N=15) gaussian_0 - gaussian_4 Gaussian (N=3) @@ -854,44 +860,38 @@ plt.close() Maximum Log Likelihood Model: gaussian_0 - sigma 6.053 + sigma 23.970 ... [51 lines of output truncated] ... - sigma 1.03 (1.03, 1.04) - centre 50.00 (49.99, 50.00) - normalization 20.93 (20.81, 21.08) + sigma 1.91 (1.13, 2.67) + centre 52.57 (52.11, 53.05) + normalization 0.13 (0.12, 0.15) gaussian_2 - sigma 10.35 (9.05, 19.85) - centre 52.24 (51.07, 52.39) - normalization 2.46 (0.46, 3.09) + sigma 14.46 (14.22, 14.59) + centre 50.06 (49.99, 50.12) + normalization 169.14 (166.98, 171.78) gaussian_3 - sigma 16.61 (16.57, 16.64) - centre 49.98 (49.94, 50.01) - normalization 206.37 (205.66, 206.90) + sigma 5.69 (5.60, 5.77) + centre 50.02 (49.99, 50.05) + normalization 57.79 (55.90, 59.21) gaussian_4 - sigma 4.02 (2.66, 5.25) - centre 48.04 (47.53, 48.52) - normalization 0.45 (0.35, 0.60) + sigma 1.01 (1.00, 1.01) + centre 50.00 (50.00, 50.01) + normalization 20.32 (20.16, 20.50) instances - <>:28: SyntaxWarning: invalid escape sequence '\s' - <>:28: SyntaxWarning: invalid escape sequence '\s' - /tmp/ipykernel_20726/2299532750.py:28: SyntaxWarning: invalid escape sequence '\s' - plt.ylabel("Normalized Residuals ($\sigma$)") - - -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_2.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_1.png) -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_3.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_2.png) @@ -977,6 +977,23 @@ print(model.info) sigma UniformPrior [54], lower_limit = 0.0, upper_limit = 25.0 +This time the model itself changes. The `centre` is now shared across the five `Gaussian`'s instead of being +independent, so there are 11 free parameters and 1 shared prior where before there were 15 and none. + +That is the difference between the two approaches. Tuning the priors narrowed where the search looks, whereas +assuming a shared `centre` shrinks the model itself, because four of its parameters have stopped existing. + + +```python +af.ModelPlotter(model).figure() +``` + + + +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_45_0.png) + + + We now repeat the model-fit using this updated model where the `centre` of each `Gaussian` is the same. You should again note that the run time of the fit is significantly faster than the previous fits @@ -1001,39 +1018,18 @@ print("The search has finished run - you may now continue the notebook.") The non-linear search has begun running. This Jupyter notebook cell with progress once the search has completed - this could take a few minutes! - 2026-07-11 16:24:09,915 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. - - - 2026-07-11 16:24:09,924 - root - INFO - Output to hard-disk disabled, input a search name to enable. - - - 2026-07-11 16:24:09,924 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). - - - 2026-07-11 16:24:09,930 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores - - - 2026-07-11 16:24:09,961 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search - - + 2026-09-15 00:27:50,332 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. + 2026-09-15 00:27:50,333 - autofit.non_linear.search.abstract_search - INFO - On-the-fly updates of the maximum likelihood model are disabled. Set `updates: iterations_per_quick_update` in config/general.yaml to a finite number of iterations to enable them. + 2026-09-15 00:27:50,351 - root - INFO - Output to hard-disk disabled, input a search name to enable. + 2026-09-15 00:27:50,353 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). + 2026-09-15 00:27:50,365 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores + 2026-09-15 00:27:50,438 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search ~/venv/PyAuto/lib/python3.12/site-packages/dynesty/dynesty.py:194: UserWarning: Specifying slice option while using rwalk sampler does not make sense warnings.warn('Specifying slice option while using rwalk sampler' - - - 2891it [00:16, 175.08it/s, +50 | bound: 355 | nc: 1 | ncall: 15423 | eff(%): 19.131 | loglstar: -inf < 183.400 < inf | logz: 128.070 +/- 0.872 | dlogz: 0.001 > 0.059] - - - - - 2026-07-11 16:24:27,006 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. - - - 2026-07-11 16:24:27,245 - root - INFO - Removing search internal folder. - - - 2026-07-11 16:24:27,341 - root - INFO - Search complete, returning result - - + 3811it [00:39, 96.54it/s, +50 | bound: 475 | nc: 1 | ncall: 19896 | eff(%): 19.455 | loglstar: -inf < 111.574 < inf | logz: 37.693 +/- 1.007 | dlogz: 0.001 > 0.059] + 2026-09-15 00:28:31,346 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. + 2026-09-15 00:28:32,168 - root - INFO - Removing search internal folder. + 2026-09-15 00:28:32,261 - root - INFO - Search complete, returning result The search has finished run - you may now continue the notebook. @@ -1057,7 +1053,7 @@ plt.errorbar( plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -1067,16 +1063,16 @@ plt.close() residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") -plt.ylabel("Normalized Residuals ($\sigma$)") +plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() plt.clf() plt.close() ``` - Bayesian Evidence 128.06963951 - Maximum Log Likelihood 183.39954200 + Bayesian Evidence 37.69286198 + Maximum Log Likelihood 111.57406728 model Collection (N=11) gaussian_0 - gaussian_4 Gaussian (N=3) @@ -1084,44 +1080,38 @@ plt.close() Maximum Log Likelihood Model: gaussian_0 - gaussian_4 - centre 49.999 + centre 50.002 ... [42 lines of output truncated] ... - gaussian_0 + gaussian_0 - gaussian_2 normalization 0.00 (0.00, 0.00) - sigma 18.66 (16.22, 21.61) + gaussian_0 + sigma 8.16 (3.55, 13.13) gaussian_1 - normalization 99.47 (91.47, 107.33) - sigma 12.09 (11.90, 12.30) + normalization 205.17 (204.08, 206.36) + sigma 16.80 (16.70, 16.88) gaussian_2 - normalization 129.93 (124.16, 136.62) - sigma 19.42 (19.11, 19.79) + sigma 15.45 (8.11, 20.96) gaussian_3 - normalization 49.48 (47.84, 51.33) - sigma 5.33 (5.26, 5.43) + normalization 20.70 (20.57, 20.86) + sigma 1.02 (1.01, 1.02) gaussian_4 - normalization 20.24 (20.11, 20.40) - sigma 1.01 (1.01, 1.02) + normalization 72.72 (71.62, 73.86) + sigma 6.24 (6.18, 6.29) instances - <>:28: SyntaxWarning: invalid escape sequence '\s' - <>:28: SyntaxWarning: invalid escape sequence '\s' - /tmp/ipykernel_20726/2299532750.py:28: SyntaxWarning: invalid escape sequence '\s' - plt.ylabel("Normalized Residuals ($\sigma$)") - - -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_2.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_1.png) -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_3.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_2.png) @@ -1198,6 +1188,21 @@ print(model.info) sigma UniformPrior [69], lower_limit = 0.0, upper_limit = 25.0 +The model is back where it started: five `Gaussian`'s with independent parameters and 15 free parameters. +Approaches 1 and 2 each changed the model or its priors, whereas this third approach changes neither, because +searching parameter space more thoroughly is a property of the search and not of the model. + + +```python +af.ModelPlotter(model).figure() +``` + + + +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_55_0.png) + + + __Search__ We again use the nested sampling algorithm Dynesty to fit the model to the data, but now increase the number of live @@ -1236,39 +1241,18 @@ print("The search has finished run - you may now continue the notebook.") The non-linear search has begun running. This Jupyter notebook cell with progress once the search has completed - this could take a few minutes! - 2026-07-11 16:24:27,728 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. - - - 2026-07-11 16:24:27,740 - root - INFO - Output to hard-disk disabled, input a search name to enable. - - - 2026-07-11 16:24:27,741 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). - - - 2026-07-11 16:24:27,748 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores - - - 2026-07-11 16:24:27,986 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search - - + 2026-09-15 00:28:33,065 - autofit.non_linear.search.abstract_search - INFO - Starting non-linear search with 1 cores. + 2026-09-15 00:28:33,066 - autofit.non_linear.search.abstract_search - INFO - On-the-fly updates of the maximum likelihood model are disabled. Set `updates: iterations_per_quick_update` in config/general.yaml to a finite number of iterations to enable them. + 2026-09-15 00:28:33,083 - root - INFO - Output to hard-disk disabled, input a search name to enable. + 2026-09-15 00:28:33,086 - root - INFO - Starting new Dynesty non-linear search (no previous samples found). + 2026-09-15 00:28:33,098 - autofit.non_linear.initializer - INFO - Generating initial samples of model using JAX LH Function cores + 2026-09-15 00:28:33,435 - autofit.non_linear.initializer - INFO - Initial samples generated, starting non-linear search ~/venv/PyAuto/lib/python3.12/site-packages/dynesty/dynesty.py:194: UserWarning: Specifying slice option while using rwalk sampler does not make sense warnings.warn('Specifying slice option while using rwalk sampler' - - - 53790it [04:36, 194.73it/s, +300 | bound: 1197 | nc: 1 | ncall: 274714 | eff(%): 19.711 | loglstar: -inf < 181.059 < inf | logz: 2.310 +/- 0.729 | dlogz: 0.001 > 0.309] - - - - - 2026-07-11 16:29:15,166 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. - - - 2026-07-11 16:29:21,761 - root - INFO - Removing search internal folder. - - - 2026-07-11 16:29:21,853 - root - INFO - Search complete, returning result - - + 41278it [05:42, 120.67it/s, +300 | bound: 933 | nc: 1 | ncall: 211999 | eff(%): 19.640 | loglstar: -inf < 185.495 < inf | logz: 48.341 +/- 0.632 | dlogz: 0.001 > 0.309] + 2026-09-15 00:34:24,612 - autofit.non_linear.search.updater - INFO - Creating latent samples by drawing 100 from the PDF. + 2026-09-15 00:34:29,444 - root - INFO - Removing search internal folder. + 2026-09-15 00:34:29,509 - root - INFO - Search complete, returning result The search has finished run - you may now continue the notebook. @@ -1292,7 +1276,7 @@ plt.errorbar( plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -1302,22 +1286,16 @@ plt.close() residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") -plt.ylabel("Normalized Residuals ($\sigma$)") +plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() plt.clf() plt.close() ``` - <>:28: SyntaxWarning: invalid escape sequence '\s' - <>:28: SyntaxWarning: invalid escape sequence '\s' - /tmp/ipykernel_20726/2299532750.py:28: SyntaxWarning: invalid escape sequence '\s' - plt.ylabel("Normalized Residuals ($\sigma$)") - - - Bayesian Evidence 2.31047815 - Maximum Log Likelihood 181.05850686 + Bayesian Evidence 48.34113705 + Maximum Log Likelihood 185.49517157 model Collection (N=15) gaussian_0 - gaussian_4 Gaussian (N=3) @@ -1325,23 +1303,23 @@ plt.close() Maximum Log Likelihood Model: gaussian_0 - centre 49.241 + centre 49.607 ... [51 lines of output truncated] ... centre 50.00 (49.99, 50.00) - normalization 20.40 (20.23, 20.62) - sigma 1.01 (1.01, 1.02) + normalization 19.98 (19.79, 20.16) + sigma 1.00 (0.99, 1.00) gaussian_2 - centre 78.99 (72.92, 86.73) - normalization 5.85 (4.39, 7.93) - sigma 20.31 (16.88, 23.43) + centre 50.02 (49.96, 50.09) + normalization 55.91 (52.53, 58.98) + sigma 5.51 (5.38, 5.64) gaussian_3 - centre 49.89 (49.85, 49.93) - normalization 53.26 (49.45, 57.62) - sigma 5.49 (5.35, 5.68) + centre 50.56 (50.33, 50.82) + normalization 123.80 (118.58, 129.69) + sigma 13.18 (12.84, 13.50) gaussian_4 - centre 51.10 (50.62, 52.22) - normalization 71.07 (65.35, 76.17) - sigma 11.90 (11.45, 12.53) + centre 37.85 (36.72, 39.34) + normalization 2.13 (1.43, 3.09) + sigma 4.74 (3.73, 5.78) instances @@ -1350,13 +1328,13 @@ plt.close() -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_2.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_1.png) -![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_3.png) +![png](tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_2.png) @@ -1468,8 +1446,3 @@ representation of the data? These are all questions you should be asking yourself before beginning your model-fitting task, but they will become easier to answer as you gain experience with model-fitting and **PyAutoFit**. - - -```python - -``` diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_17_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_17_0.png new file mode 100644 index 0000000..3468da0 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_17_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_23_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_23_0.png deleted file mode 100644 index dee7b8b..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_23_0.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png index 774c6d1..9462a3f 100644 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_25_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_0.png new file mode 100644 index 0000000..0300758 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_1.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_1.png deleted file mode 100644 index 391b9f7..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_27_1.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_29_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_29_0.png new file mode 100644 index 0000000..2b63ab8 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_29_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_0.png new file mode 100644 index 0000000..3468da0 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_2.png deleted file mode 100644 index f7a47ec..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_2.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_3.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_3.png deleted file mode 100644 index ce5d040..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_35_3.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_1.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_1.png new file mode 100644 index 0000000..31f9aec Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_1.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_2.png new file mode 100644 index 0000000..afc670e Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_39_2.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_2.png deleted file mode 100644 index 708562f..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_2.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_3.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_3.png deleted file mode 100644 index 3c86380..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_43_3.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_45_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_45_0.png new file mode 100644 index 0000000..4b257da Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_45_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_1.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_1.png new file mode 100644 index 0000000..179b9eb Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_1.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_2.png new file mode 100644 index 0000000..2690466 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_49_2.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_2.png deleted file mode 100644 index c4c6edc..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_2.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_3.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_3.png deleted file mode 100644 index 07e4c95..0000000 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_53_3.png and /dev/null differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_55_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_55_0.png new file mode 100644 index 0000000..3468da0 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_55_0.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_1.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_1.png new file mode 100644 index 0000000..257901d Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_1.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_2.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_2.png new file mode 100644 index 0000000..6ee2ae2 Binary files /dev/null and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_61_2.png differ diff --git a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_7_0.png b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_7_0.png index 8513174..0281abd 100644 Binary files a/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_7_0.png and b/markdown/chapter_1_introduction/tutorial_4_why_modeling_is_hard_files/tutorial_4_why_modeling_is_hard_7_0.png differ diff --git a/notebooks/chapter_1_introduction/tutorial_4_why_modeling_is_hard.ipynb b/notebooks/chapter_1_introduction/tutorial_4_why_modeling_is_hard.ipynb index a49399d..e9c4d70 100644 --- a/notebooks/chapter_1_introduction/tutorial_4_why_modeling_is_hard.ipynb +++ b/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", diff --git a/notebooks/chapter_1_introduction/tutorial_6_gradients.ipynb b/notebooks/chapter_1_introduction/tutorial_6_gradients.ipynb index b282e2f..463ce79 100644 --- a/notebooks/chapter_1_introduction/tutorial_6_gradients.ipynb +++ b/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", diff --git a/notebooks/chapter_1_introduction/tutorial_7_the_details.ipynb b/notebooks/chapter_1_introduction/tutorial_7_the_details.ipynb index eace05a..3a7e3e2 100644 --- a/notebooks/chapter_1_introduction/tutorial_7_the_details.ipynb +++ b/notebooks/chapter_1_introduction/tutorial_7_the_details.ipynb @@ -47,6 +47,7 @@ "- **Comparing Searches**: Running the same problems with MCMC, nested sampling and gradient descent to see how each is affected, and why comparing searches is a diagnostic in itself.\n", "- **Clipping**: Parameter combinations that are unphysical even when every individual prior is sensible, the NaN likelihoods they produce, and the resample figure of merit each search substitutes.\n", "- **NaN Diagnostics**: Reading the value-NaN and gradient-NaN counters in `search.summary`, and why a finite likelihood does not guarantee a finite gradient.\n", + "- **Hamiltonian Diagnostics**: Reading `n_divergent`, `ess_min` and `mean_acceptance` off a `BlackJAXNUTS` fit, and what each says about whether the samples can be trusted.\n", "- **Unit Cube Vs Physical**: How a search actually sees parameter space through the priors, and why that changes how it explores.\n", "- **Summary**: When these details matter, and what attending to them buys you." ] @@ -1445,8 +1446,109 @@ "the value, because by the time it runs the damage to the derivative is already recorded. Safety has to live where\n", "the `nan` is created, not where it is detected.\n", "\n", - "A good likelihood value does not mean a good gradient!\n", + "A good likelihood value does not mean a good gradient!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "__Hamiltonian Diagnostics__\n", + "\n", + "Tutorial 6 introduced `BlackJAXNUTS`, the Hamiltonian sampler which follows the gradient along a trajectory instead\n", + "of proposing a random step. Trajectories are a more elaborate machine than a random hop, and they come with their\n", + "own ways of going wrong, so NUTS reports a set of diagnostics that no other search in this chapter produces.\n", + "\n", + "They matter here for the same reason the NaN counters did. A NUTS fit that has gone wrong still returns a result\n", + "with error bars on it, and the error bars still look reasonable. The diagnostics are how you find out otherwise.\n", + "\n", + "Let's run a short NUTS fit on the single Gaussian dataset, which is well behaved, so we can see what healthy\n", + "diagnostics look like before describing what unhealthy ones mean." + ] + }, + { + "cell_type": "code", + "metadata": {}, + "source": [ + "model_nuts = af.Collection(gaussian=af.Model(Gaussian))\n", + "\n", + "model_nuts.gaussian.centre = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)\n", + "model_nuts.gaussian.normalization = af.UniformPrior(lower_limit=0.0, upper_limit=100.0)\n", + "model_nuts.gaussian.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=25.0)\n", + "\n", + "search = af.BlackJAXNUTS(\n", + " num_warmup=200, # Steps used to tune the sampler, which are then discarded.\n", + " num_samples=300, # Steps kept as samples of the posterior.\n", + ")\n", + "\n", + "print(\n", + " \"\"\"\n", + " The non-linear search has begun running.\n", + " This Jupyter notebook cell with progress once the search has completed - this could take a few minutes!\n", + " \"\"\"\n", + ")\n", + "\n", + "start = time.time()\n", + "\n", + "result_nuts = search.fit(model=model_nuts, analysis=analysis_x1_jax)\n", + "\n", + "print(f\"BlackJAXNUTS run time: {time.time() - start} seconds\")\n", + "print(\"The search has finished run - you may now continue the notebook.\")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The diagnostics live in the `samples_info` dictionary. Three are worth knowing:\n", + "\n", + "- `n_divergent`: the number of trajectories which \"diverged\", meaning the trajectory left the region the likelihood\n", + " describes instead of following it, and had to be abandoned. A handful is tolerable; many means the steps along the\n", + " trajectory are too large and the samples cannot be 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_nuts = result_nuts.samples\n", + "\n", + "print(\"Diagnostics of the Hamiltonian Monte Carlo fit:\\n\")\n", + "print(f\"Number of divergent trajectories = {samples_nuts.samples_info.get('n_divergent')}\")\n", + "print(f\"Minimum effective sample size = {samples_nuts.samples_info.get('ess_min')}\")\n", + "print(\n", + " f\"Mean acceptance rate = {samples_nuts.samples_info.get('mean_acceptance')}\"\n", + ")" + ], + "outputs": [], + "execution_count": null + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "On this dataset you should see few or no divergences, an effective sample size which is a decent fraction of the 300\n", + "samples drawn, and an acceptance rate near 0.8. That is what a healthy gradient sampler looks like.\n", "\n", + "The failure to watch for is the combination this tutorial has been building towards: divergences climbing while the\n", + "acceptance rate falls, on a model whose parameter space has one of the pathologies we have studied. A degeneracy\n", + "like the flip creates a long narrow ridge, and a trajectory following a ridge whose width changes along its length\n", + "is exactly what makes a trajectory diverge. The diagnostics do not tell you the model is badly parameterized, but\n", + "they are often the first sign of it, and unlike the result itself they do not quietly look fine." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ "__Unit Cube Vs Physical__\n", "\n", "The last detail is the most fundamental, and it changes how you think about every prior you have written.\n", diff --git a/scripts/chapter_1_introduction/images/bad_fit.png b/scripts/chapter_1_introduction/images/bad_fit.png new file mode 100644 index 0000000..2cba25a Binary files /dev/null and b/scripts/chapter_1_introduction/images/bad_fit.png differ diff --git a/scripts/chapter_1_introduction/images/bad_normalized_residual_map.png b/scripts/chapter_1_introduction/images/bad_normalized_residual_map.png new file mode 100644 index 0000000..54f48d5 Binary files /dev/null and b/scripts/chapter_1_introduction/images/bad_normalized_residual_map.png differ diff --git a/scripts/chapter_1_introduction/images/good_fit.png b/scripts/chapter_1_introduction/images/good_fit.png new file mode 100644 index 0000000..1d01c75 Binary files /dev/null and b/scripts/chapter_1_introduction/images/good_fit.png differ diff --git a/scripts/chapter_1_introduction/images/good_normalized_residual_map.png b/scripts/chapter_1_introduction/images/good_normalized_residual_map.png new file mode 100644 index 0000000..84a07be Binary files /dev/null and b/scripts/chapter_1_introduction/images/good_normalized_residual_map.png differ diff --git a/scripts/chapter_1_introduction/images/okay_fit.png b/scripts/chapter_1_introduction/images/okay_fit.png new file mode 100644 index 0000000..5ed0c26 Binary files /dev/null and b/scripts/chapter_1_introduction/images/okay_fit.png differ diff --git a/scripts/chapter_1_introduction/images/okay_normalized_residual_map.png b/scripts/chapter_1_introduction/images/okay_normalized_residual_map.png new file mode 100644 index 0000000..87667c3 Binary files /dev/null and b/scripts/chapter_1_introduction/images/okay_normalized_residual_map.png differ diff --git a/scripts/chapter_1_introduction/tutorial_4_why_modeling_is_hard.py b/scripts/chapter_1_introduction/tutorial_4_why_modeling_is_hard.py index 7e078fa..9aaf41b 100644 --- a/scripts/chapter_1_introduction/tutorial_4_why_modeling_is_hard.py +++ b/scripts/chapter_1_introduction/tutorial_4_why_modeling_is_hard.py @@ -166,8 +166,8 @@ def model_data_from(self, xvalues: np.ndarray) -> np.ndarray: To define the Analysis class for this model-fit, we need to ensure that the `log_likelihood_function` can handle an instance containing multiple 1D profiles. Below is an expanded explanation and the corresponding class definition: -The log_likelihood_function will now assume that the instance it receives consists of multiple Gaussian profiles. -For each Gaussian in the instance, it will compute the model_data and then sum these to create the overall `model_data` +The `log_likelihood_function` will now assume that the instance it receives consists of multiple Gaussian profiles. +For each Gaussian in the instance, it will compute the `model_data` and then sum these to create the overall `model_data` that is compared to the observed data. """ @@ -395,7 +395,7 @@ def model_data_from_instance(self, instance): plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -424,7 +424,7 @@ def model_data_from_instance(self, instance): capsize=2, linestyle="", ) -plt.title(f"Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Residuals") plt.show() @@ -445,7 +445,7 @@ def model_data_from_instance(self, instance): residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() @@ -608,7 +608,7 @@ def model_data_from_instance(self, instance): plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -618,7 +618,7 @@ def model_data_from_instance(self, instance): residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() @@ -725,7 +725,7 @@ def model_data_from_instance(self, instance): plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -735,7 +735,7 @@ def model_data_from_instance(self, instance): residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() @@ -840,7 +840,7 @@ def model_data_from_instance(self, instance): plt.plot(range(data.shape[0]), model_data, color="r") for model_data_1d_individual in model_data_list: plt.plot(range(data.shape[0]), model_data_1d_individual, "--") -plt.title(f"Fit (log likelihood = {result.log_likelihood})") +plt.title(f"Fit (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel("Profile normalization") plt.show() @@ -850,7 +850,7 @@ def model_data_from_instance(self, instance): residual_map = data - model_data normalized_residual_map = residual_map / noise_map plt.plot(xvalues, normalized_residual_map, color="k") -plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood})") +plt.title(f"Normalized Residuals (log likelihood = {result.log_likelihood:.2f})") plt.xlabel("x values of profile") plt.ylabel(r"Normalized Residuals ($\sigma$)") plt.show() diff --git a/scripts/chapter_1_introduction/tutorial_6_gradients.py b/scripts/chapter_1_introduction/tutorial_6_gradients.py index 94746c3..11e3d7a 100644 --- a/scripts/chapter_1_introduction/tutorial_6_gradients.py +++ b/scripts/chapter_1_introduction/tutorial_6_gradients.py @@ -587,7 +587,7 @@ def log_likelihood_from_vector(vector): advanced together in one compiled call using `vmap`, so twelve lanes cost far less than twelve times one lane. This search is given a `name` and `path_prefix`, so its results are written to the `output` folder as tutorial 5 -did, because there is a file in there we are about to read. +did, where you can inspect them once the fit completes. """ search = af.MultiStartAdam( name="tutorial_6_gradients_adam", @@ -639,29 +639,6 @@ def log_likelihood_from_vector(vector): plt.close() """ -Because this search wrote its results to hard disk, its output folder contains a file called `search.summary`: how -long the search took, how long one log likelihood evaluation took and, for gradient searches, diagnostics on how -often things went wrong. We read it back below. -""" -search_summary_path = path.join(str(search.paths.output_path), "search.summary") - -if path.exists(search_summary_path): - with open(search_summary_path) as f: - print(f.read()) - -""" -The block at the bottom, headed `Resampling Info`, contains two entries only a gradient search can report: - -- `Value-NaN Lane-Steps`: the number of times a lane stepped somewhere the log likelihood could not be computed at - all, most often because it stepped outside the priors. - -- `Gradient-NaN Lane-Steps`: the number of times the log likelihood *was* computable but its gradient was not. This - is the sneakier of the two, because such a lane does not crash or die, it simply stops moving while continuing to - look perfectly healthy. - -Both counters are usually small and harmless, but they are the vocabulary you need to diagnose a gradient fit that -has gone quietly wrong. Tutorial 7 explains where they come from and what to do about them. - __Markov Chain Monte Carlo (MCMC)__ In tutorial 3 we used `Emcee`, whose walkers propose a step, compute the likelihood there and accept or reject the @@ -691,18 +668,19 @@ def log_likelihood_from_vector(vector): print(result.info) """ -Now the gradient-aware alternative, `BlackJAXNUTS`, which is Hamiltonian Monte Carlo. The physical picture behind it -is genuinely helpful. +Now the gradient-aware alternative, `BlackJAXNUTS`, which is Hamiltonian Monte Carlo. It is still MCMC, with walkers +moving through parameter space and proposals accepted or rejected, but the walker now knows which way to step, +because autodiff hands it the gradient at every point it visits. -Imagine the likelihood surface turned upside down, so its peak becomes a valley, and place a ball on the resulting -landscape. Give it a random flick and let it roll: it accelerates down slopes, coasts up the other side and travels -a long way while staying in regions the landscape favours. That trajectory is computed from the gradient at each -moment, which is what autodiff hands us for free. Where the ball stops becomes the next sample, and because it -travelled a long, informed distance rather than a small random hop, consecutive samples are far less similar. +That changes how far one proposal can usefully travel. Instead of a single small random hop, the walker follows the +gradient along a trajectory of many small steps, staying in the regions the likelihood favours the whole way, and +ends up somewhere genuinely far from where it started. Consecutive samples are therefore far less correlated than +`Emcee`'s. -"NUTS" stands for the No U-Turn Sampler, which solves the awkward choice here: how long to let the ball roll. Roll -too briefly and you wasted the gradient; roll too long and the ball curves back on itself. NUTS stops the trajectory -when it starts doubling back. Being a gradient method, it needs the JAX analysis. +The one thing left to choose is how long to follow that trajectory. Stop too early and the gradient was wasted; +carry on too long and the path curves back on itself and returns to where it began. "NUTS" stands for the No U-Turn +Sampler, which watches for that doubling back and ends the trajectory there. Being a gradient method, it needs the +JAX analysis. """ search = af.BlackJAXNUTS( num_warmup=200, # Steps used to tune the sampler, which are then discarded. @@ -725,29 +703,6 @@ def log_likelihood_from_vector(vector): print(result.info) -""" -Hamiltonian sampling comes with its own diagnostics, stored in the `samples_info` dictionary. Three are worth -knowing: - -- `n_divergent`: the number of trajectories which "diverged", meaning the ball flew off to infinity instead of - following the landscape. A handful is tolerable; many means the steps are too large and the samples cannot be - trusted. - -- `ess_min`: the "effective sample size" of the worst constrained parameter. Consecutive samples are correlated, so - 300 samples are worth fewer than 300 independent draws, and this says how many they are worth. - -- `mean_acceptance`: the fraction of proposed trajectories accepted, which for NUTS should sit high, around the 0.8 - the warm up phase tunes towards. A low value means the sampler is struggling. -""" -samples = result.samples - -print("Diagnostics of the Hamiltonian Monte Carlo fit:\n") -print(f"Number of divergent trajectories = {samples.samples_info.get('n_divergent')}") -print(f"Minimum effective sample size = {samples.samples_info.get('ess_min')}") -print( - f"Mean acceptance rate = {samples.samples_info.get('mean_acceptance')}" -) - """ Because NUTS maps out the posterior, we can plot the Probability Density Functions of its samples with `corner.py`, wrapped via the `aplt.corner_cornerpy` function, exactly as we did for `Emcee` in tutorial 5. @@ -892,7 +847,7 @@ def log_likelihood_from_vector(vector): r""" __Wrap Up__ -This tutorial took the one sentence tutorial 3 used to describe how an MLE search moves and unpacked it: +Tutorial 3 described how an MLE search moves in a single sentence. This tutorial unpacked that sentence: 1. **Gradients**: the gradient of the log likelihood is a vector with one entry per free parameter, evaluated at a point, pointing in the direction the likelihood increases fastest. @@ -907,8 +862,8 @@ def log_likelihood_from_vector(vector): compiled call, making it far harder to trap in a local maximum. 4. **MCMC**: `Emcee` proposes random steps and accepts or rejects them, never asking which way is up. -`BlackJAXNUTS` rolls a ball across the landscape using the gradient to shape its trajectory, producing far less -correlated samples for the same number of steps. +`BlackJAXNUTS` follows the gradient along a trajectory before taking each sample, producing far less correlated +samples for the same number of steps. 5. **Nested sampling**: `Nautilus` replaces low likelihood live points with higher likelihood ones drawn from the priors, a procedure with no direction of travel for a gradient to inform. JAX speeds it up by evaluating batches of diff --git a/scripts/chapter_1_introduction/tutorial_7_the_details.py b/scripts/chapter_1_introduction/tutorial_7_the_details.py index 5545883..cda17d7 100644 --- a/scripts/chapter_1_introduction/tutorial_7_the_details.py +++ b/scripts/chapter_1_introduction/tutorial_7_the_details.py @@ -42,6 +42,7 @@ - **Comparing Searches**: Running the same problems with MCMC, nested sampling and gradient descent to see how each is affected, and why comparing searches is a diagnostic in itself. - **Clipping**: Parameter combinations that are unphysical even when every individual prior is sensible, the NaN likelihoods they produce, and the resample figure of merit each search substitutes. - **NaN Diagnostics**: Reading the value-NaN and gradient-NaN counters in `search.summary`, and why a finite likelihood does not guarantee a finite gradient. +- **Hamiltonian Diagnostics**: Reading `n_divergent`, `ess_min` and `mean_acceptance` off a `BlackJAXNUTS` fit, and what each says about whether the samples can be trusted. - **Unit Cube Vs Physical**: How a search actually sees parameter space through the priors, and why that changes how it explores. - **Summary**: When these details matter, and what attending to them buys you. """ @@ -1061,7 +1062,80 @@ def model_data_from_instance(self, instance): the `nan` is created, not where it is detected. A good likelihood value does not mean a good gradient! +""" + +r""" +__Hamiltonian Diagnostics__ + +Tutorial 6 introduced `BlackJAXNUTS`, the Hamiltonian sampler which follows the gradient along a trajectory instead +of proposing a random step. Trajectories are a more elaborate machine than a random hop, and they come with their +own ways of going wrong, so NUTS reports a set of diagnostics that no other search in this chapter produces. + +They matter here for the same reason the NaN counters did. A NUTS fit that has gone wrong still returns a result +with error bars on it, and the error bars still look reasonable. The diagnostics are how you find out otherwise. + +Let's run a short NUTS fit on the single Gaussian dataset, which is well behaved, so we can see what healthy +diagnostics look like before describing what unhealthy ones mean. +""" +model_nuts = af.Collection(gaussian=af.Model(Gaussian)) + +model_nuts.gaussian.centre = af.UniformPrior(lower_limit=0.0, upper_limit=100.0) +model_nuts.gaussian.normalization = af.UniformPrior(lower_limit=0.0, upper_limit=100.0) +model_nuts.gaussian.sigma = af.UniformPrior(lower_limit=0.0, upper_limit=25.0) + +search = af.BlackJAXNUTS( + num_warmup=200, # Steps used to tune the sampler, which are then discarded. + num_samples=300, # Steps kept as samples of the posterior. +) + +print( + """ + The non-linear search has begun running. + This Jupyter notebook cell with progress once the search has completed - this could take a few minutes! + """ +) + +start = time.time() + +result_nuts = search.fit(model=model_nuts, analysis=analysis_x1_jax) + +print(f"BlackJAXNUTS run time: {time.time() - start} seconds") +print("The search has finished run - you may now continue the notebook.") + +""" +The diagnostics live in the `samples_info` dictionary. Three are worth knowing: + +- `n_divergent`: the number of trajectories which "diverged", meaning the trajectory left the region the likelihood + describes instead of following it, and had to be abandoned. A handful is tolerable; many means the steps along the + trajectory are too large and the samples cannot be trusted. + +- `ess_min`: the "effective sample size" of the worst constrained parameter. Consecutive samples are correlated, so + 300 samples are worth fewer than 300 independent draws, and this says how many they are worth. + +- `mean_acceptance`: the fraction of proposed trajectories accepted, which for NUTS should sit high, around the 0.8 + the warm up phase tunes towards. A low value means the sampler is struggling. +""" +samples_nuts = result_nuts.samples +print("Diagnostics of the Hamiltonian Monte Carlo fit:\n") +print(f"Number of divergent trajectories = {samples_nuts.samples_info.get('n_divergent')}") +print(f"Minimum effective sample size = {samples_nuts.samples_info.get('ess_min')}") +print( + f"Mean acceptance rate = {samples_nuts.samples_info.get('mean_acceptance')}" +) + +""" +On this dataset you should see few or no divergences, an effective sample size which is a decent fraction of the 300 +samples drawn, and an acceptance rate near 0.8. That is what a healthy gradient sampler looks like. + +The failure to watch for is the combination this tutorial has been building towards: divergences climbing while the +acceptance rate falls, on a model whose parameter space has one of the pathologies we have studied. A degeneracy +like the flip creates a long narrow ridge, and a trajectory following a ridge whose width changes along its length +is exactly what makes a trajectory diverge. The diagnostics do not tell you the model is badly parameterized, but +they are often the first sign of it, and unlike the result itself they do not quietly look fine. +""" + +""" __Unit Cube Vs Physical__ The last detail is the most fundamental, and it changes how you think about every prior you have written. diff --git a/workspace_index.json b/workspace_index.json index 77e2ebe..3bcadab 100644 --- a/workspace_index.json +++ b/workspace_index.json @@ -175,6 +175,7 @@ "Comparing Searches", "Clipping", "NaN Diagnostics", + "Hamiltonian Diagnostics", "Unit Cube Vs Physical", "Summary" ],