# An importance sampling method for Feldman-Cousins confidence intervals

Lukas Berns<EMAIL_ADDRESS>Tohoku University, Japan

###### Abstract

In various high-energy physics contexts, such as neutrino-oscillation
experiments, several assumptions underlying the typical asymptotic confidence
interval construction are violated, such that one has to resort to
computationally expensive methods like the Feldman-Cousins method for
obtaining confidence intervals with proper statistical coverage. By
construction, the computation of intervals at high confidence levels requires
fitting millions or billions of pseudo-experiments, while wasting most of the
computational cost on overly precise intervals at low confidence levels. In
this work, a simple importance sampling method is introduced which reuses
pseudo-experiments produced for all tested parameter values in a single
mixture distribution. This results in a significant error reduction on the
estimated critical values, especially at high confidence levels, and
simultaneously yields a correct interpolation of these critical values between
the parameter values at which the pseudo-experiments were produced. The
theoretically calculated performance is demonstrated numerically using a
simple example from the analysis of neutrino oscillations. The relationship to
similar techniques applied in statistical mechanics and $p$-value computations
is discussed.

##  I Introduction

An essential part of any experiment is the statistical analysis to extract
information about the model parameters, such as physics constants, from the
measurement outcome. As measurements inherently include statistical
fluctuations, one often reports these constraints in the form of confidence
intervals (or confidence regions in higher dimensions). These are intervals
over the parameter space calculated from the observed data, which are
constructed in such a way that for any true value of the parameters, at least
a pre-defined percentage of the possible experimental outcomes would produce
an interval that covers the true parameter value. The pre-defined percentage
over possible experimental outcomes is called the confidence level (CL).

