Examples/iminuit

examples/iminuit/README.md

@@ -0,0 +1,18 @@
+# iminuit Examples
+
+Each sub-directory contains a self-contained example. The order in
+which the examples are to appear is specified in `order.json` (an
+array of directory names in the expected order).
+
+In each example directory you'll find:
+
+* `config.toml` - must conform to the specification outlined here:
+  https://docs.pyscript.net/latest/user-guide/configuration/ This is
+  parsed and ultimately turned into a JSON representation as part of
+  the package's API object.
+* `setup.py` - Python code for contextual and environmental setup,
+  NOT SEEN BY THE END USER, but is run before the `code.py` code is
+  evaluated. Allows us to create useful (IPython) shims, avoid
+  repeating boilerplate and whatnot.
+* `code.py` - the actual code added to the editor which forms the
+  practical example of using the package.

examples/iminuit/least_squares_fit/code.py

@@ -0,0 +1,70 @@
+"""
+A first taste of iminuit: fit a noisy decay curve with a least-squares
+cost function. Minuit will find the best parameters and estimate the
+uncertainties on each one.
+
+Docs: https://scikit-hep.org/iminuit/
+"""
+from IPython.core.display import display, HTML
+
+import numpy as np
+import matplotlib.pyplot as plt
+from iminuit import Minuit
+from iminuit.cost import LeastSquares
+
+rng = np.random.default_rng(2024)
+
+
+heading("Fitting an exponential decay")
+note(
+    "We have noisy measurements of a signal that we believe decays "
+    "exponentially: <code>y = A * exp(-k * t) + c</code>. We want to "
+    "recover the amplitude, decay rate, and baseline."
+)
+
+
+# The model we believe describes the data. The first argument is the
+# independent variable; the rest are the parameters Minuit will fit.
+def decay_model(t, amplitude, rate, baseline):
+    return amplitude * np.exp(-rate * t) + baseline
+
+
+# Generate synthetic measurements with known "true" parameters.
+true_amplitude, true_rate, true_baseline = 5.0, 0.7, 1.2
+t_data = np.linspace(0, 6, 40)
+y_error = np.full_like(t_data, 0.2)
+y_data = (
+    decay_model(t_data, true_amplitude, true_rate, true_baseline)
+    + rng.normal(0, y_error)
+)
+
+# LeastSquares is a ready-made cost function: it knows how to compare
+# the model to data given per-point uncertainties.
+cost = LeastSquares(t_data, y_data, y_error, decay_model)
+
+# Pass starting values for each parameter by name.
+minuit = Minuit(cost, amplitude=1.0, rate=0.1, baseline=0.0)
+minuit.migrad()   # run the MIGRAD minimizer
+minuit.hesse()    # estimate parameter uncertainties
+
+note("Fit results, with one-sigma uncertainties:")
+for name, value, error in zip(
+    minuit.parameters, minuit.values, minuit.errors,
+):
+    note(f"<code>{name} = {value:.3f} &plusmn; {error:.3f}</code>")
+
+note(f"Reduced chi-squared: <code>{minuit.fval / (len(t_data) - 3):.2f}</code>")
+
+# Plot data and fitted curve.
+fig, ax = plt.subplots(figsize=(8, 4))
+ax.errorbar(t_data, y_data, yerr=y_error, fmt="o",
+            color="gray", label="Data", markersize=4)
+t_smooth = np.linspace(t_data.min(), t_data.max(), 200)
+ax.plot(t_smooth, decay_model(t_smooth, *minuit.values),
+        color="crimson", linewidth=2, label="Best fit")
+ax.set_xlabel("t")
+ax.set_ylabel("y")
+ax.set_title("Exponential decay fit with iminuit")
+ax.legend()
+fig.tight_layout()
+display(fig, append=True)

examples/iminuit/least_squares_fit/config.toml

@@ -0,0 +1 @@
+packages = ["iminuit", "numpy", "matplotlib"]

examples/iminuit/least_squares_fit/setup.py

@@ -0,0 +1,40 @@
+"""
+Shim IPython's display API onto PyScript so example code written in a
+Jupyter/IPython idiom runs unmodified in the browser.
+"""
+
+import sys
+import types
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+ipython = types.ModuleType("IPython")
+core = types.ModuleType("IPython.core")
+core_display = types.ModuleType("IPython.core.display")
+core_display.display = display
+core_display.HTML = HTML
+ipython.core = core
+core.display = core_display
+ipython.get_ipython = lambda: None
+ipython.display = core_display
+sys.modules["IPython"] = ipython
+sys.modules["IPython.core"] = core
+sys.modules["IPython.core.display"] = core_display
+sys.modules["IPython.display"] = core_display
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)

examples/iminuit/minos_and_profile/code.py

