Notebook originel: minimization.ipynb

Minimization

Illustrate chi2-minimization procedures using different librairies.

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

Generative model & Dataset

def model(x, a, b, c):
    return a + b * (x - c)**2

n = 21
pmod = a, b, c = 1, 2, 0  # true values
p0 = 0, 1, 0              # initial guess

rng = np.random.default_rng(seed=1)  # Random Number Generator with a fixed seed (for reproducibility)

data_x = np.linspace(-1, 1, n)
data_y = model(data_x, *pmod)             # Noiseless model

data_yerr = rng.lognormal(mean=-1, sigma=0.5, size=n)  # Log-normal distributed errors
data_y += data_yerr * rng.normal(size=n)               # Add Gaussian noise

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"Model")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"Model: {pmod}");

scipy.optimize.curve_fit

scipy.optimize.curve_fit is the standard minimization procedure of scipy. It is easy to use, but not necessarily the most functionnal nor the more robust.

NB: See curve_fit_utils for some useful addition tools.

from scipy.optimize import curve_fit
from scipy.stats import distributions

pfit, pcov, infodict, msg, ier = curve_fit(
    model,             # model(x, *params)
    data_x, data_y,    # data x and y
    p0=p0,             # initial guesses
    sigma=data_yerr, absolute_sigma=True,  # absolute (not relative) error on y
    full_output=True)  # output additional infodict, msg, ier
print(f"Status {ier}: {msg}")  # ALWAYS check the fit status
if ier in [1, 2, 3, 4]:  # Successful fit
    perr = np.diag(pcov)**0.5  # Diagonal errors on parameters
    print(f"Results: {pfit} ± {perr}")
else:
    print(f"Fit did not converge [{ier}]: {msg}")

Status 1: Both actual and predicted relative reductions in the sum of squares
  are at most 0.000000
Results: [ 1.09715969  1.55898978 -0.02976324] ± [0.11931261 0.24954372 0.04049687]

chi2 = np.sum(((data_y - model(data_x, *pfit)) / data_yerr)**2)  # Final chi2
ndof = len(data_x) - len(pfit)                                   # Degrees of Freedom
pvalue = distributions.chi2.sf(chi2, ndof)  # p-value of the final chi2 / DoF
zscore = distributions.norm.isf(pvalue)     # equivalent z-score

print(f"χ²/DoF = {chi2:.2f}/{ndof}, p-value={pvalue:.3g} ({zscore:.1f} σ)")

χ²/DoF = 25.76/18, p-value=0.105 (1.3 σ)

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"{pmod}")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
ax.plot(data_x, model(data_x, *pfit),
        label=f"{np.round(pfit, 2)} ± {np.round(perr, 2)}")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"χ²/DoF={chi2:.2f}/{ndof} ({zscore:.1f} σ)");

astropy.modeling

astropy.modeling offers some convenient model-building tools, but is based on standard scipy minimization routines.

from astropy import modeling, table, __version__
print("astropy version:", __version__)

astropy version: 5.2.2

@modeling.models.custom_model  # Convert model to astropy
def amodel(x, a=0, b=0, c=0):
    return model(x, a, b, c)

fitter = modeling.fitting.LevMarLSQFitter()  # Levenberg-Marquardt least-squares
fit = fitter(amodel(), data_x, data_y, weights=1/data_yerr)  # chi2 (weighted lsq)

# Parameter errors (https://github.com/astropy/astropy/issues/7202)
errors = amodel().copy()
modeling.fitting.fitter_to_model_params(errors, np.diag(fitter.fit_info['param_cov']**0.5))

/tmp/ipykernel_615550/446977705.py:3: RuntimeWarning: invalid value encountered in sqrt
  modeling.fitting.fitter_to_model_params(errors, np.diag(fitter.fit_info['param_cov']**0.5))

table.Table([fit.param_names, fit.parameters, errors.parameters],
             names=['Names', 'Value', 'Error'])

Table length=3

Names	Value	Error
str1	float64	float64
a	1.0971596940951862	0.14273697695016538
b	1.5589897790695104	0.2985360115620135
c	-0.029763242123276447	0.048447521510153586

chi2 = np.sum(((data_y - fit(data_x))/data_yerr)**2)
ndof = len(data_x) - len(fit.parameters)
pvalue = distributions.chi2.sf(chi2, ndof)
zscore = distributions.norm.isf(pvalue)

print(f"χ²/DoF = {chi2:.2f}/{ndof}, p-value={pvalue:.3g} ({zscore:.1f} σ)")

χ²/DoF = 25.76/18, p-value=0.105 (1.3 σ)

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"{pmod}")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
ax.plot(data_x, fit(data_x),
        label=f"a={np.round(fit.parameters, 2)}±{np.round(errors.parameters, 2)}")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"χ²/DoF={chi2:.2f}/{ndof:.0f} ({zscore:.1f} σ)");

iminuit

iminuit is a jupyter-friendly Python interface to the Minuit2 C++ library used in particle physics. It has more functionalities than scipy, and is generally more robust.

