  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.