@@ -0,0 +1,58 @@
+# ---------------------------------------------------------------------
+# Asymmetric errors via MINOS, plus a 1D likelihood profile plot.
+# ---------------------------------------------------------------------
+
+heading("Going beyond symmetric error bars")
+note(
+    "HESSE gives a symmetric Gaussian approximation to the parameter "
+    "uncertainties. MINOS does better when the likelihood is skewed: "
+    "it walks along the cost function until it climbs by a fixed amount, "
+    "yielding asymmetric ± intervals."
+)
+
+
+# A small calibration dataset: signal counts vs. exposure time.
+def power_model(t, normalization, exponent):
+    return normalization * t ** exponent
+
+
+t_data = np.linspace(0.5, 5.0, 12)
+true_norm, true_exp = 3.0, 1.4
+y_error = 0.5 + 0.1 * t_data
+y_data = power_model(t_data, true_norm, true_exp) + rng.normal(0, y_error)
+
+cost = LeastSquares(t_data, y_data, y_error, power_model)
+minuit = Minuit(cost, normalization=1.0, exponent=1.0)
+minuit.limits["normalization"] = (0, None)
+minuit.migrad()
+minuit.hesse()
+minuit.minos()  # asymmetric intervals for all parameters
+
+note("HESSE (symmetric) errors:")
+display(minuit.errors.to_dict(), append=True)
+
+note("MINOS (asymmetric) errors:")
+for name in minuit.parameters:
+    me = minuit.merrors[name]
+    note(
+        f"<code>{name}: {minuit.values[name]:.3f} "
+        f"(+{me.upper:.3f} / {me.lower:.3f})</code>"
+    )
+
+# Visualize the 1D profile of the cost around the exponent parameter.
+# `mnprofile` re-minimizes the other parameters at each point.
+exponent_grid, fcn_values, _ = minuit.mnprofile("exponent", size=40)
+fmin = minuit.fval
+
+fig, ax = plt.subplots(figsize=(8, 4))
+ax.plot(exponent_grid, fcn_values - fmin, color="purple", linewidth=2)
+ax.axhline(1.0, color="gray", linestyle="--",
+           label="Δχ² = 1 (1σ contour)")
+ax.axvline(minuit.values["exponent"], color="black",
+           linestyle=":", label="Best fit")
+ax.set_xlabel("exponent")
+ax.set_ylabel("χ² − χ²_min")
+ax.set_title("Profile likelihood for the 'exponent' parameter")
+ax.legend()
+fig.tight_layout()
+display(fig, append=True)

examples/iminuit/minos_and_profile/config.toml

@@ -0,0 +1 @@
+packages = ["iminuit", "numpy", "matplotlib"]

examples/iminuit/minos_and_profile/setup.py

@@ -0,0 +1,27 @@
+"""Setup for the Minos / profile example (no IPython shim needed)."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from iminuit import Minuit
+from iminuit.cost import LeastSquares
+
+rng = np.random.default_rng(123)

examples/iminuit/order.json

@@ -0,0 +1,5 @@
+[
+    "least_squares_fit",
+    "unbinned_likelihood",
+    "minos_and_profile"
+]

examples/iminuit/unbinned_likelihood/code.py

@@ -0,0 +1,61 @@
+# ---------------------------------------------------------------------
+# Unbinned maximum-likelihood fit of a Gaussian, with parameter limits.
+# ---------------------------------------------------------------------
+
+import numpy as np
+import matplotlib.pyplot as plt
+from iminuit import Minuit
+from iminuit.cost import UnbinnedNLL
+from scipy.stats import norm
+
+rng = np.random.default_rng(7)
+
+heading("Estimating a Gaussian's parameters from samples")
+note(
+    "Imagine a sensor reports 800 measurements that we expect to be "
+    "Gaussian-distributed. UnbinnedNLL fits the probability density "
+    "directly to the individual samples — no histogram needed."
+)
+
+# True (unknown to the fitter) population parameters.
+true_mu, true_sigma = 2.5, 0.8
+samples = rng.normal(true_mu, true_sigma, size=800)
+
+
+# An unbinned NLL needs a *normalized* probability density function.
+def gaussian_pdf(x, mu, sigma):
+    return norm.pdf(x, mu, sigma)
+
+
+cost = UnbinnedNLL(samples, gaussian_pdf)
+
+minuit = Minuit(cost, mu=0.0, sigma=1.0)
+
+# Constrain sigma to be positive: parameter limits are set like a dict.
+minuit.limits["sigma"] = (1e-3, None)
+
+minuit.migrad()
+minuit.hesse()
+
+note("Best-fit parameters:")
+display(minuit.values.to_dict(), append=True)
+note("One-sigma errors:")
+display(minuit.errors.to_dict(), append=True)
+note(
+    f"True values were mu = {true_mu}, sigma = {true_sigma}. "
+    f"Minimum valid: <strong>{minuit.valid}</strong>."
+)
+
+# Overlay the fitted PDF on a histogram of the samples.
+fig, ax = plt.subplots(figsize=(8, 4))
+ax.hist(samples, bins=40, density=True, color="lightsteelblue",
+        edgecolor="white", label="Samples")
+x_grid = np.linspace(samples.min(), samples.max(), 200)
+ax.plot(x_grid, gaussian_pdf(x_grid, *minuit.values),
+        color="darkblue", linewidth=2, label="Fitted Gaussian")
+ax.set_xlabel("x")
+ax.set_ylabel("density")
+ax.set_title("Unbinned maximum-likelihood fit")
+ax.legend()
+fig.tight_layout()
+display(fig, append=True)

examples/iminuit/unbinned_likelihood/config.toml

@@ -0,0 +1 @@
+packages = ["iminuit", "numpy", "matplotlib", "scipy"]

examples/iminuit/unbinned_likelihood/setup.py

@@ -0,0 +1,20 @@
+"""Setup for the unbinned likelihood example (no IPython shim needed)."""
+import js
+from pyscript import window, HTML, display as _display
+
+js.alert = window.alert
+
+
+def display(*args, **kwargs):
+    return _display(
+        *args, **kwargs, target=__pyscript_display_target__,
+    )
+
+
+def heading(text, level=2):
+    display(HTML(f"<h{level}>{text}</h{level}>"), append=True)
+
+
+def note(text):
+    display(HTML(f"<p>{text}</p>"), append=True)
+