from iminuit import Minuit, __version__
print("iminuit version:", __version__)
from iminuit.cost import LeastSquares  # standard LSQ
import jacobi  # Fast numerical derivatives (uncertainty propagation)
print("jacobi version:", jacobi.__version__)

iminuit version: 2.21.3
jacobi version: 0.8.1

chi2fn = LeastSquares(data_x, data_y, data_yerr, model)  # chi2 objective function
m = Minuit(chi2fn, a=pmod[0], b=pmod[1], c=pmod[2])      # initial guess
m.migrad()  # finds minimum of chi2fn function
m.hesse()   # accurately computes uncertainties

Migrad
FCN = 25.76 (χ²/ndof = 1.4)		Nfcn = 86
EDM = 2.8e-07 (Goal: 0.0002)
Valid Minimum		No Parameters at limit
Below EDM threshold (goal x 10)		Below call limit
Covariance	Hesse ok	Accurate	Pos. def.	Not forced

	Name	Value	Hesse Error	Minos Error-	Minos Error+	Limit-	Limit+	Fixed
0	a	1.10	0.12
1	b	1.56	0.25
2	c	-0.03	0.04

	a	b	c
a	0.0142	-0.023 (-0.773)	-0.0004 (-0.093)
b	-0.023 (-0.773)	0.0623	0.0023 (0.228)
c	-0.0004 (-0.093)	0.0023 (0.228)	0.00164

for param in m.parameters:
    print(f"{param}: {m.values[param]:+.3f} ± {m.errors[param]:.3f}")

a: +1.097 ± 0.119
b: +1.559 ± 0.250
c: -0.030 ± 0.040

pvalue = distributions.chi2.sf(m.fval, m.ndof)
zscore = distributions.norm.isf(pvalue)

print(f"χ²/DoF = {m.fval:.2f}/{m.ndof}, p-value={pvalue:.3g} ({zscore:.1f} σ)")

χ²/DoF = 25.76/18.0, p-value=0.105 (1.3 σ)

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"{pmod}")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
ax.plot(data_x, model(data_x, *m.values),
        label=f"a={np.round(m.values, 2)}±{np.round(m.errors, 2)}")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"χ²/DoF={m.fval:.2f}/{m.ndof:.0f} ({zscore:.1f} σ)");

Confidence interval

fig, axs = m.draw_mnmatrix(cl=[1, 2, 3])  # CI at 1, 2, 3 sigmas
# Add initial (true) parameters
axs[1][0].plot([m.init_params['a'].value], [m.init_params['b'].value], 'k*')
axs[2][0].plot([m.init_params['a'].value], [m.init_params['c'].value], 'k*')
axs[2][1].plot([m.init_params['b'].value], [m.init_params['c'].value], 'k*');

Uncertainty propagation

Propagate uncertainties (and correlations) on parameters to uncertainties on model using 1st-order (linear) propagation (with jacobi library).

NB: see also luprox based on JAX.

%%timeit

y, ycov = jacobi.propagate(lambda p: model(data_x, *p), m.values, m.covariance)
yerr = np.diag(ycov)**0.5

3.19 ms ± 640 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

y, ycov = jacobi.propagate(lambda p: model(data_x, *p), m.values, m.covariance)
yerr_p = np.diag(ycov)**0.5  # Errors on model

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"Model {pmod}")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
l, = ax.plot(data_x, model(data_x, *m.values),
             label=f"Fit {np.round(m.values, 2)}±{np.round(m.errors, 2)}")
ax.fill_between(data_x, y - yerr_p, y + yerr_p,
                facecolor=l.get_color(), alpha=0.5, label="Uncertainty")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"χ²/DoF={m.fval:.2f}/{m.ndof:.0f}");

Bootstrap uncertainties

Estimate uncertainties on model using bootstrapping (random sampling of parameters).

%%timeit

p_boot = rng.multivariate_normal(m.values, m.covariance, size=500)
# standard deviation of bootstrapped curves
yerr = np.std([ model(data_x, *p) for p in p_boot ], axis=0)

2.18 ms ± 82.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

p_boot = rng.multivariate_normal(m.values, m.covariance, size=500)
# standard deviation of bootstrapped curves
yerr_b = np.std([ model(data_x, *p) for p in p_boot ], axis=0)

fix, ax = plt.subplots()
ax.plot(data_x, model(data_x, *pmod), 'k:', label=f"Model {pmod}")
ax.errorbar(data_x, data_y, data_yerr, fmt="o", label="Observations")
l, = ax.plot(data_x, model(data_x, *m.values),
             label=f"Fit {np.round(m.values, 2)}±{np.round(m.errors, 2)}")
ax.fill_between(data_x, y - yerr_b, y + yerr_b,
                facecolor=l.get_color(), alpha=0.5, label="Uncertainty")
ax.legend()
ax.set(xlabel="x", ylabel="y", title=f"χ²/DoF={m.fval:.2f}/{m.ndof:.0f}");

fix, ax = plt.subplots()
ax.plot(data_x, yerr_p, label="Error prop.")
ax.plot(data_x, yerr_b, label="Bootstrap")
ax.set(xlabel="x", ylabel="Model error", title="Comparison of uncertainties")
ax.legend();

Cette page a été générée à partir de minimization.ipynb.