For the rest of this paper we shall use the following notation: $x$ denotes
the experimental outcome, which can be a vector of many observations within
the single experiment. $\theta$ denotes the model parameters, which can
contain one or higher dimensional continuous degrees of freedom, and may
contain discrete degrees of freedom as well. $p(x\mid\theta)$ denotes the
probability distribution function for the experimental outcomes given some
model parameters. $p(x\mid\theta)$ seen as a function of $\theta$ for a given
experimental outcome is called the likelihood function and denoted
$L(\theta\mid x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=p(x\mid\theta)$.
The parameter value for which the likelihood is maximized is denoted
$\hat{\theta}(x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathop{\mathrm{arg~{}max}}_{\theta}L(\theta\mid
x)$, and the difference of the log-likelihood at some parameter value to the
maximum likelihood is denoted as $\Delta\chi^{2}(\theta\mid x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=-2\log
L(\theta\mid x)/L(\hat{\theta}(x)\mid x)$. The confidence level is denoted
$1-\alpha$.

In many cases a useful theorem by Wilks [1] can be applied, which greatly
simplifies the construction of such confidence intervals. The theorem says
that in the asymptotic limit, $\Delta\chi^{2}(\theta\mid x)$ evaluated at the
true parameter value is distributed as a chi-squared distribution with $k$
degrees of freedom, where $k$ is the dimension of the parameter space
$\theta$, which has to be continuous. The theorem holds under suitable
conditions which ensure that a maximum likelihood value can be found in the
neighborhood of the true parameter value with a quadratic Taylor expansion of
the likelihood. Given this asymptotic distribution, one can thus construct a
confidence interval by all values of $\theta$ that satisfy
$\Delta\chi^{2}(\theta\mid x)\leq\Delta\chi^{2}_{c}$, where the critical value
$\Delta\chi^{2}_{c}$ is easily computed from the quantile function of the chi-
squared distribution.

Due to the necessary assumptions, confidence intervals based on Wilks’ theorem
are not suitable if the number of observations is small, or the parameter
space is unsuitable because of physical boundaries (such as $\theta\geq 0$),
discrete degrees of freedom, or periodicities that cannot be captured by the
quadratic expansion. Neutrino oscillation experiments for example suffer from
all of these deficiencies, for which we will present an example later. In this
situation, one has to resort to actually producing ensembles of pseudo-
experiments for selected parameter values to study the distribution of a
suitable statistic to be used for the construction of the confidence interval.

A commonly used method is the Feldman-Cousins (FC) method [2], where for each
pseudo-experiment $x^{\prime}$ generated assuming a true value $\theta_{t}$,
the $\Delta\chi^{2}(\theta_{t}\mid x^{\prime})$ value at the true parameter
value is computed to obtain its distribution. Then the critical value
$\Delta\chi^{2}_{c}$ is obtained by the empirical $1-\alpha$ percentile of
this distribution. Since the distribution of $\Delta\chi^{2}(\theta_{t}\mid
x^{\prime})$ will in general be different for each true parameter value, the
critical values are now a function of the true value at which they are
computed, which we denote as $\Delta\chi^{2}_{c}(\theta_{t})$. Finally, the
confidence interval for the actually observed data $x$ is constructed by
choosing $\Delta\chi^{2}(\theta\mid x)\leq\Delta\chi^{2}_{c}(\theta)$. In
practice, it is only possible to compute $\Delta\chi^{2}_{c}(\theta)$ at
selected parameter values, which need to be interpolated, for example
linearly, in order to compute the confidence intervals.

The Feldman-Cousins method is very inefficient for obtaining high-CL
intervals, because by definition, only a small fraction of pseudo-experiments
contribute to the quantile computation. For example, in particle physics the
threshold for “discovery” is commonly chosen at $\alpha=5.7\times 10^{-7}$
(the “$5\sigma$” threshold), in which case only one in 1.7 million pseudo-
experiments would (by definition) have a $\Delta\chi^{2}(\theta_{t}\mid
x^{\prime})$ value larger than the critical value. As a result, one easily
ends up with millions of pseudo-experiments to be fitted in order to obtain
the necessary critical values, while simultaneously “wasting” most of this
computation time for over-precise critical values at lower CL. In practice, FC
confidence intervals are often computed only up to $2\sigma$
($\alpha=4.6\times 10^{-2}$) or $3\sigma$ CL ($\alpha=2.7\times 10^{-3}$) for
such reasons.

In this work, we show that it is actually extremely easy to introduce an
alternative sampling distribution that generates high-CL pseudo-experiments
much more frequently: one simply reuses the pseudo-experiments generated at
the values of the parameters in the form of a mixture distribution. By
appropriate reweighting, this results in an exponential reduction in the
errors on critical values for high CL. The method also introduces a method for
correctly interpolating the critical values between the subset of true
parameter values, thus removing the need of naive interpolation methods that
are commonly employed.

The paper is organized as follows. First, we review the conventional FC
method. Next we define the new method, deriving it from a discussion of an
ideal importance sampling distribution. Bounds for the importance sampling
weights are calculated, which are used to calculate the reduction of errors on
the estimated critical values compared to the conventional FC method. The
ability to interpolate critical values and the calculation of errors and other
other diagnostics are discussed. Next, a toy example from the analysis of
neutrino oscillations is used to compare the two computation methods and the
improvement is checked against the theoretical upper bounds from the previous
section. Finally, we discuss the relationship to similar techniques in
statistical mechanics and $p$-value calculations, the relationship to Bayesian
marginalized likelihoods, and the limit of applicability in the presence of
nuisance parameters.

##  II Critical values in the conventional Feldman-Cousins method

To prepare the notation, we briefly review the computation of critical values
in the conventional Feldman-Cousins method. First, we make a choice of $S$
points in the parameter space, which we denote $\theta_{s}$ with $s$ going
from $1$ to $S$. At each $\theta_{s}$, we now generate an ensemble of
$n_{\mathrm{exp}}$ pseudo-experiments $\\{x\\}_{s}$ by sampling from
$p(x\mid\theta_{s})$. While all pseudo-experiments are assumed to live in the
same space, the $s$ suffix on the curly brackets representing the ensemble
keeps track of the distribution that generated the experiments. For each
pseudo-experiment $x\in\\{x\\}_{s}$, we now compute
$\Delta\chi^{2}(\theta_{s}\mid x)$ and find the $1-\alpha$ quantile
$\Delta\chi^{2}_{c,s}$ through any suitable estimator. For example, one may
simply sort the $\Delta\chi^{2}(\theta_{s}\mid x)$ values and take the
$\lfloor\alpha\times n_{\mathrm{exp}}\rfloor$ largest value as
$\Delta\chi^{2}_{c,s}$, in which case we have

$\sum_{x\in\\{x\\}_{s}}I\bigl{(}\Delta\chi^{2}(\theta_{s}\mid
x)\geq\Delta\chi^{2}_{c,s}\bigr{)}=\lfloor\alpha\times
n_{\mathrm{exp}}\rfloor.$ (1)

Here, $I(\cdot)$ is the indicator function returning $1$ if the logical
statement in the parentheses is true, and $0$ otherwise. $\lfloor\cdot\rfloor$
denotes the floor function.

Finally, the critical value function $\Delta\chi^{2}_{c}(\theta)$ is obtained
by some interpolation scheme. For example, one may set
$\Delta\chi^{2}_{c}(\theta_{s})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}_{c,s}$
and linearly interpolate for any $\theta$ values in between. To reduce the
interpolation error, one typically has to either manually or automatically [3]
adjust the choice of sampling parameter values $\\{\theta\\}_{S}$ in an
iterative scheme.

The asymptotic variance on the critical values is proportional to the binomial
error $\alpha(1-\alpha)/n_{\mathrm{exp}}$, so high-CL ($\alpha\ll 1$)
generally means that one needs $n_{\mathrm{exp}}\gg 1/\alpha$ for reliable
critical values. Since the whole process is repeated for all $S$ points in the
parameter space, the total number of generated (and fitted) pseudo-experiments
is $S\times n_{\mathrm{exp}}\gg S/\alpha$.

##  III The mixture Feldman-Cousins method

###  III.1 Definition

Our new method, which we shall refer to as the “mixture Feldman-Cousins”
method, differs from the conventional method mainly in the reuse of _all_
generated pseudo-experiments for the critical-value computation of _each_
target parameter space point $\theta_{t}$ with an additional weight

$w(x\mid\theta_{t})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{p(x\mid\theta_{t})}{\frac{1}{S}\sum_{s=1}^{S}p(x\mid\theta_{s})}=\frac{1}{\frac{1}{S}\sum_{s=1}^{S}\exp[-\frac{1}{2}\\{\Delta\chi^{2}(\theta_{s}\mid
x)-\Delta\chi^{2}(\theta_{t}\mid x)\\}]}.$ (2)

The value in the denominator is the sampling probability distribution of
$x\in\\{x\\}_{\mathrm{mix}}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\bigcup_{s=1}^{S}\\{x\\}_{s}$,
which is the mixture distribution of $p(x\mid\theta_{s})$ for all $\theta_{s}$
values. Since the weights are based on the sampling probabilities which are
nothing but the likelihood function, they are computable using the same
procedure that calculates the $\Delta\chi^{2}(\theta\mid x)$ for each pseudo-
experiment. Due to taking the difference of two $\Delta\chi^{2}$ values, the
contribution from the minimum $\chi^{2}$ at $\hat{\theta}(x)$, as well as any
$\theta$-independent offsets (e.g. the $n!$ factor in the Poisson likelihood)
vanish in the denominator and hence do not need to be known accurately. While
in the conventional method one only needs to compute
$\Delta\chi^{2}(\theta_{s}\mid x)$ for the $\theta_{s}$ value at which the
pseudo-experiment was generated, here we need it for all $\theta_{s^{\prime}}$
(including $s^{\prime}\neq s$) and $\theta_{t}$.

Now we can define the critical value $\Delta\chi^{2}_{c,t}$ as the
$w$-weighted $1-\alpha$ quantile of $\Delta\chi^{2}(\theta_{t}\mid x)$ for
$x\sim\\{x\\}_{\mathrm{mix}}$, for example

$\sum_{x\in\\{x\\}_{\mathrm{mix}}}w(x\mid\theta_{t})I\bigl{(}\Delta\chi^{2}(\theta_{t}\mid
x)\geq\Delta\chi^{2}_{c,t}\bigr{)}\lesssim\alpha\times Sn_{\mathrm{exp}},$ (3)

where the $\lesssim$ is meant to represent that we take the smallest
$\Delta\chi^{2}_{c,t}$ that satisfies the inequality.

###  III.2 Derivation

In order to obtain more pseudo-experiments at large $\Delta\chi^{2}$ values,
which would yield more precise high-CL critical values, we use an importance
sampling approach: instead of directly sampling from the target distribution
$p(x\mid\theta_{t})$, we sample from a different distribution and weight the
sampled the toys by the ratio of probability distributions to calculate the
relevant quantities under the target distribution (the critical values). The
question therefore becomes: what is the ideal sampling distribution to
generate the desired pseudo-experiments?

Note that it is important to find a sampling distribution that is as close as
possible to the target distribution apart from generating high
$\Delta\chi^{2}$ pseudo-experiments with higher probability. In particular, if
each experiment $x$ consists of $m$ measurements, the experiments are points
in an $m$-dimensional space and there are $m$ dimensions in which we can
stretch or shrink the sampling distribution. Instead of thinking about
estimating quantiles, let’s think of estimating the probability density
$p(Y(x)\mid\theta_{t})$ using histograms for $Y(x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}(\theta_{t}\mid
x)$. When using reweighting, in addition to the binomial error
$n_{\mathrm{exp}}p(1-p)$ for the number of pseudo-experiments falling into a
bin, there will be an additional contribution due to the variance of weights:
given the estimator

$\hat{P}_{b}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{n_{\mathrm{exp}}}\sum_{i=1}^{n_{\mathrm{exp}}}w(x_{i})I(y_{b}\leq
Y(x_{i})<y_{b+1})$ (4)

and using $I(\cdot)^{2}=I(\cdot)$ we get

$\displaystyle\mathbb{E}[\hat{P}_{b}]$
$\displaystyle=\pi_{b}\mathbb{E}_{b}[w]$ (5)
$\displaystyle\mathrm{\mathop{Var}}[\hat{P}_{b}]$
$\displaystyle=\frac{1}{n_{\mathrm{exp}}}\left(\pi_{b}(1-\pi_{b})\mathbb{E}_{b}[w]^{2}+\pi_{b}\mathrm{\mathop{Var}}_{b}[w]\right)$
(6) $\displaystyle\pi_{b}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{E}[I(y_{b}\leq
Y(x)<y_{b+1})]$ (7) $\displaystyle\mathbb{E}_{b}[w^{k}]$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{\pi_{b}}\mathbb{E}[w(x)^{k}I(y_{b}\leq
Y(x)<y_{b+1})]$ (8) $\displaystyle\mathrm{\mathop{Var}}_{b}[w]$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\mathbb{E}_{b}[w^{2}]-\mathbb{E}_{b}[w]^{2}$
(9)

where in Eq. (6) the first term is the usual binomial error due to the number
of pseudo-experiments falling into the bin, and the second is the additional
term due to the variance of weights among pseudo-experiments falling into the
bin.

We therefore want to increase the number ($n_{\mathrm{exp}}\pi_{b}$) of
pseudo-experiments falling into a high-$y(x)$ bin to reduce the binomial
error, while at the same time keeping the weight-variance within the bin as
small as possible. This means the ideal case of $0$-variance would be for the
weights to depend on $x$ through $y(x)$ alone. Or equivalently, since the
weights are the ratio of the target and sampling distribution, we want to use
a sampling distribution that differs from the target distribution only by a
functional factor of $\Delta\chi^{2}(\theta_{t}\mid x)$.

The key idea is to think about the meaning of a high
$\Delta\chi^{2}(\theta_{t}\mid x)$ value. The likelihood $L(\theta\mid x)$ is
the probability to sample the given pseudo-experiment $x$ from $\theta$. A
high $\Delta\chi^{2}(\theta_{t}\mid x)=-2\log L(\theta_{t}\mid
x)/L(\hat{\theta}(x)\mid x)$ means there exists a value $\hat{\theta}(x)$
where it’s more likely to sample the given pseudo-experiment than at the
“target” $\theta_{t}$ value. Thus by using pseudo-experiments generated at
$\theta\neq\theta_{t}$, we can more efficiently obtain ones with high
$\Delta\chi^{2}(\theta_{t}\mid x)$.

The naive choice of simply using pseudo-experiments generated at some
$\theta^{\prime}$ ($\neq\theta_{t}$) weighted by the ratio of sampling
probabilities $p(x\mid\theta_{t})/p(x\mid\theta^{\prime})$ however will do
worse than before. This is because $\hat{\theta}(x)$ depends on the pseudo-
experiment $x$, such that for some pseudo-experiments it may be more
preferable to sample $x$ from $p(x\mid\theta_{t})$ than from
$p(x\mid\theta^{\prime})$, resulting in an exponentially large (often
unbounded) variance of weights.

The solution is simple: by using a mixture distribution
$p_{\mathrm{sample}}(x)=\frac{1}{S}\sum_{s=1}^{S}p(x\mid\theta_{s})$ over a
set $\\{\theta\\}_{S}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\\{\theta_{1},\theta_{2},\cdots,\theta_{S}\\}$
which includes $\theta_{t}$, we can guarantee the weights to be bounded from
above, because $p_{\mathrm{sample}}(x)\geq\frac{1}{S}p(x\mid\theta_{t})$ and
hence

$w(x\mid\theta_{t})\leq S.$ (10)

###  III.3 Bounds on pseudo-experiment weights for a good grid

If we choose the grid $\\{\theta\\}_{S}$ dense _and_ wide enough (a _good_
grid) such that we may assume to have a good minimum $\hat{\theta}_{S}(x)$ on
$\\{\theta\\}_{S}$ in the sense of

$\displaystyle\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid x)$
$\displaystyle\leq\begin{cases}\epsilon&\text{if\;}\Delta\chi^{2}(\theta_{t}\mid
x)\leq\Delta\chi^{2}_{\mathrm{max}}\\\ \Delta\chi^{2}(\theta_{t}\mid
x)&\text{otherwise}\end{cases}$ (11)
$\displaystyle\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid x)$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\min_{s}\Delta\chi^{2}(\theta_{s}\mid
x)$ (12)

