一、Least-squares fits
import numpy as np
def model(x, a, b):
return a + b * x ** 2
rng = np.random.default_rng(4)
truth = 1, 2
x = np.linspace(0, 1, 20)
yt = model(x, *truth)
ye = 0.4 * x**5 + 0.1
y = rng.normal(yt, ye)
plt.plot(x, yt, ls="--", label="truth")
plt.errorbar(x, y, ye, fmt="ok", label="data")
plt.legend(loc="upper left")
c = cost.LeastSquares(x, y, ye, model)
m1 = Minuit(c, a=0, b=0)
m1.migrad()
m1.visualize()
plt.plot(c.x, model(c.x, *truth), ls="--", label="truth")
def model2(x_y, a, bx, by):
x, y = x_y
return a + bx * x + by * y
# generate a regular grid in x and y
x = np.linspace(-1, 1, 10)
y = np.linspace(-1, 1, 10)
X, Y = np.meshgrid(x, y)
x = X.flatten()
y = Y.flatten()
# model truth
Z = model2((x, y), 1, 2, 3)
# add some noise
rng = np.random.default_rng(1)
Zerr = 1
Z = rng.normal(Z, Zerr)
plt.scatter(x, y, c=Z)
plt.colorbar()
二、Robust least-squares
c.y[3] = 3
m3 = Minuit(c, a=0, b=0)
m3.migrad()
m3.visualize()
plt.plot(c.x, model(c.x, 1, 2), ls="--", label="truth")
from resample.bootstrap import variance as bvar
def fit(x, y, ye):
c = cost.LeastSquares(x, y, ye, model, loss="soft_l1")
m = Minuit(c, a=0, b=0)
m.migrad()
return m.values
berr = bvar(fit, c.x, c.y, c.yerror, size=1000, random_state=1) ** 0.5
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
for i, axi in enumerate(ax):
axi.errorbar(0, m1.values[i], m1.errors[i], fmt="o")
axi.errorbar(1, m3.values[i], m3.errors[i], fmt="o")
axi.errorbar(2, m3.values[i], berr[i], fmt="o")
axi.set_xticks(np.arange(3), ("no outlier", "Minuit, soft_l1", "bootstrap"))