12.11. cost functions comparison: exercise 12.4#
12.11.1. Exercise 12.4#
Starting from samples of pseudo-random numbers generated
according to a Gaussian probability density function,
using the technique of the toy experiments
study the distribution of the minimum of the cost function
used by iminuit in Gaussian fits.
Does it match a \(\chi^2\) distribution for any number of events, for the cases of maximum likelihood, extended maximum likelihood, and least squares?
from matplotlib import pyplot as plt
import numpy as np
from myrand import generate_TCL_ms
from iminuit import Minuit
from scipy.stats import norm
from math import floor, ceil, sqrt
# for the binning choice
def sturges (N) :
return floor (1 + np.log2 (N))
# parameters valid for all tests
N_toys = 5000
sample_size = 500
bins = sturges (sample_size)
12.11.1.1. binned likelihood#
from iminuit.cost import BinnedNLL
def mod_signal_bin_LL (bin_edges, mu, sigma) :
''' fitting model for binned likelihood'''
return norm.cdf (bin_edges, mu, sigma)
means_bin_LL = []
sigmas_bin_LL = []
Q2_bin_LL = []
N_dof_LL = 0
for iToy in range (N_toys) :
subsample = generate_TCL_ms (1., 0.7, sample_size)
bin_content, bin_edges = np.histogram (subsample, bins,
range = (floor (min (subsample)), ceil (max (subsample))))
my_cost_func_bin_LL = BinnedNLL (bin_content, bin_edges, mod_signal_bin_LL)
my_minuit_bin_LL = Minuit (
my_cost_func_bin_LL,
mu = np.mean (subsample),
sigma = np.std (subsample),
)
my_minuit_bin_LL.limits['sigma'] = (0, None)
my_minuit_bin_LL.migrad ()
if not my_minuit_bin_LL.valid : continue
means_bin_LL.append (my_minuit_bin_LL.values['mu'])
sigmas_bin_LL.append (my_minuit_bin_LL.values['sigma'])
Q2_bin_LL.append (my_minuit_bin_LL.fval)
if (iToy == 0) : N_dof_LL = my_minuit_bin_LL.ndof
12.11.1.2. extended binned likelihood#
from iminuit.cost import ExtendedBinnedNLL
def mod_signal_bin_ext_LL (bin_edges, N_signal, mu, sigma) :
''' fitting model for extended binned likelihood'''
return N_signal * norm.cdf (bin_edges, mu, sigma)
means_bin_ext_LL = []
sigmas_bin_ext_LL = []
Q2_bin_ext_LL = []
N_dof_ext_LL = 0
for iToy in range (N_toys) :
subsample = generate_TCL_ms (1., 0.7, sample_size)
bin_content, bin_edges = np.histogram (subsample, bins,
range = (floor (min (subsample)), ceil (max (subsample))))
my_cost_func_bin_ext_LL = ExtendedBinnedNLL (bin_content, bin_edges, mod_signal_bin_ext_LL)
my_minuit_bin_ext_LL = Minuit (
my_cost_func_bin_ext_LL,
N_signal = len (subsample),
mu = np.mean (subsample),
sigma = np.std (subsample),
)
my_minuit_bin_ext_LL.limits['N_signal', 'sigma'] = (0, None)
my_minuit_bin_ext_LL.migrad ()
if not my_minuit_bin_ext_LL.valid : continue
means_bin_ext_LL.append (my_minuit_bin_ext_LL.values['mu'])
sigmas_bin_ext_LL.append (my_minuit_bin_ext_LL.values['sigma'])
Q2_bin_ext_LL.append (my_minuit_bin_ext_LL.fval)
if (iToy == 0) : N_dof_ext_LL = my_minuit_bin_ext_LL.ndof
12.11.1.3. least squares#
from iminuit.cost import LeastSquares
def func_approx (x, N_events, mean, sigma, bin_width) :
return N_events * norm.pdf (x, mean, sigma) * bin_width
means_LS = []
sigmas_LS = []
Q2_bin_LS = []
N_dof_LS = 0
for iToy in range (N_toys) :
subsample = generate_TCL_ms (1., 0.7, sample_size)
bin_content, bin_edges = np.histogram (subsample, bins,
range = (floor (min (subsample)), ceil (max (subsample))))
bin_centres = 0.5 * (bin_edges[1:] + bin_edges[:-1])
sigma_y = [max (sqrt (num), 1.) for num in bin_content]
least_squares = LeastSquares (bin_centres, bin_content, sigma_y, func_approx)
my_minuit_LS = Minuit (least_squares,
N_events = sample_size,
mean = np.mean (subsample),
sigma = np.std (subsample),
bin_width = bin_centres[1] - bin_centres[0]
)
my_minuit_LS.fixed["bin_width", "N_events"] = True
my_minuit_LS.migrad ()
my_minuit_LS.hesse ()
if not my_minuit_LS.valid : continue
means_LS.append (my_minuit_LS.values['mean'])
sigmas_LS.append (my_minuit_LS.values['sigma'])
Q2_bin_LS.append (my_minuit_LS.fval)
if (iToy == 0) : N_dof_LS = my_minuit_LS.ndof
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[4], line 26
19 my_minuit_LS = Minuit (least_squares,
20 N_events = sample_size,
21 mean = np.mean (subsample),
22 sigma = np.std (subsample),
23 bin_width = bin_centres[1] - bin_centres[0]
24 )
25 my_minuit_LS.fixed["bin_width", "N_events"] = True
---> 26 my_minuit_LS.migrad ()
27 my_minuit_LS.hesse ()
28 if not my_minuit_LS.valid : continue
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/iminuit/minuit.py:789, in Minuit.migrad(self, ncall, iterate, use_simplex)
787 t = mutil._Timer(self._fmin)
788 with t:
--> 789 fm = _robust_low_level_fit(
790 self._fcn,
791 self._last_state,
792 replace_none(ncall, 0),
793 self._strategy,
794 self._tolerance,
795 self._precision,
796 iterate,
797 use_simplex,
798 )
800 self._last_state = fm.state
802 self._fmin = mutil.FMin(
803 fm,
804 "Migrad",
(...)
809 t.value,
810 )
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/iminuit/minuit.py:2940, in _robust_low_level_fit(fcn, state, ncall, strategy, tolerance, precision, iterate, use_simplex)
2938 if precision is not None:
2939 migrad.precision = precision
-> 2940 fm = migrad(ncall, tolerance)
2941 strategy = MnStrategy(2)
2942 migrad = MnMigrad(fcn, fm.state, strategy)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/iminuit/cost.py:632, in Cost.__call__(self, *args)
617 def __call__(self, *args: float) -> float:
618 """
619 Evaluate the cost function.
620
(...)
630 float
631 """
--> 632 r = self._value(args)
633 if self.verbose >= 1:
634 print(args, "->", r)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/iminuit/cost.py:2289, in LeastSquares._value(self, args)
2287 def _value(self, args: Sequence[float]) -> float:
2288 y, ye = self._masked.T[self._ndim :]
-> 2289 ym = self._pred(args)
2290 return self._cost(y, ye, ym)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/iminuit/cost.py:2277, in LeastSquares._pred(self, args)
2275 def _pred(self, args: Sequence[float]) -> NDArray:
2276 x = self._masked.T[0] if self._ndim == 1 else self._masked.T[: self._ndim]
-> 2277 ym = self._model(x, *args)
2278 return _normalize_output(ym, "model", self._ndata())
Cell In[4], line 4, in func_approx(x, N_events, mean, sigma, bin_width)
3 def func_approx (x, N_events, mean, sigma, bin_width) :
----> 4 return N_events * norm.pdf (x, mean, sigma) * bin_width
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py:2100, in rv_continuous.pdf(self, x, *args, **kwds)
2098 putmask(output, (1-cond0)+np.isnan(x), self.badvalue)
2099 if np.any(cond):
-> 2100 goodargs = argsreduce(cond, *((x,)+args+(scale,)))
2101 scale, goodargs = goodargs[-1], goodargs[:-1]
2102 place(output, cond, self._pdf(*goodargs) / scale)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/scipy/stats/_distn_infrastructure.py:580, in argsreduce(cond, *args)
576 newargs = [newargs, ]
578 if np.all(cond):
579 # broadcast arrays with cond
--> 580 *newargs, cond = np.broadcast_arrays(*newargs, cond)
581 return [arg.ravel() for arg in newargs]
583 s = cond.shape
File <__array_function__ internals>:200, in broadcast_arrays(*args, **kwargs)
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/numpy/lib/stride_tricks.py:546, in broadcast_arrays(subok, *args)
542 if all(array.shape == shape for array in args):
543 # Common case where nothing needs to be broadcasted.
544 return args
--> 546 return [_broadcast_to(array, shape, subok=subok, readonly=False)
547 for array in args]
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/numpy/lib/stride_tricks.py:546, in <listcomp>(.0)
542 if all(array.shape == shape for array in args):
543 # Common case where nothing needs to be broadcasted.
544 return args
--> 546 return [_broadcast_to(array, shape, subok=subok, readonly=False)
547 for array in args]
File /opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/numpy/lib/stride_tricks.py:349, in _broadcast_to(array, shape, subok, readonly)
346 raise ValueError('all elements of broadcast shape must be non-'
347 'negative')
348 extras = []
--> 349 it = np.nditer(
350 (array,), flags=['multi_index', 'refs_ok', 'zerosize_ok'] + extras,
351 op_flags=['readonly'], itershape=shape, order='C')
352 with it:
353 # never really has writebackifcopy semantics
354 broadcast = it.itviews[0]
KeyboardInterrupt:
12.11.1.4. compare the five estimates#
from scipy.stats import chi2
def plot_histo (ax, sample, ndof, label, color) :
xMin = floor (min (sample))
xMax = ceil (max (sample))
N_bins = sturges (len (sample))
bin_edges = np.linspace (xMin, xMax, N_bins)
ax.set_xlabel ('Q2')
ax.set_ylabel ('events per bin')
ax.hist (sample,
bins = bin_edges,
color = color,
label = label + ', ' + str (ndof),
)
vertical_limits = ax.get_ylim ()
ax.plot ([ndof, ndof], vertical_limits, color = 'black')
ax.legend ()
x = np.arange (xMin, xMax, 0.001)
bin_width = (xMax - xMin) / N_bins
ax.plot (x, bin_width * len (sample) * chi2.pdf (x, df=ndof))
ax.text (0.66 * xMax, len (sample) / 20., str (sum (sample) / len (sample)))
fig, axes = plt.subplots (nrows = 3, ncols = 1)
axes[0].set_title ('compare cost function minimum distributions', size=14)
plot_histo (axes[0], Q2_bin_LS, N_dof_LS, 'LS', 'red')
plot_histo (axes[1], Q2_bin_LL, N_dof_LL, 'LL', 'orange')
plot_histo (axes[2], Q2_bin_ext_LL, N_dof_ext_LL, 'ext LL', 'yellow')