for all $x$, we can put a much stricter bound on the weights than Eq. (10).
Here, $\epsilon\lesssim 1$ will be smaller for denser spacing of
$\\{\theta\\}_{S}$ and $\Delta\chi^{2}_{\mathrm{max}}$ will be larger for a
wider range covered by $\\{\theta\\}_{S}$. Since this additional condition can
deal with the case of $\theta_{t}\notin\\{\theta\\}_{S}$ as well, let us
define a symbol $C$ which is $1$ if $\theta_{t}\in\\{\theta\\}_{S}$ and $0$
otherwise. Note that to guarantee Eq. (11) under $C=0$ one generally needs to
have parameter values in $\\{\theta\\}_{S}$ that surround $\theta_{t}$
sufficiently well. For example, with a 1-dimensional continuous $\theta$
parameter, one needs $\min_{s}\theta_{s}\leq\theta_{t}\leq\max_{s}\theta_{s}$.

First, we focus on the pseudo-experiments with $\Delta\chi^{2}(\theta_{t}\mid
x)\leq\Delta\chi^{2}_{\mathrm{max}}$, which are our primary interest, and note
that

$\frac{p(x\mid\hat{\theta}_{S}(x))}{p(x\mid\theta_{t})}=\exp\big{[}\tfrac{1}{2}\\{\Delta\chi^{2}(\theta_{t}\mid
x)-\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid
x)\\}\big{]}\geq\exp\big{[}\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid
x)-\tfrac{\epsilon}{2}\big{]}.$ (13)

The sum of probability ratios is now bounded from below by

$\sum_{s=1}^{S}\frac{p(x\mid\theta_{s})}{p(x\mid\theta_{t})}\geq
C\times\frac{p(x\mid\theta_{t})}{p(x\mid\theta_{t})}+\frac{p(x\mid\hat{\theta}_{S}(x))}{p(x\mid\theta_{t})}\geq
C+\exp\big{[}\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid
x)-\tfrac{\epsilon}{2}\big{]}$ (14)

because $\\{\theta\\}_{S}$ includes both $\hat{\theta}_{S}(x)$ and (if $C=1$)
$\theta_{t}$. This means for any pseudo-experiment with
$\epsilon<\Delta\chi^{2}(\theta_{t}\mid x)\leq\Delta\chi^{2}_{\mathrm{max}}$,
it is more likely to be sampled in $Sn_{\mathrm{exp}}$ samples from
$\frac{1}{S}\sum_{s=1}^{S}p(x\mid\theta_{s})$ than in $n_{\mathrm{exp}}$
samples from the target distribution $p(x\mid\theta_{t})$. The sum of
probability ratios is further bounded from above by

$\sum_{s=1}^{S}\frac{p(x\mid\theta_{s})}{p(x\mid\theta_{t})}=C+\sum_{s(\neq
t)}\frac{L(\theta_{s}\mid x)}{L(\theta_{t}\mid x)}\leq
C+(S-C)\exp\big{[}\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid x)\big{]}$ (15)

because $L(\theta_{s}\mid x)\leq L(\hat{\theta}(x)\mid x)$. This means the
weights are bounded by

$\frac{S}{C+(S-C)\exp\big{[}\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid
x)\big{]}}\leq
w(x\mid\theta_{t})\leq\frac{S}{C+\exp\big{[}\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid
x)-\tfrac{\epsilon}{2}\big{]}}.$ (16)

We see that the bounds depend on the pseudo-experiments through
$\Delta\chi^{2}(\theta_{t}\mid x)$ only, and also note that for sufficiently
large $\Delta\chi^{2}(\theta_{t}\mid x)$ the ratio of upper $w_{\mathrm{max}}$
to lower bound $w_{\mathrm{min}}$ converges to

$\frac{w_{\mathrm{max}}}{w_{\mathrm{min}}}\to(S-C)e^{\epsilon/2},$ (17)

which indicates a small relative variance of weights as long as the number of
grid points $S$ is not a very large number.

For pseudo-experiments with $\Delta\chi^{2}(\theta_{t}\mid x)$ above the
threshold $\Delta\chi^{2}_{\mathrm{max}}$, we have

$\frac{p(x\mid\hat{\theta}_{S}(x))}{p(x\mid\theta_{t})}=\exp\big{[}\tfrac{1}{2}\\{\Delta\chi^{2}(\theta_{t}\mid
x)-\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid x)\\}\big{]}\geq 1$ (18)

by Eq. (11), and hence an upper bound on the weights

$w(x\mid\theta_{t})\leq\frac{S\times p(x\mid\theta_{t})}{C\times
p(x\mid\theta_{t})+p(x\mid\hat{\theta}_{S}(x))}\leq\frac{S}{C+1}.$ (19)

###  III.4 Critical value estimator performance with a good grid

Since quantiles (the critical values) are just the inverse function of the
cumulative distribution function (CDF), we can estimate the relative reduction
of the quantile estimation variation by the reduction of the CDF estimation
variance. The relationship for an observable $y\sim f(y)$ is given by
$\mathrm{\mathop{Var}}[\hat{y}(P)]=f(y)^{-2}\mathrm{\mathop{Var}}[\hat{P}(y)]$
where $\hat{y}(P)$ is the quantile function estimator, $\hat{P}(y)$ the CDF
estimator, and $f(y)$ the probability distribution function.

Following Eq. (3), and using the shorthand notation $Y(x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}(\theta_{t}\mid
x)$, our CDF estimator is

$\hat{P}(y\mid\theta_{t})=\frac{1}{Sn_{\mathrm{exp}}}\sum_{x\in\\{x\\}_{\mathrm{mix}}}w(x\mid\theta_{t})I(Y(x)\geq
y)$ (20)

with $x\sim\frac{1}{S}\sum_{s=1}^{S}p(x\mid\theta_{s})$. This is an unbiased
estimator for the target CDF $P(y\mid\theta_{t})$

$\displaystyle\mathbb{E}\big{[}\hat{P}(y\mid\theta_{t})\big{]}$
$\displaystyle=\mathbb{E}\big{[}w(x\mid\theta_{t})I(Y(x)\geq y)\big{]}$ (21)
$\displaystyle=\mathbb{E}\bigl{[}I(Y(x)\geq y)\bigm{|}\theta_{t}\bigr{]}$ (22)
$\displaystyle=P(y\mid\theta_{t}).$ (23)

