Probabilistic treatment of the uncertainty from the finite size of weighted Monte Carlo data

TL;DR

This paper introduces a probabilistic framework that accounts for the finite size of Monte Carlo samples in high-energy physics, reducing bias in parameter estimation by replacing standard likelihoods with new generalized probability distributions.

Contribution

It develops analytic generalized probability distributions that incorporate Monte Carlo statistical uncertainties, applicable to various likelihood functions in high-energy physics analyses.

Findings

Reduces bias in parameter estimation from finite Monte Carlo samples

Provides analytic formulas involving Lauricella functions and Dirichlet averages

Demonstrates improved accuracy in a toy Monte Carlo normalization problem

Abstract

Parameter estimation in HEP experiments often involves Monte-Carlo simulation to model the experimental response function. A typical application are forward-folding likelihood analyses with re-weighting, or time-consuming minimization schemes with a new simulation set for each parameter value. Problematically, the finite size of such Monte Carlo samples carries intrinsic uncertainty that can lead to a substantial bias in parameter estimation if it is neglected and the sample size is small. We introduce a probabilistic treatment of this problem by replacing the usual likelihood functions with novel generalized probability distributions that incorporate the finite statistics via suitable marginalization. These new PDFs are analytic, and can be used to replace the Poisson, multinomial, and sample-based unbinned likelihoods, which covers many use cases in high-energy physics. In the limit of infinite statistics, they reduce to the respective standard probability distributions. In the general case of arbitrary Monte Carlo weights, the expressions involve the fourth Lauricella function $F_D$, for which we find a new finite-sum representation in a certain parameter setting. The result also represents an exact form for Carlson's Dirichlet average $R_n$ with $n>0$, and thereby an efficient way to calculate the probability generating function of the Dirichlet-multinomial distribution, the extended divided difference of a monomial, or arbitrary moments of univariate B-splines. We demonstrate the bias reduction of our approach with a typical toy Monte Carlo problem, estimating the normalization of a peak in a falling energy spectrum, and compare the results with previously published methods from the literature.