[Template Fitting] Is it ok to have large amount of empty bins in the template?

[zjkjsd]

Aloha iminuit community and experts,

I am a student working on R(D) and R(D*) measurements at Belle II. The current plan of signal extraction is through a 2d template fitting ( $|\vec{p_D}| + |\vec{p_\ell}|$ vs. $M_{miss}^2$ ). In my practice using iminuit template fitting, templates are constructed for 6 different signal-like B meson decay modes from signal MC.

These 2d histograms are flattened and then used as templates. However, as seen in the picture, the 2 fitting quantities ( $|\vec{p_D}| + |\vec{p_\ell}|$ vs. $M_{miss}^2$ ) are correlated and there is a large number of empty bins in each template. I thought these empty bins might introduce unnecessary degrees of freedom if the templates are used as is. It might be better if these empty bins can be removed before using them as templates.

So I tested this idea by fitting templates to an independently produced MC sample (test sample). My code snippet and fitting results with and without empty bins are attached below. The empty bins are defined in the test sample.

from iminuit import Minuit
from iminuit.cost import poisson_chi2, Template
from IPython.display import display

xedge = np.linspace(0, len(data), len(data)+1)
init_yields = [12000]*6
c4 = Template(data, xedge, temp, method="da")
m4 = Minuit(c4, *init_yields, name=model.config.samples)
m4.limits = (0, None)
m4.migrad()
m4.hesse()
m4.minos()

Above is the result with empty bins and below is without. The reduced $\chi^2$ changed from 0.4 to 1.2 as expected but the uncertainties increased quite a lot. Now I am not sure if the empty bins should be removed or not. Could you share some insights? Any comments or suggestions are much appreciated!

(ps. I don't know why the Name column in the 2nd block has a lot of extra spaces...)
Mahalo,
Boyang

[HDembinski]

tl;dr: You should keep the empty bins and it makes sense that the fit with empty bins included gives you smaller uncertaintes.

Empty bins in the data are completely ok, no need to cut them. Problematic are empty bins in the templates, more precisely, situations in which all templates have zero entries in some bin but the data bin is non-zero. Such bins have to be discarded in the likelihood function, because it is not possible to draw any information from them (all predictions happen to be empty, so there is no way to estimate the amplitudes).

If at least one template is filled for a given bin and the data bin is empty, the Poisson-based template fit can draw information from that. For the Poisson distribution we can compute the probability to observe nothing for a given expected value. This is why your fit with empty bins included gives you more precise estimates for the yields.

[HDembinski]

That looks like a bug, the Template class should not modify the temp array.

[HDembinski]

I can reproduce this. This must be a bug somewhere in the Template class. This is very strange, because for triggering the bug, it is enough to do this

Template(test_data, xedge, temp, method="da")
c = Template(test_data, xedge, temp, method="da")
m = Minuit(c, *init_yields, name=sample_names)
m.limits = (0, None)
m.simplex().migrad(ncall=1000000)
m.hesse()

[HDembinski]

The bug is fixed in #856, I will release a new iminuit version with this fix later today. I am sorry that this bug slipped through, despite all the testing that I do.

[zjkjsd]

Thank you very much for your help!

[HDembinski]

Thank you for providing your data and code, so I could reproduce the problem. If I find the time, I will use this to investigate the numerical issues, perhaps there is a way to make the likelihood more stable using Kahan summation.

[HDembinski]

Yes, you need to set limits. If you don't do that, the fit may chose a combination of amplitudes where one or more bins get an expected count that is negative, which is mathematically not possible. It is your responsibility to set the parameter limits so that this never happens.

This tutorial explains the problem in more detail.
https://iminuit.readthedocs.io/en/stable/notebooks/scipy_and_constraints.html

Unless your are in the rase case where you need to fit an interference pattern, the amplitudes of components are always positive.

In a future version of iminuit, limits for the amplitudes may be set automatically by default, so that you have to explicitly unset them if you really want negative amplitudes.

[zjkjsd]

Thank you for your response! I wonder why we should set limits for template fitting but not for the example in the tutorial. Is it because the template fitting is essentially a least-square fitting? (I am not sure if this is true).

Another question: I want to apply a constraint:
$$\frac{yield(D\tau\nu)}{yield(D^\ast\tau\nu)}=\frac{yield(D\ell\nu)}{yield(D^\ast\ell\nu)}* \frac{efficiency(\tau\nu)}{efficiency(\ell\nu)}$$
The efficiency of $\tau\nu$ and $\ell\nu$ can be calculated from signal MC. Could you please share some insights about how to implement this constraint? Thank you very much for your help!

[HDembinski]

In a template fit with two templates, you compute the expected rate mu in bin i as mu[i] = a * t1[i] + b * t2[i], where t1 and t2 are the templates. Mathematically, mu[i] must be positive, because the expectation in a bin cannot be negative. t1 and t2 are always non-negative in every bin (bins where both t1 and t2 are zero are discarded by the likelihood), because they are derived from histograms. So there are two ways to make mu[i] positive everywhere. 1) You use a minimizer in which you can straight-forwardly define the constraint "mu[i] must be positive everywhere". Minuit cannot do this, but Scipy can. Note that one of a and b can be negative in this case, but not both. 2) You apply a stronger constraint, namely that a and b must be both positive. That is a constraint that Minuit supports.

[zjkjsd]

Thank you so much for explaining! It is clearer to me now! I am trying to follow the scipy example in a template fit but am stuck on setting up the constraint. I would appreciate it if you could take a look at it in the code below. And perhaps give me some hints on the yield efficiency constraint in my last question if you have time, thank you very much!

from scipy.optimize import NonlinearConstraint
c3 = Template(data, xedge, temp, method="da") # 'jsc', 'asy'
m3 = Minuit(c3, *init_yields, name=samples_name)
m3.limits = (None, None)
m3.simplex().scipy(constraints=NonlinearConstraint(lambda m3: c3.prediction(m3.values)[0], 0, np.inf))
m3.hesse()

[HDembinski]

The function in the constraint must accept the same number of arguments as your cost function. This should work:
m3.simplex().scipy(constraints=NonlinearConstraint(lambda *args: c3.prediction(*args)[0], 0, np.inf))

[zjkjsd]

Thank you again! This got my scipy trial to work! There is a small typo in the code above, it should be:
m3.simplex().scipy(constraints=NonlinearConstraint(lambda *args: c3.prediction(args)[0], 0, np.inf)).
Note the prediction(args). For users who come across this question in the future, this is because the c3.prediction() only takes a Sequence[float] argument.

[HDembinski]

Off-topic: @zjkjsd You may like this PR #858