where $\mathbb{E}[\,\cdot\,]$ means to take the expectation with
$x\sim\frac{1}{S}\sum_{s=1}^{S}p(x\mid\theta_{s})$, and
$\mathbb{E}[\;\cdot\mid\theta_{t}]$ to take the expectation with $x\sim
p(x\mid\theta_{t})$. Now defining

$y_{\mathrm{max}}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\max\\{\Delta\chi^{2}_{\mathrm{max}},y\\},$
(24)

the variance from a single pseudo-experiment is

$\displaystyle\mathrm{\mathop{Var}}\big{[}w(x\mid\theta_{t})I(Y(x)\geq
y)\big{]}$ (25)
$\displaystyle=\mathbb{E}\big{[}w(x\mid\theta_{t})^{2}I(Y(x)\geq
y)^{2}\big{]}-\mathbb{E}\big{[}w(x\mid\theta_{t})I(Y(x)\geq y)\big{]}^{2}$
(26) $\displaystyle=\mathbb{E}\big{[}w(x\mid\theta_{t})I(Y(x)\geq
y)\mid\theta_{t}\big{]}-P(y\mid\theta_{t})^{2}$ (27)
$\displaystyle\leq\mathbb{E}\biggl{[}\frac{S\times I(y\leq Y(x)\leq
y_{\mathrm{max}})}{C+\exp\big{[}\tfrac{1}{2}\big{(}Y(x)-\epsilon\big{)}\big{]}}\biggm{|}\theta_{t}\biggr{]}+\frac{S}{C+1}\mathbb{E}\bigl{[}I(Y(x)\geq
y_{\mathrm{max}})\bigm{|}\theta_{t}\bigr{]}-P(y\mid\theta_{t})^{2}$ (28)
$\displaystyle\leq\frac{S}{C+\exp\big{[}\tfrac{1}{2}(y-\epsilon)\big{]}}\mathbb{E}\left[I(y\leq
Y(x)\leq
y_{\mathrm{max}})\mid\theta_{t}\right]+\frac{S}{C+1}P(y_{\mathrm{max}}\mid\theta_{t})-P(y\mid\theta_{t})^{2}$
(29)
$\displaystyle=S\times\left[\frac{P(y\mid\theta_{t})-P(y_{\mathrm{max}}\mid\theta_{t})}{C+\exp\big{[}\tfrac{1}{2}(y-\epsilon)\big{]}}+\frac{P(y_{\mathrm{max}}\mid\theta_{t})}{C+1}\right]-P(y\mid\theta_{t})^{2}$
(30)

where in going from the second to the third line we used
$I(\cdot)^{2}=I(\cdot)$, and in going to the fourth line we used the upper
bound from Eq. (16) and Eq. (19), and in going to the fifth line we used
$Y(x)\geq y$ from the argument of the indicator function. The variance of the
CDF estimator is therefore

$\displaystyle\mathrm{\mathop{Var}}\big{[}\hat{P}(y\mid\theta_{t})\big{]}$
$\displaystyle=\frac{1}{Sn_{\mathrm{exp}}}\mathrm{\mathop{Var}}\big{[}w(x\mid\theta_{t})I(\Delta\chi^{2}(\theta_{t}\mid
x)\geq y)\big{]}$ (31)
$\displaystyle\leq\frac{1}{n_{\mathrm{exp}}}\left(\frac{P(y\mid\theta_{t})-P(y_{\mathrm{max}}\mid\theta_{t})}{C+\exp\big{[}\tfrac{1}{2}(y-\epsilon)\big{]}}+\frac{P(y_{\mathrm{max}}\mid\theta_{t})}{C+1}-\frac{P(y\mid\theta_{t})^{2}}{S}\right)$
(32)

where we note that the $S$ factors in the first two terms were cancelled
thanks to being able to reuse the pseudo-experiments generated at all $S$
values for the CDF estimation of each $\theta_{t}$ value.

For reference, the variance on the CDF estimator in the conventional FC method
(denoted in the following equations by “conv”) is given by the binomial error

$\mathrm{\mathop{Var}}[\hat{P}_{\mathrm{conv}}(y)\mid\theta_{t}]=\frac{1}{n_{\mathrm{exp}}}\left(P(y\mid\theta_{t})-P(y\mid\theta_{t})^{2}\right),$
(33)

so the variance on the estimated critical values $\hat{y}(P\mid\theta_{t})$ in
the mixture-FC method is smaller by the factor

$\displaystyle\gamma$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{\mathrm{\mathop{Var}}[\hat{y}(P\mid\theta_{t})]}{\mathrm{\mathop{Var}}[\hat{y}_{\mathrm{conv}}(P\mid\theta_{t})\mid\theta_{t}]}$
(34)
$\displaystyle=\frac{\mathrm{\mathop{Var}}\big{[}\hat{P}(y\mid\theta_{t})\big{]}}{\mathrm{\mathop{Var}}\big{[}\hat{P}_{\mathrm{conv}}(y\mid\theta_{t})\mid\theta_{t}\big{]}}$
(35)
$\displaystyle\leq\frac{A(y)+B(y)P(y_{\mathrm{max}}\mid\theta_{t})/P-\frac{1}{S}P}{1-P}$
(36) $\displaystyle A(y)$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{C+\exp\big{[}\tfrac{1}{2}(y-\epsilon)\big{]}}$
(37) $\displaystyle B(y)$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{C+1}-A(y).$
(38)

where $y$ is the true $P$-quantile satisfying $P(y\mid\theta_{t})=P$. The
typical functional shape of the upper bound is shown in Fig. 1a. Let us first
consider the case of $P\ll P(y_{\mathrm{max}}\mid\theta_{t})$. For the small
$P\leq 1/2$ values one is typically interested in, the mixture model method
obtains more precise critical values than the conventional method (i.e.
$\gamma\leq 1$) for all $y\geq\epsilon$ if $\theta_{t}\in\\{\theta\\}_{S}$
($C=1$), or all $y\geq\epsilon+2\log 2$ if $\theta_{t}\notin\\{\theta\\}_{S}$
($C=0$). As $y$ increases, the relative variance first decreases linearly, and
for $y\gtrsim 2+\epsilon$ it starts to decrease exponentially as
$\gamma\lesssim\exp(-y/2)$. As $y$ further increases toward
$\Delta\chi^{2}_{\mathrm{max}}$, and $P\ll P(y_{\mathrm{max}}\mid\theta_{t})$
fails to hold anymore, the $B(y)P(y_{\mathrm{max}}\mid\theta_{t})/P$ term
becomes dominant, which saturates to $P(y_{\mathrm{max}}\mid\theta_{t})/P=1$
for $y\geq\Delta\chi^{2}_{\mathrm{max}}$. Hence the improvement flattens out
to $\gamma\leq\frac{1}{C+1}$ for $y\geq\Delta\chi^{2}_{\mathrm{max}}$, which
is still at least as good as the conventional FC method. By choosing suitable
parameter points $\\{\theta\\}_{S}$ and thus a suitable
$\Delta\chi^{2}_{\mathrm{max}}$, critical values of the desired precision can
be calculated. As the exponential reduction in variance cancels the typically
exponential dependence of the CDF on the test-statistic ($\exp(-y/2)$ in the
case of a chi-squared distribution), the relative error on the estimated CDF
becomes approximately flat over a wide range of test-statistic values (Fig.
1b), which is much more efficient than for the conventional FC where low-CL
become over-precise with more pseudo-experiments, while high-CL still suffer
from large errors.

(a) Upper bound on $\gamma$

(b) Relative error on estimated CDF

Figure 1:  (a) Example functional shape of upper bound on the ratio of
estimated critical value variance from Eq. (36). The red line indicates the
error contribution from pseudo-experiments with $y\leq
Y(x)<\Delta\chi^{2}_{\mathrm{max}}$ (first term with $A(y)$) which is
responsible for the exponential reduction of total uncertainty until the
contribution from pseudo-experiments with $y\geq\Delta\chi^{2}_{\mathrm{max}}$
(second term with $B(y)$ shown by green line) takes over for very high CL
critical values. (b) The relative error on the calculated CDF estimator
$\hat{P}(y\mid\theta_{t})$ assuming $n_{\mathrm{exp}}=10,000$ pseudo-
experiments at each sampling value $\theta_{s}$. A reference 10% error
threshold is indicated by the dotted line. The example used for both plots is
constructed assuming $\epsilon=1$, $S=10$, $C=1$,
$\Delta\chi^{2}_{\mathrm{max}}=35$ and the true $Y(x)$-CDF is assumed to be
chi-squared with 1 degree of freedom. In (a), the exponential growth factor
for the green line depends on the assumed CDF, unlike the red line whose decay
factor is given by Eq. (36).

