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

Template fits: comparing two chi-square distributed test statistics

from iminuit import cost, Minuit
import numpy as np
from matplotlib import pyplot as plt
from scipy.stats import chi2
from IPython.display import display

xr = (0, 2)  # xrange
rng = np.random.default_rng(1)

nmc = 1000
trials = 1000

data = {}
data2 = {}

first = True
for trial in range(trials):
    for bins in (20, 200,):
        xdata = rng.normal(1, 0.1, size=1000)
        ydata = rng.exponential(size=len(xdata))
        xmix = np.append(xdata, ydata)
        xmix = xmix[(xr[0] < xmix) & (xmix < xr[1])]

        n, xe = np.histogram(xmix, bins=bins, range=xr)

        x = rng.normal(1, 0.1, size=nmc)
        y = rng.exponential(size=nmc)
        t = [
            np.histogram(x, bins=bins, range=xr)[0],
            np.histogram(y, bins=bins, range=xr)[0],
        ]
        c = cost.Template(n, xe, t)
        m = Minuit(c, 1, 1)
        m.migrad()
        assert m.valid
        assert m.accurate
        data.setdefault(bins, []).append(m.fmin.fval)
        data2.setdefault(bins, []).append(np.nansum(c.pulls(m.values) ** 2))
        # display one example fit
        if first:
            display(m)
    first = False

for key in tuple(data):
    val = data[key]
    data[key] = np.array(val)
    val = data2[key]
    data2[key] = np.array(val)

for bins in data:
    plt.figure()
    plt.title(f"bins = {bins}")
    plt.hist(chi2(bins-2).cdf(data[bins]), bins=10, range=(0, 1), label="cost function")
    plt.hist(chi2(bins-2).cdf(data2[bins]), bins=10, range=(0, 1), alpha=0.5, label="sum of pulls squared")
    plt.legend()