[Template Fitting] Does my fitting result indicate a bias?

[zjkjsd]

Aloha Dr. Dembinski,

In my new practice of template fitting. I constructed a template T and a test dataset D from the same MC sample. I used all the available data to construct T and randomly sampled 50% of the data in each component for D. One might assume that the fitting results would be consistent with the true values of D. However, my results are off from the true values by about 1 $\sigma$. To further investigate this behavior, I performed the fitting 20000 times, with each iteration drawing the test set from a Poisson distribution defined by D (as I learned from the tutorials. Is this the so-called Toy MC?).

@joblib.delayed
def compute_toys(itry,init_yields):
    rng = np.random.default_rng(itry)
    toy = rng.poisson(data)
    cost = Template(toy, edges, temp, method="da")
    m = Minuit(cost, *init_yields, name=samples_name)
    m.limits = (0, None)
    m.migrad();
    return list(m.values) if m.valid else None

# run 20000 toys
ntoys=20000
toys_migrad = joblib.Parallel(-1)(compute_toys(i,init_yields) for i in range(ntoys))
toys_migrad = [v for v in toys_migrad if v is not None]

The true value of D is
0 $D\tau\nu$ = 17474
1 $D^\ast\tau\nu$ = 8124
2 $D^{\ast\ast}\tau\nu$ = 14310
3 $D\ell\nu$ = 12236
4 $D^\ast\ell\nu$ = 8181
5 $D^{\ast\ast}\ell\nu$ = 9371

We can see the means of fitted $D\tau\nu$, $D^\ast\ell\nu$ and $D^{\ast\ast}\ell\nu$ are off by 1 $\sigma$ and the standard deviation of $D^\ast\tau\nu$ is quite large. Is this result a cause for concern?

[zjkjsd]

Here are the template and test datasets.
template_test.zip

# load template
file_path = 'templates1.json'
with open(file_path, encoding="utf-8") as file:
    ws = json.load(file)
electron_channel = ws['channels'][0]
sample_names = [s['name'] for s in electron_channel['samples']]
templates = np.array([s['data'] for s in electron_channel['samples']])

# load test_data
file_path2 = 'templates_subset1.json'
with open(file_path2, encoding="utf-8") as file:
    ws2 = json.load(file)
electron_channel2 = ws2['channels'][0]
data = np.array([s['data'] for s in electron_channel2['samples']])
counts=np.sum(data,axis=0).reshape(41,47)

# remove empty bins in the template and test_data
indices_non_zero = np.where(np.sum(templates,axis=0).reshape(41,47) != 0) # mask for removing empty bins in the template
temp = [i.reshape(41,47)[indices_non_zero] for i in templates]
data = counts[indices_non_zero]
edges = np.linspace(0, len(test_data), len(test_data)+1)
init_yields = [10000]*6

[HDembinski]

If you look into our paper on template fitting, you will see that these fits can be biased if you are not in the asymptotic limit. Try to make sure that the number of counts per bin is above 10 everywhere to avoid bias.

[zjkjsd]

Thank you! This is very helpful, I will try it out and come back!

[HDembinski]

To clarify, the template fit in iminuit without weights can handle bins with few entries ok, but not when you use weights. Our method is stil one of the best performing, but the low statistic case with weights fundamentally cannot be bias-free.

You can workaround the problem by using adaptive binning: make bins wider where the density is low and keep them narrow where the density is sufficient. In higher dimensions, you can try to generate such a binning automatically using a KD-tree, which is implemented in scipy, but it is not trivial to turn such a tree into a histogram and I am not aware of a library which provides histograms based on KDTrees.

I suggest to estimate the bias using Monte-Carlo simulation and then correct the final results obtained from data based on that.