###  III.5 Interpolation

While for the conventional FC method one can only compute the critical values
at the parameter value $\theta_{s}$ where the pseudo-experiments were
generated, in the mixture-FC method it is sufficient to guarantee that the
target parameter value $\theta_{t}$ is sufficiently close and surrounded by
the sampling points $\\{\theta\\}_{S}$ such that condition Eq. (11) holds.
Considering that for a typical setup the toys to be generated are the same as
those used in the conventional FC method, this means that the mixture-FC
method not only reduces the uncertainty on the critical values at the sampling
points $\\{\theta\\}_{S}$, but also allows interpolating the critical values
between these points with similar performance.

###  III.6 Diagnostics and error estimation

As the mixture-FC method exploits the relationship of the $\Delta\chi^{2}$
statistic to the probability of sampling pseudo-experiments, it is essential
that the calculation of $\Delta\chi^{2}$ matches the process used to generate
the pseudo-experiments. It is for example not allowed to sample from a poisson
random number generator while using an approximation like Pearson’s $\chi^{2}$
for $\Delta\chi^{2}$. A simple diagnostic is to calculate the average weight
across all pseudo-experiments in $\\{x\\}_{\mathrm{mix}}$ and check that this
is equal to $1$ up to statistical fluctuations. Since the same pseudo-
experiments will be used for all target parameter values $\theta_{t}$ (which
the weights are a function of), the statistical fluctuations of these average
weights will be correlated for different $\theta_{t}$ values.

To estimate the error of the computed critical values, we recommend using
resampling methods such as the non-parametric bootstrap [4] or jackknife [5]
instead of simple methods like binomial errors, in order to capture not only
the statistical fluctuations in the number of pseudo-experiments that fall
into a range of $\Delta\chi^{2}$ values, but also the statistical fluctuations
in their weights.

##  IV Example with a single cyclic parameter

We consider a simple example that uses a binned-Poisson model, inspired by the
search for CP violation in a long-baseline neutrino oscillation experiment,
here in particular the T2K experiment [6]. The model has a single angular
parameter called the “CP violation phase” $\delta_{\mathrm{CP}}\in[-\pi,\pi]$
which is constrained by $B=10$ Poisson-distributed observations
$n_{b}\sim\mathrm{Poisson}(\lambda_{b})$ with the predicted event rate

$\displaystyle\lambda_{b}(\delta_{\mathrm{CP}})$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=10\times(1-\phi_{b}^{2})\times\left(1-\frac{1}{4}\sin(\delta_{\mathrm{CP}}+\phi_{b})\right)$
(39) $\displaystyle\phi_{b}$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{b-5.5}{10}$
(40)

for each bin with index $b=1,2,\cdots,10$ (Fig. 2). The main feature of this
model is that one is mostly sensitive to $\sin\delta_{\mathrm{CP}}$ through
the overall normalization of approximately 100 total observations
($\sum_{b}n_{b}$), and weakly sensitive to the $\cos\delta_{\mathrm{CP}}$
component through the “shape” of the observations as a function of $b$ (meant
to represent bins of increasing neutrino energy). Deviations from Wilks’
theorem are caused by $\sin\delta_{\mathrm{CP}}$ having physical boundaries at
$\pm 1$ (resulting in reduced critical values around
$\sin\delta_{\mathrm{CP}}=\pm 1$), the sign of $\cos\delta_{\mathrm{CP}}$
acting as an effectively discrete degree of freedom (resulting in increased
critical values at some $\sin\delta_{\mathrm{CP}}\neq\pm 1$ values), as well
as the Poisson nature of the observations. In an actual experiment, one would
have further continuous discrete physics parameters degenerate with
$\delta_{\mathrm{CP}}$ as well as various systematic uncertainties treated as
nuisance parameters. For simplicity and clarity however we focus on
$\delta_{\mathrm{CP}}$ alone, which for continuity with the earlier sections
will be referred to as $\theta=(\delta_{\mathrm{CP}})$, and the observations
as $x=(n_{1},n_{2},\cdots,n_{10})$.

Figure 2: Predicted number of events $\lambda_{b}(\delta_{\mathrm{CP}})$ for
each bin index $b$ as used in the example.

We generate $n_{\mathrm{exp}}=10,000$ pseudo-experiments at each of $S=16$
values of $\theta$ evenly distributed in the parameter range $[-\pi,\pi]$. We
first focus on the target value of $\theta_{t}=-\pi/2$. Fig. 3a shows the
distribution of $\Delta\chi^{2}(\theta_{t}\mid x)$ obtained for pseudo-
experiments $x$ sampled from different $\theta_{s}$ values. In the
conventional FC method, only those generated at $\theta_{s}=\theta_{t}$ are
used, which correspond to the black histogram which falls off quickly for
large $\Delta\chi^{2}(\theta_{t}\mid x)$. In the mixture-FC method we further
make use of the pseudo-experiments generated at all other $\theta_{s}$ values,
of which $\theta_{s}=0$ and the other extreme of $\theta_{s}=\pi/2$ are shown
by the red and green histograms respectively. Clearly, the pseudo-experiments
sampled from the shifted $\theta_{s}$ values have a significantly higher
fraction of large $\Delta\chi^{2}(\theta_{t}\mid x)$ values. At the same time,
one can see one of the problems arising from using only the pseudo-experiments
generated at $\theta_{s}=\pi/2$, in that one would need to apply very large
weights for the small $\Delta\chi^{2}(\theta_{t}\mid x)$ region where
$\theta_{s}=\pi/2$ has a very small sampling probability. The mixture of
pseudo-experiments generated at all 16 $\theta$ values however, shown by the
blue histogram, is able to generate more pseudo-experiments for all
$\Delta\chi^{2}(\theta_{t}\mid x)$ values, with the difference in slope
compared to the black target histogram showing the exponential increase is
pseudo-experiments for larger $\Delta\chi^{2}(\theta_{t}\mid x)$ values. This
is even clearer to see in Fig. 3b where the mixture distribution was
reweighted using the assigned weights. Good agreement with the target
distribution as simulated by the conventional FC method is seen, and the total
number of unweighted pseudo-experiments in the mixture-FC method exceeds the
theoretical lower bound.

(a) $\Delta\chi^{2}_{t}$ at various $\theta_{s}$

(b) Estimated $\Delta\chi^{2}_{t}$ distributions at $\theta_{t}$

Figure 3: $\Delta\chi^{2}_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}(\theta_{t}\mid
x)$ distributions with target $\theta_{t}=-\pi/2$ for (a) various sampling
parameter values $\theta_{s}$, and (b) comparison of estimated distributions
at $\theta_{t}$ obtained using standard FC and mixture-FC methods. In both
plots, error bars indicate $1\sigma$ binomial confidence intervals. Vertical
dashed lines indicate $1,2,3,4,5\sigma$ confidence level critical values
obtained by the mixture-FC method. (a) Error bars are omitted for bins with
zero entries for clarity. (b) The red “weighted mixture” histogram is also
drawn with boxes representing the error from number of pseudo-experiments in
each bin and their weight variance, but these errors are smaller than the line
width and not visible. The “theoretical lower limit” on the total number of
pseudo-experiments in the mixture distribution is obtained by multiplying the
lower bound in Eq. (14) to the red “weighted mixture” histogram assuming $C=1$
and $\epsilon=0.3$.

(a) Distribution of weights

(b) Distribution of $\Delta\chi^{2}(\hat{\theta}_{S}\mid x)$

Figure 4: Distribution of (a) weights $w(x\mid\theta_{t})$ and (b) minimum
$\Delta\chi^{2}$ values over $\\{\theta\\}_{S}$ for each pseudo-experiment as
function of $\Delta\chi^{2}_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}(\theta_{t}\mid
x)$. (a) Dashed blue lines indicate theoretical limits assuming $\epsilon=0.3$
and $C=0$. (b) The dashed line indicates $\epsilon=0.3$. Figure 5: The sample
mean of mixture-FC pseudo-experiment weights $w(x\mid\theta_{t})$ as a
function of the target $\theta_{t}$ value. The error bands indicate the
$1\sigma$ standard error on the mean, which are correlated between different
target $\theta_{t}$ values.

