Python--iminuit相关学习笔记（五）

Draw error bands–(学习笔记，仅供参考)

import numpy as np
from numba_stats import norm
from iminuit import Minuit
from iminuit.cost import ExtendedUnbinnedNLL
import matplotlib.pyplot as plt

rng = np.random.default_rng(1)
x = rng.normal(size=100)

w, xe = np.histogram(x, bins=100, range=(-5, 5))

werr = w ** 0.5
cx = 0.5 * (xe[1:] + xe[:-1])
dx = np.diff(xe)
d = w / dx
derr = werr / dx


def model(x, par):
    return par[0], par[0] * norm.pdf(x, par[1], par[2])

cost = ExtendedUnbinnedNLL(x, model)

m = Minuit(cost, (1, 0, 1))
m.migrad()
m.hesse()

plt.errorbar(cx, d, derr, fmt="o", label="data", zorder=0)

plt.plot(cx, model(cx, m.values)[1], lw=3,
         label="fit")
plt.legend(frameon=False,
           title=f"$n = {m.values[0]:.2f} +/- {m.errors[0]:.2f}$\n"
                 f"$\mu = {m.values[1]:.2f} +/- {m.errors[1]:.2f}$\n"
                 f"$\sigma = {m.values[2]:.2f} +/- {m.errors[2]:.2f}$")

from jacobi import propagate


y, ycov = propagate(lambda p: model(cx, p)[1], m.values, m.covariance)


plt.errorbar(cx, d, derr, fmt="o", label="data", zorder=0)

plt.plot(cx, y, lw=3, label="fit")

yerr_prop = np.diag(ycov) ** 0.5
plt.fill_between(cx, y - yerr_prop, y + yerr_prop, facecolor="C1", alpha=0.5)

plt.legend(frameon=False,
           title=f"$n = {m.values[0]:.2f} +/- {m.errors[0]:.2f}$\n"
                 f"$\mu = {m.values[1]:.2f} +/- {m.errors[1]:.2f}$\n"
                 f"$\sigma = {m.values[2]:.2f} +/- {m.errors[2]:.2f}$")

rng = np.random.default_rng(1)

par_b = rng.multivariate_normal(m.values, m.covariance, size=1000)

y_b = [model(cx, p)[1] for p in par_b]
yerr_boot = np.std(y_b, axis=0)

plt.errorbar(cx, d, derr, fmt="o", label="data", zorder=0)

plt.plot(cx, y, lw=3, label="fit")

plt.fill_between(cx, y - yerr_boot, y + yerr_boot, facecolor="C1", alpha=0.5)

plt.legend(frameon=False,
           title=f"$n = {m.values[0]:.2f} +/- {m.errors[0]:.2f}$\n"
                 f"$\mu = {m.values[1]:.2f} +/- {m.errors[1]:.2f}$\n"
                 f"$\sigma = {m.values[2]:.2f} +/- {m.errors[2]:.2f}$")

par_b = rng.multivariate_normal(m.values, m.covariance, size=50)
y_b = [model(cx, p)[1] for p in par_b]
yerr_boot_50 = np.std(y_b, axis=0)

plt.errorbar(cx, d, derr, fmt="o", label="data", zorder=0)

plt.plot(cx, y, lw=3, label="fit")

plt.fill_between(cx, y - yerr_boot_50, y + yerr_boot_50, facecolor="C1", alpha=0.5)

plt.legend(frameon=False,
           title=f"$n = {m.values[0]:.2f} +/- {m.errors[0]:.2f}$\n"
                 f"$\mu = {m.values[1]:.2f} +/- {m.errors[1]:.2f}$\n"
                 f"$\sigma = {m.values[2]:.2f} +/- {m.errors[2]:.2f}$")

fig, ax = plt.subplots(1, 2, figsize=(12, 5))

plt.sca(ax[0])
plt.plot(cx, y - yerr_prop, "-C0", label="prop")
plt.plot(cx, y + yerr_prop, "-C0", label="prop")
plt.plot(cx, y - yerr_boot, "--C1", label="boot[1000]")
plt.plot(cx, y + yerr_boot, "--C1", label="boot[1000]")
plt.plot(cx, y - yerr_boot_50, ":C2", label="boot[50]")
plt.plot(cx, y + yerr_boot_50, ":C2", label="boot[50]")
plt.legend()
plt.semilogy();

plt.sca(ax[1])
plt.plot(cx, yerr_boot / yerr_prop, label="boot[1000] / prop")
plt.plot(cx, yerr_boot_50 / yerr_prop, label="boot[50] / prop")
plt.legend()
plt.axhline(1, ls="--", color="0.5", zorder=0)
for delta in (-0.1, 0.1):
    plt.axhline(1 + delta, ls=":", color="0.5", zorder=0)
plt.ylim(0.5, 1.5)

Fits with shared parameters

from iminuit import Minuit
from iminuit.cost import UnbinnedNLL
from iminuit.util import describe
from matplotlib import pyplot as plt
import numpy as np
from numba_stats import norm

rng = np.random.default_rng(1)

width = 2.0
data1 = rng.normal(0, width, size=1000)
data2 = rng.normal(5, width, size=1000)

def pdf1(x, μ_1, σ):
    return norm.pdf(x, μ_1, σ)

def pdf2(x, μ_2, σ):
    return norm.pdf(x, μ_2, σ)

lh = UnbinnedNLL(data1, pdf1) + UnbinnedNLL(data2, pdf2)

def plot(cost, xe, minuit, ax, **style):
    signature = describe(cost)
    data = cost.data

    values = minuit.values[signature]
    errors = minuit.errors[signature]

    cx = (xe[1:] + xe[:-1]) / 2

    ym = np.diff(norm.cdf(xe, *values)) * np.sum(w)
    t = []
    for n, v, e in zip(signature, values, errors):
        t.append(f"${n} = {v:.3f} ± {e:.3f}$")
    ax.plot(cx, ym, label="\n".join(t), **style)

m = Minuit(lh, μ_1=1, μ_2=2, σ=1)

fig, ax = plt.subplots(1, 2, figsize=(14, 5))

hists = [np.histogram(lhi.data, bins=50) for lhi in lh]

for lhi, (w, xe), axi in zip(lh, hists, ax):
    cx = (xe[1:] + xe[:-1]) / 2
    axi.errorbar(cx, w, np.sqrt(w), fmt="ok", capsize=0, zorder=0)
    plot(lhi, xe, m, axi, ls="--")

m.migrad()


for lhi, (w, xe), axi in zip(lh, hists, ax):
    plot(lhi, xe, m, axi)
    axi.legend()