(a) $1$–$5\sigma$ CL

(b) $1\sigma$ CL

(c) $2\sigma$ CL

(d) $3\sigma$ CL

Figure 6: $1,2,3,4,5\sigma$ confidence level critical values (given as
$\sqrt{\Delta\chi^{2}_{c}}$) obtained from the same set of pseudo-experiments
with the standard FC (black error bars) and mixture-FC method (red error
bands). In both cases the error bars/bands indicate the $1\sigma$ error on the
csritical values obtained with binomial/bootstrap errors for the
standard/mixture FC method respectively. Dashed lines indicate critical values
by Wilks’ theorem, which are not valid here, but still drawn for reference.
(a) The black error bars for $4$ and $5\sigma$ have been slightly offset to
prevent overlap. (b,c,d) The gray error bands indicate the linear
interpolation of the error bar end-points in the standard FC method. Figure 7:
Estimated relative errors on the CDF estimator $\hat{P}(y\mid\theta_{t})$ with
target $\theta_{t}=-\pi/2$. For the standard FC the standard error from the
binomial distribution is shown (black solid line), where the more precise CDF
estimate from the mixture-FC method was used in computing these errors. For
the mixture-FC method the bootstrap error estimate (red solid line) is well
below the theoretical upper limit of Eq. (36) calculated assuming
$\epsilon=0.3$ and $\Delta\chi^{2}_{\mathrm{max}}=32$ (blue dashed line).

We now check some of the diagnostics for the mixture-FC method. The
distribution of importance sampling weights $w(x\mid\theta_{t})$ are shown in
Fig. 4a which are found to be mostly a function of
$\Delta\chi^{2}(\theta_{t}\mid x)$ with small additional variance. The weights
are found to be well contained by the theoretical bounds from Eq. (16), which
were drawn assuming a $\epsilon=0.3$ value by looking at the
$\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid x)$ distributions in Fig. 4b. The sum
of weights is found to be consistent with $1$ (Fig. 5).

Next, we look at the critical values. Fig. 6 shows the critical values as
function of the (true/target) parameter value $\theta$ using both the standard
FC method (black error bars) and the mixture-FC method. Despite using the same
set of pseudo-experiments, the critical values obtained with the mixture-FC
method have significantly smaller uncertainty especially at higher CL, and
also provide access to details of the functional shape between the 16 sampling
values of $\theta$.

For the $1\sigma$ critical values (Fig. 6b) we see that despite the relatively
fine spacing of sampling values, the interpolation error as indicated by the
non-overlap of red and gray error bands next to the $\theta=\pm\pi/2$ values
is larger than the size of the binomial error band in the conventional method.
As these binomial error bands do not capture the interpolation error, their
smallness can be misleading, and renders the interpolation feature of the
mixture-FC method very useful.

For the $2\sigma$ (Fig. 6c) and $3\sigma$ critical values (Fig. 6d) we see
good consistency between the two methods while also noting the significantly
smaller errors in the mixture-FC calculation. For $3\sigma$ CL (Fig. 6d) we
see the errors in the conventional method are already so large that some of
the features of the critical values are not recognizable, such as the bumps at
$\theta=0,\pi$ and the asymmetry of critical values for a flip of the
$\sin\delta_{\mathrm{CP}}$ sign, caused by Poisson statistics.

For $4\sigma$ and higher CL the conventional FC method is unable to determine
the critical values except for a lower limit. The mixture-FC method on the
other hand still produces critical values with comparable relative error sizes
to the lower CL critical values.

The estimated relative errors are plotted in Fig. 7 and are consistent with
the typical shape from theoretical arguments (Fig. 1b). To draw the upper
bound from $\gamma$ in Eq. (36), we conservatively assume
$\Delta\chi^{2}_{\mathrm{max}}=32$ based on Fig. 4b, i.e. we will only assume
$\Delta\chi^{2}(\hat{\theta}_{S}(x)\mid x)\leq\epsilon$ up to
$\Delta\chi^{2}(\theta_{t}\mid x)\leq 32$. In this example, the actual
mixture-FC error estimated with the bootstrap is smaller than the theoretical
upper limit from $\gamma$ by about factor 2 for $\Delta\chi^{2}_{t}<16$. This
can be interpreted as more than one sampling value $\theta_{s}$ contributing
to the sampling of each pseudo-experiment, rather than the assumption in the
theoretical upper limit that only $\hat{\theta}_{S}(x)$ would contribute. For
$\Delta\chi^{2}_{t}>16=\Delta\chi^{2}_{\mathrm{max}}/2$ on the other hand the
theoretical upper limit starts to increase significantly, whereas the actual
error estimated with the bootstrap only grows slowly. This can be interpreted
as our choice of $\Delta\chi^{2}_{\mathrm{max}}=32$ being overly conservative:
with the present example, the chosen sampling grid $\\{\theta\\}_{S}$ appears
to be effective up to significantly higher $\Delta\chi^{2}$ values. This is
partly due to the convenient situation of having a parameter
$\theta=(\delta_{\mathrm{CP}})$ with a bounded parameter space
$\delta_{\mathrm{CP}}\in[-\pi,\pi]$.

##  V Discussion

###  V.1 Relation to techniques in statistical mechanics

The presented method is similar in spirit to the “Multiple Histogram
Reweighting” (“multi-histogram”) method [7] in statistical mechanics, where
statistical ensembles are simulated for various parameter values and combined
by reweighting to the desired parameter value. In the multi-histogram method,
the ensembles are combined with an additional per-ensemble weight, which is
adjusted to minimize the overall error on the variable to be estimated. A
similar per-ensemble weighting could be applied in the presented mixture-FC
method as well, where these additional weights would be allowed to depend on
the target $\Delta\chi^{2}_{t}\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\Delta\chi^{2}(\theta_{t}\mid
x)$ value as well, in order to reduce the variance on the critical value
estimator as much as possible.

One difference to the multi-histogram method however, is that because we do
not resort to Markov-Chain Monte-Carlo techniques to sample the pseudo-
experiments, the sampling distribution of pseudo-experiments $x$ at each
parameter value $\theta$ is known exactly including the normalization
constant. Hence the iterative procedure that is required at the end of the
multi-histogram method to self-consistently determine these normalization
constants (the free energies) is not necessary in the mixture-FC method.

###  V.2 Relation to the marginal distribution

The sampling distribution constructed as a mixture over several parameter
values $\\{\theta\\}_{S}$ can be considered a marginal probability
distribution with prior
$\pi(\theta)=\frac{1}{S}\sum_{s=1}^{S}\delta(\theta-\theta_{s})$, where
$\delta(\cdot)$ is the Dirac delta function. Additional per-ensemble weights
as discussed in the previous paragraph would correspond to an alternative
prior
$\pi(\theta\mid\Delta\chi^{2}_{t})=\sum_{s=1}^{S}r_{s}(\Delta\chi^{2}_{t})\,\delta(\theta-\theta_{s})$
where $r_{s}(\Delta\chi^{2}_{t})$ can be optimized to reduce errors subject to
the condition $\sum_{s}r_{s}(\Delta\chi^{2}_{t})=1$ for all
$\Delta\chi^{2}_{t}$. One can even generalize the discussion to continuous
priors $\pi(\theta\mid\Delta\chi^{2}_{t})$, where in order to preserve the
arguments on efficiency reduction, we would need to extend the single-point
condition from Eq. 11 to a condition on a finite-size region on
$\pi(\theta\mid\Delta\chi^{2}_{t})$.

Unlike in the conventional FC method, where one needs a large number of
pseudo-experiments at each target parameter value, it can be preferable in the
mixture-FC method to generate less pseudo-experiments at each sampling value,
but instead increase the number of considered sampling points $S$. If
$Sn_{\mathrm{exp}}$ is held fixed, this results in a reduction of the variance
of critical values by reducing the variance in weights bounded from above by
$\exp(\epsilon/2)$.

Given this relation to the marginal distribution, let us now consider the
computation of $\Delta\chi^{2}(x\mid\theta_{t})=-2\log
L(x\mid\theta_{t})/L(x\mid\hat{\theta}(x))$ as being approximated by $-2\log
L(x\mid\theta_{t})/L_{m}(x)$, where in the denominator, the profiling
operation was replaced by a marginalization over $\theta$ with some prior over
$\theta$. We have therefore a simple likelihood ratio test between
$p(x\mid\theta_{t})$ and $p_{m}(x)\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\int\mathrm{d}\theta\,\pi(\theta)p(x\mid\theta)$
and it now becomes evident that in order to efficiently generate pseudo-
experiments with small $p$-values under the null hypothesis
$p(x\mid\theta_{t})$, one should simply generate the pseudo-experiments from
the alternative hypothesis $p_{m}(x)$, which is what is being done in the
mixture-FC method.

In practice, it will be easier to use the discrete “prior” over
$\\{\theta\\}_{S}$ as was discussed in the text, because unless the likelihood
is gaussian, the numerical integration required for marginalization usually
increases the computational cost and complexity. This relation to the
profiling/marginalization similarities can nevertheless be exploited to
motivate an ideal spacing of $\\{\theta\\}_{S}$ values. Out of the well-known
objective priors, Jeffreys’ prior [8] is known to produce a prior that would
be uniform in the parameterization in which the likelihood is gaussian, if
such a parameterization exists. Since profiling and marginalization with a
uniform prior over a gaussian likelihood produce equivalent results up to a
constant offset, Jeffreys’ prior can be considered a good candidate for
choosing the $\\{\theta\\}_{S}$ values at which to generate pseudo-
experiments. For example, in the CP-violation analysis that was discussed in
the earlier section, it would be more suitable to choose a uniform spacing of
parameter values not in $\delta_{\mathrm{CP}}$ but in
$\sin\delta_{\mathrm{CP}}$ with equal probabilities for the sign of
$\cos\delta_{\mathrm{CP}}$, since the dominant constraint is due to the total
number of events $N\sim\mathrm{Poisson}(\lambda=A+B\sin\delta_{\mathrm{CP}})$
for some constants $A$ and $B$, resulting in an approximately gaussian
likelihood over $\sin\delta_{\mathrm{CP}}$.

###  V.3 Nuisance parameters

Because the significant error-reduction in the mixture-FC method exploits the
specific relation of the $\Delta\chi^{2}(\theta_{t}\mid x)$ statistic to the
distribution that generates the pseudo-experiments, one cannot assume all
features to directly translate to an analysis with nuisance parameters or
“systematic” parameters as they are often called in physics. Especially for
the commonly used methods of profile-FC [9] or posterior Highland-Cousins
methods [10], where the space of nuisance parameters from which to generate
the pseudo-experiments is significantly reduced based on constraints by the
observed experimental data, it is possible to have situations where the
straightforward application of the mixture-FC method does not yield the
exponential reduction of errors on the estimated critical values given by Eq.
36. One should therefore not rely on these to estimate the number of required
pseudo-experiments.

In a relatively general setting, when the target distribution is directly a
part of the mixture distribution (so $C=1$), one can show that even in the
worst case, the variance on the critical values only increases very slightly
compared to the conventional method, by a factor $1/(1-P(y))$ (see Appendix
A). This factor is negligible considering that for high CL we have $P(y)\ll
1$. The weights are bounded from above by a similar limit, which is important
for well defined importance sampling behavior. The naive application of the
mixture-FC method to Feldman-Cousins confidence intervals is therefore still
worth a try. In fact, certain situations may yield near-exponential reduction
of errors as in the case without nuisance parameters, but due the lack of
theoretical guarantees it is suggested to carefully study the distribution of
weights and the reliability of bootstrap error estimates in this situation.

Because one cannot guarantee an exponential reduction of errors in a setting
with nuisance parameters, the ability to interpolate critical values will be
more interesting in this setting. Here it is important that the pseudo-
experiments generated between neighboring $\theta_{s}$ values (and suitable
values of nuisance parameters) sufficiently overlap in the space of pseudo-
experiments. Otherwise, the mismatch between pseudo-experiment generation and
the statistical model behind the test statistic may quickly result in a large
spread of weight values, that would make both the estimated critical values as
well as their error estimates unreliable. This is because with nuisance
parameters, there are significantly more dimensions in which the pseudo-
experiments can differ, even if they have similar values for the test
statistic.

In one specific situation however, all properties discussed in earlier
sections are directly applicable despite the presence of nuisance parameters.
This is when using the prior Highland-Cousins method in conjunction with a
marginal-$\Delta\chi^{2}$ statistic, where it is essential to use the same
prior distribution $\pi(\eta)$ for the nuisance parameters $\eta$ in both
cases. This is because here the effect of nuisance parameters is entirely
absorbed by the probability model to generate the pseudo-experiments, in the
sense of $p(x\mid\theta)=\int\mathrm{d}\eta\,\pi(\eta)\,p(x\mid\theta,\eta)$,
such that as far as the mixture-FC method is concerned, no nuisance parameters
exist.

More detailed discussions with examples and possible modifications to the
sampling distributions for pseudo-experiments will be discussed in a separate
publication.

###  V.4 Relation to similar techniques for statistical inference

Very similar importance sampling techniques have been used for the calculation
of $p$-values under a null hypothesis with a likelihood ratio statistic. For
example, Woodroofe [11] discusses the case with a continuous prior over the
parameter of interest. In our notation,

$p_{\mathrm{sample}}(x)=\int\mathrm{d}\theta\,\pi(\theta)p(x\mid\theta)$ (41)

with only a lower bound on the weights

$w(x\mid\theta_{t})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{p(x\mid\theta_{t})}{p_{\mathrm{sample}}(x)}\geq\frac{p(x\mid\theta_{t})}{p(x\mid\hat{\theta}(x))}=\exp\bigl{[}-\tfrac{1}{2}\Delta\chi^{2}(\theta_{t}\mid
x)\bigr{]}$ (42)

given, rather than an upper bound which would be essential for showing small
errors on the estimated $p$-values. An asymptotic formula for the weights
using the Saddle-point method is also given.

Ref. [12] (Sect. 5.6) describes a method developed in the search for the Higgs
boson by the ATLAS experiment [13]. They point out the difficulty of
performing the integral over the continuous prior in Woodroofe’s method and
instead use a set of discrete points $\\{\theta\\}_{S}$ including
$\theta_{t}$, as we used for the mixture-FC method (with $C=1$). The choice of
weight function however is different in that a pseudo-experiment is used only
if the

$\omega_{s}(x\mid\theta_{t})\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{p(x\mid\theta_{t})}{p(x\mid\theta_{s})}=\exp\bigl{[}-\tfrac{1}{2}\left(\Delta\chi^{2}(\theta_{t}\mid
x)-\Delta\chi^{2}(\theta_{s}\mid x)\right)\bigr{]}$ (43)

value at the parameter value $\theta_{s}$ from which the pseudo-experiment was
sampled from is the smallest among all other values in $\\{\theta\\}_{S}$ —
i.e.
$\omega_{s}(x\mid\theta_{t})=\min_{s^{\prime}}\omega_{s^{\prime}}(x\mid\theta_{t})$
— and discarded otherwise. If the pseudo-experiment is used, it is weighted by
$\omega_{s}(x\mid\theta_{t})$. Then by combining the pseudo-experiments
sampled from all $\\{\theta\\}_{S}$ values with their weights, the desired
distribution $p(x\mid\theta_{t})$ is attained with higher probability to
sample pseudo-experiments of large $\Delta\chi^{2}(\theta_{t}\mid x)$. Since
$\theta_{t}\in\\{\theta\\}_{S}$, this procedure ensures that
$w(x\mid\theta_{t})\leq 1$ for well-behaved weights.

One downside of this vetoing technique, as explained by the authors, is that
the spacing of $\\{\theta\\}_{S}$ must not be too dense in order not to reduce
the efficiency of the method with a high vetoing probability. The mixture-FC
method does not have this problem, because the weights are computed using the
actual sampling probability, which is the sum of probabilities over
$\\{\theta\\}_{S}$, and no vetoing is necessary. While the claimed benefit of
the vetoing technique is its independence from the exact normalization of the
sampling probability distribution — due to only using the probability _ratios_
$\omega_{s}$ — the same is true for the choice of weights in the mixture-FC
method, whose weights from Eq. (2) can be written as

$w(x\mid\theta_{s})=\frac{1}{\frac{1}{S}\sum_{s=1}^{S}\bigl{[}\omega_{s}(x\mid\theta_{t})\bigr{]}^{-1}}.$
(44)

For the problem of finding $p$-values under a null hypothesis with a
likelihood-ratio statistic, the relevant part of the mixture-FC can therefore
be regarded a slight improvement to the method by Ref. [12]. Furthermore, we
have explicitly shown that under suitable conditions, which for a typical
setup requires the absence of nuisance parameters, the variance on the
estimated $p$-values is reduced exponentially for large values of the test
statistic.

Finally, we note some of the differences of computing $p$-values to the FC
confidence interval construction in the context of importance sampling. When
computing $p$-values, we are typically interested in the distribution of the
test statistic under a _single_ null hypothesis. In contrast, in the FC method
we need the test statistic distribution for _all_ plausible parameter values,
which in practice is achieved by computing them for a finite set
$\\{\theta\\}_{S}$, and interpolating in between. The FC construction
therefore benefits from the ability to interpolate critical values with
importance sampling, which is not always of interest in the computation of
$p$-values. In addition, the pseudo-experiments sampled from different
parameter values as required for the construction of the mixture distribution
are already available even in the conventional FC method, making the
transition to the mixture-FC method straightforward.

##  VI Summary

We presented a new method to compute critical values for Feldman-Cousins
confidence intervals. The method is a simple extension of the conventional
method in that the same sets of pseudo-experiments generated at different
parameter values are simply combined with suitable weights. We showed that
this results in a significant reduction of the errors on the critical values,
with exponential reduction for high confidence level critical values, at
almost no additional computational cost. The method was further shown to
enable accurate interpolation of critical values between the parameter values
at which the pseudo-experiments were generated. The theoretically calculated
performance was confirmed using a simple example for the analysis of neutrino
oscillations. While the exponential reduction of errors is currently only
guaranteed for analyses without nuisance parameters, the general technique is
applicable to any analysis making use of the Feldman-Cousins method.

*

##  Appendix A Analysis of critical value variances for generic mixtures

Let us denote the target distribution of pseudo-experiments at $\theta_{t}$ by
$p_{t}(x)$. In a setting with nuisance parameters $\eta$ with probability
distribution $p(x\mid\theta,\eta)$, this could for example be
$p\big{(}x\mid\theta_{t},\hat{\hat{\eta}}(\theta_{t}\mid
x_{\mathrm{obs}})\big{)}$ for the profile-FC method, or
$\int\mathrm{d}\eta\,\pi(\eta\mid
x_{\mathrm{obs}},\theta_{t})\,p(x\mid\theta_{t},\eta)$ in the posterior HC
method, with $\hat{\hat{\eta}}(\theta\mid
x_{\mathrm{obs}})=\mathop{\mathrm{arg~{}min}}_{\eta}\chi^{2}(\theta,\eta\mid
x_{\mathrm{obs}})$ the profile best-fit values and $\pi(\eta\mid
x_{\mathrm{obs}},\theta)$ the posterior distribution for nuisance parameters
conditioned by the target $\theta$ value for a fit to the observed data
$x_{\mathrm{obs}}$. The other pseudo-experiments are sampled from $p_{a}(x)$,
whose distribution we don’t explicitly specify here, but could for example be
a mixture over different $\theta$ and $\eta$ values. The mixture of $N_{t}$
pseudo-experiments sampled from $p_{t}(x)$ and $N_{a}$ pseudo-experiments
sampled from $p_{a}(x)$ weighted by
$w(x)=(N_{t}+N_{a})p_{t}(x)/\big{[}N_{t}p_{t}(x)+N_{a}p_{a}(x)\big{]}$ can be
evaluated analogously to the main text and using the estimators

$\displaystyle\hat{P}(y)$ $\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{N_{t}+N_{a}}\sum_{i=1}^{N_{t}+N_{a}}w(x_{i})I\big{(}Y(x_{i})\geq
y\big{)}$ (45) $\displaystyle\hat{P}_{\mathrm{conv}}(y)$
$\displaystyle\mathrel{\hbox
to0.0pt{\raisebox{1.29167pt}{$\cdot$}\hss}\raisebox{-1.29167pt}{$\cdot$}}=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}I\big{(}Y(x_{i}^{(t)})\geq
y\big{)}$ (46) $\displaystyle\mathbb{E}[\hat{P}(y)]$
$\displaystyle=\mathbb{E}[\hat{P}_{\mathrm{conv}}(y)]=P(y)$ (47) yield a
variance reduction of $\displaystyle\gamma$
$\displaystyle=\frac{\mathrm{\mathop{Var}}[\hat{P}(y)]}{\mathrm{\mathop{Var}}[\hat{P}_{\mathrm{conv}}(y)]}$
(48)
$\displaystyle\leq\frac{\frac{1}{N_{t}}P(y)-\frac{1}{N_{t}+N_{a}}P(y)^{2}}{\frac{1}{N_{t}}P(y)-\frac{1}{N_{t}}P(y)^{2}}$
(49) $\displaystyle\leq\frac{1}{1-P(y)}$ (50) where $\displaystyle x_{i}$
$\displaystyle\sim\frac{N_{t}p_{t}(x)+N_{a}p_{a}(x)}{N_{t}+N_{a}}$ (51)
$\displaystyle x_{i}^{(t)}$ $\displaystyle\sim p_{t}(x).$ (52)

###### Acknowledgements.

We would like to thank Christophe Bronner and Louis Lyons for useful
discussions and connecting us to Kyle Cranmer, whom we would like to thank for
introducing Refs. [11, 12]. This research was supported by JSPS KAKENHI Grant
Number 19J22440.

## References

  * Wilks [1938] S. S. Wilks, The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses, Annals Math. Statist. 9, 60 (1938).
  * Feldman and Cousins [1998] G. J. Feldman and R. D. Cousins, A Unified approach to the classical statistical analysis of small signals, Phys. Rev. D 57, 3873 (1998), arXiv:physics/9711021 .
  * Li  _et al._ [2020] L. Li, N. Nayak, J. Bian, and P. Baldi, Efficient neutrino oscillation parameter inference using Gaussian processes, Phys. Rev. D 101, 012001 (2020), arXiv:1811.07050 [physics.data-an] .
  * Efron [1979] B. Efron, Bootstrap Methods: Another Look at the Jackknife, The Annals of Statistics 7, 1 (1979).
  * Miller [1974] R. G. Miller, The Jackknife–A Review, Biometrika 61, 1 (1974).
  * Abe  _et al._ [2011] K. Abe _et al._ (T2K), The T2K Experiment, Nucl. Instrum. Meth. A 659, 106 (2011), arXiv:1106.1238 [physics.ins-det] .
  * Ferrenberg and Swendsen [1989] A. M. Ferrenberg and R. H. Swendsen, Optimized Monte Carlo analysis, Phys. Rev. Lett. 63, 1195 (1989).
  * Jeffreys [1946] H. Jeffreys, An invariant form for the prior probability in estimation problems, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 186, 453 (1946), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1946.0056 .
  * Acero  _et al._ [2022] M. A. Acero _et al._ (NOvA), The Profiled Feldman-Cousins technique for confidence interval construction in the presence of nuisance parameters, (2022), arXiv:2207.14353 [hep-ex] .
  * Cousins and Highland [1992] R. D. Cousins and V. L. Highland, Incorporating systematic uncertainties into an upper limit, Nucl. Instrum. Meth. A 320, 331 (1992).
  * Woodroofe [2010] M. Woodroofe, Importance sampling and error probabilities, http://www.sfu.ca/~lockhart/richard/banff2010/woodroofe.pdf (2010), Banff Discovery Workshop.
  * Cranmer [2014] K. Cranmer, Practical Statistics for the LHC, in  _2011 European School of High-Energy Physics_ (2014) pp. 267–308, arXiv:1503.07622 [physics.data-an] .
  * Aad  _et al._ [2012] G. Aad _et al._ (ATLAS), Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC, Phys. Lett. B 716, 1 (2012), arXiv:1207.7214 [hep-ex] .