#  Statistical significance testing for mixed priors:
a combined Bayesian and frequentist analysis

Jakob Robnik<EMAIL_ADDRESS>Uroš Seljak<EMAIL_ADDRESS>Physics Department, University of California at Berkeley,
Berkeley, California 94720, USA  Lawrence Berkeley National Laboratory,
Berkeley, California 94720, USA

###### Abstract

In many hypothesis testing applications, we have mixed priors, with well-
motivated informative priors for some parameters but not for others. The
Bayesian methodology uses the Bayes factor and is helpful for the informative
priors, as it incorporates Occam’s razor via multiplicity or trials factor in
the Look Elsewhere Effect. However, if the prior is not known completely, the
frequentist hypothesis test via the false positive rate is a better approach,
as it is less sensitive to the prior choice. We argue that when only partial
prior information is available, it is best to combine the two methodologies by
using the Bayes factor as a test statistic in the frequentist analysis. We
show that the standard frequentist likelihood-ratio test statistic corresponds
to the Bayes factor with a non-informative Jeffrey’s prior. We also show that
mixed priors increase the statistical power in frequentist analyses over the
likelihood ratio test statistic. We develop an analytic formalism that does
not require expensive simulations using a statistical mechanics approach to
hypothesis testing in Bayesian and frequentist statistics. We introduce the
counting of states in a continuous parameter space using the uncertainty
volume as the quantum of the state. We show that both the p-value and Bayes
factor can be expressed as energy versus entropy competition. We present
analytic expressions that generalize Wilks’ theorem beyond its usual regime of
validity and work in a non-asymptotic regime. In specific limits, the
formalism reproduces existing expressions, such as the p-value of linear
models and periodograms. We apply the formalism to an example of exoplanet
transits, where multiplicity can be more than $10^{7}$. We show that our
analytic expressions reproduce the p-values derived from numerical
simulations.

##  I Introduction

The nature of scientific discovery typically proceeds via falsification of the
null hypothesis via a test that is guided by an alternative hypothesis we wish
to compare to. In many cases, we can write the alternative as an extension of
the null hypothesis, such that, for example there is a parameter of the
alternative hypothesis whose value is fixed under the null hypothesis. In
these situations, the null hypothesis is well specified in terms of the prior
distributions of its parameters, while we may have little or no idea what the
prior distribution of parameters of the alternative hypothesis is. In
particular, we may have little or no idea what the prior distribution of the
amplitude of the new effect should be. However, we may be performing the
search over many parameters, some of which have a well-specified prior (so-
called coordinate parameters). In this paper, we are interested in these mixed
prior situations.

Bayesian methodology of hypothesis testing compares the ratio of marginal
likelihoods of the two hypotheses to form a Bayes factor jeffreyOriginalBook .
If the prior distribution for the alternative is known, this is a valid
methodology that yields optimal results. If the prior is only partially known,
the resulting answer is sensitive to the features of the model that we have
little control on. In particular, we can always arbitrarily down-weight the
alternative hypothesis in the Bayes factor by choosing a very broad prior for
its amplitude parameters. This may lead to rejecting the alternative when
using a Bayes factor. This is not justifiable when the prior is poorly known:
the end result can be a missed discovery opportunity. The opposite situation
can also occur, where the chosen prior and the corresponding Bayes factor are
too optimistic for the alternative hypothesis.

This is of most relevance when testing a new model (such as a new theory or a
new phenomenon) that has never been observed before: a test of a new model
should therefore not depend on the reported Bayes factor alone. An alternative
approach is that of a frequentist hypothesis testing, which is defined in
terms of the false positive rate of the null hypothesis (p-value, or Type I
error). It is independent of the validity of the assumed prior for the
alternative hypothesis: we can rule out the null hypothesis even in the
absence of a well developed alternative hypothesis. However, to do so we need
a test statistic, and the best test statistic is the one developed based on
expected properties of alternative hypothesis.

While frequentist false positive rate of the null hypothesis is a well-defined
approach to hypothesis testing, the Bayesian approach has its advantages too.
There are parameters other than amplitude that have well-specified priors, for
which Bayesian hypothesis testing has Occam’s razor built-in
OccamsRazorBayesian , and automatically accounts for effects such as the Look
Elsewhere Effect LEE1 . For example, if we scan for a signal with a large
template bank, we must account for the trials factor, and the corresponding
p-value significance is increased Miller1981 ; multhyptest . There is no
single established frequentist procedure to do this, while Bayes factor
automatically accounts for it. In this case, the dependence on the prior in
Bayes factor analysis is an asset. In many realistic situations, the
alternative hypothesis has several parameters, some with well-specified priors
and some without. This paper aims to address these situations from both
Bayesian and frequentist perspectives and propose a solution that takes
advantage of both methodologies.

One typical example is searching for a signal in a time-series, such as a
localized feature (e.g., planet transit in stellar flux): the signal could be
anywhere in the time series, and we must pay the multiplicity penalty for
looking at many time stamps, and Bayes factor automatically accounts for it.
On the other hand, while the signal can have any amplitude, and we may not
know what amplitudes are reasonable, it is unreasonable to pay a high
multiplicity penalty for a very broad amplitude prior. The Bayesian
methodology does not have a good way to handle this prior dependence, but
frequentist methodology does. We would like to develop a method that preserves
the positive features of the two schools and removes the undesired features.

The main themes of this paper are:

  * •

The Bayes factor is the test statistic with the highest power and should be
used even in frequentist analyses, assuming some of the priors are informative
(known).

  * •

For many applications, some of the priors are known, and some are unknown.
This mixed prior information requires an analysis that combines Bayes factor
with frequentist p-value.

  * •

For mixed prior problems with some unknown priors frequentist p-value or Type
I error (false positive rate) evaluated on the Bayes factor is a better way to
summarize the significance of the alternative hypothesis than the Bayes factor
itself. In many situations, this can be done analytically without the
simulations.

The Neyman-Pearson lemma guarantees that the likelihood ratio is the highest
power test statistic for a simple alternative hypothesis with no free
parameters. Similarly, the Bayes factor is the test statistic with the highest
power (it minimizes the Type II error (false negative rate) at a fixed Type I
error) if the alternative hypothesis has multiple parameters with some prior
information. This is because the Bayes factor summarizes optimally the prior
information we have about the alternative hypothesis. It does not mean that
the estimated Type II error is correct, as the prior we assumed could be
wrong. It does however mean that we cannot find a better test statistic on
Type II error without knowing a better prior.

In this paper, we will offer an interpretation of the results in terms of
statistical mechanics, using concepts such as counting of states, uncertainty
quantification, and entropy. These connections have been explored before in
the context of Bayesian statistics jeffreypriorvolumes , here we extend it to
the frequentist statistics, and show that many of the familiar statistical
mechanics concepts can be applied to frequentist hypothesis testing.

The main goal of this paper is to develop analytic methods for evaluation of
p-value with the Bayes factor test statistic from the frequentist and Bayesian
perspectives. Our formalism enables evaluating the p-value without running
expensive simulations, and reproduces many of the commonly used expressions in
specific limits.111We do not discuss the decision theory once the p-value or
the Bayes factor have been established. In frequentist methodology, the
decision to accept or reject the null hypothesis is made with respect to some
predetermined p-value (e.g. 0.05) to guarantee frequentist coverage, but
choosing the specific value is problem and domain-dependent and its discussion
is beyond the scope of this paper.

The outline of the paper is as follows: in Section II we define the Bayesian
hypothesis testing. In Section III we introduce the statistical mechanics
approach to the Bayesian hypothesis testing. In Section IV we will develop an
analytic formalism for computing the false positive rate, both from the
Bayesian and frequentist perspective. In section V we will apply this
formalism to practical examples, notably an exoplanet transit search.

##  II Bayesian hypothesis testing

We are given some data $x$ and want to test them against competing hypotheses.
We will assume there is a single null hypothesis $\mathcal{H}_{0}$, which
assumes there is no discovery in the data, and a collection of the alternative
hypotheses $\mathcal{H}_{1}$, all predicting some new signal. For example,
when we are looking for a planet transit in the time-series data, then the
null hypothesis would be that stellar variability and noise alone are
responsible for the flux variations, while the alternative hypothesis would
predict an exoplanet transits cause some signal in the data. There are
multiple alternative hypotheses because we do not know planet’s properties,
such as its period, phase, or amplitude, so the alternative hypothesis is not
simple, and has parameters we need to vary.

The Bayesian approach to hypothesis testing is to examine the ratio of the
marginal likelihoods, i.e. the Bayes factor

$B_{01}=B_{10}^{-1}\equiv\frac{p(\hbox{\boldmath{$x$}}|\mathcal{H}_{0})}{p(\hbox{\boldmath{$x$}}|\mathcal{H}_{1})},$
(1)

where the marginal integral is

$p(\hbox{\boldmath{$x$}}|\mathcal{H}_{1})=\int
p(\hbox{\boldmath{$x$}}|\hbox{\boldmath{$z$}})p(\hbox{\boldmath{$z$}})d\hbox{\boldmath{$z$}}.$
(2)

Here $p(\hbox{\boldmath{$z$}})$ is the prior for the alternative hypothesis
parameters and $p(\hbox{\boldmath{$x$}}|\hbox{\boldmath{$z$}})$ is the data
likelihood under a general $z$. A typical situation is that the null
hypothesis corresponds to some specific values of $z$, such as $z_{1}=0$,
where $z_{1}$ is the amplitude of the signal for the alternative hypothesis.

We are interested in the posterior odds between the competing hypothesis which
follow directly from the Bayes factors by
$B_{01}\,P(\mathcal{H}_{0})/P(\mathcal{H}_{1})$, where $P(\mathcal{H}_{0})$
and $P(\mathcal{H}_{1})$ are the prior odds. We typically assume the prior
odds of the two hypotheses to be equal, in which case the Bayes factor is also
the posterior odds, and in the following, we will focus on the Bayes factor
only.

The Bayes factor is an integral over all possible alternative hypotheses, but
usually a relatively small range of parameters, where the data likelihood
peaks, dominates the integral. In this case it suffices to apply a local
integration at the peak of the posterior mass, which is often (but not always)
at the maximum a posteriori (MAP) peak, where posterior density peaks. If the
location of the highest posterior mass peak under $\mathcal{H}_{1}$ is
$\widehat{\boldsymbol{z}}$, this gives the Bayes factor as LEE1

$B_{10}=\frac{p(\hbox{\boldmath{$x$}}|\widehat{\boldsymbol{z}})}{p(\hbox{\boldmath{$x$}}|\mathcal{H}_{0})}\,p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})V_{\rm
post}(\widehat{\boldsymbol{z}}).$ (3)

The Bayes factor depends on the likelihood ratio of the optimal model
$p(\hbox{\boldmath{$x$}}|\widehat{\boldsymbol{z}})/p(\hbox{\boldmath{$x$}}|\mathcal{H}_{0})$,
prior of the MAP under $\mathcal{H}_{1}$
$p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})$, and the posterior volume
$V_{\rm post}(\widehat{\boldsymbol{z}})$, which is a result of the local
integration of the posterior ratio
$p(\boldsymbol{z}|\boldsymbol{x})/p(\widehat{\boldsymbol{z}}|\boldsymbol{x})$
around the peak. It can be approximated using the Laplace approximation as
$V_{\rm post}\approx(2\pi)^{d/2}\sqrt{\det\hbox{\boldmath{$\Sigma$}}}$, where
$\Sigma$ is the covariance matrix, given by the inverse Hessian,

$(\hbox{\boldmath{$\Sigma$}}^{-1})_{ij}=-\partial_{i}\partial_{j}\ln[p(\hbox{\boldmath{$x$}}|\hbox{\boldmath{$z$}})p(\hbox{\boldmath{$z$}})]_{z=\widehat{z}}.$
(4)

Here $\partial_{i}$ is the derivative with respect to $z_{i}$. The choice of
prior is subjective, and when the prior is known it is informative and should
be used. When the prior is not known, an option for a non-informative prior is
Jeffrey’s prior, which is given by the square root of the expectation of the
Hessian determinant (Fisher information), and is invariant under
reparametrization. However, Jeffrey’s prior also has the property that it
exactly cancels out the parameter dependence of $V_{\rm
post}(\widehat{\boldsymbol{z}})$ in the Bayes factor, making it dependent only
on the likelihood ratio.

We can generalize these concepts to non-Gaussian posteriors. As we will show,
the Laplace approximation can be poor if the posterior is very non-Gaussian,
but the local Bayes factor integration is still well-defined, in which case
Equation (3) can be viewed as the definition of $V_{\rm post}$ (for an example
of this see Equation (37)). If the prior is
$p(\hbox{\boldmath{$z$}}|\mathcal{H}_{1})\propto
1/V_{\text{post}}(\hbox{\boldmath{$z$}})$ the Bayes factor is directly
proportional to the likelihood ratio. Such prior choice may be called non-
informative, as only the likelihood is used to test the two hypotheses. This
can be viewed as a generalization of the Jeffrey’s prior to the non-Gaussian
posteriors.

###  II.1 Gaussian likelihood

We will here make some further assumptions which simplify the analysis and
will be useful in the examples we will present. Let us assume the data is a
real vector $\boldsymbol{x}\in\mathbb{R}^{n}$ whose likelihood is Gaussian:

$p(\boldsymbol{x}|\mathcal{H}_{i})=\mathcal{N}(\bf{x};\,\boldsymbol{\mu}^{(i)},\,\Sigma).$
(5)

The hypotheses $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$ differ only in their
prediction for the data mean $\mu$.

Since we are assuming the null hypothesis has no free parameters, we can
redefine
$\hbox{\boldmath{$\Sigma$}}^{-1/2}(\boldsymbol{x}-\boldsymbol{\mu}^{(0)})\xrightarrow[]{}\boldsymbol{x}$
to make the likelihood under $\mathcal{H}_{0}$ a standard Gaussian.

$\mathcal{H}_{1}$ on the other hand predicts some feature in the data:

$\boldsymbol{\mu}^{(1)}=\boldsymbol{m}(\boldsymbol{z})=z_{1}\boldsymbol{M}(\boldsymbol{z}_{>1})$
(6)

where $\boldsymbol{z}=(z_{1},\,z_{2},\,...\,z_{d})$ is parametrizing different
options in $\mathcal{H}_{1}$. We have split the amplitude of the signal
$z_{1}$ from the remaining non-linear parameters
$\boldsymbol{z}_{>1}=(z_{2},\,...\,z_{d})$. By redefining the amplitude
$z_{1}$, we can normalize the template:
$\boldsymbol{M}(\boldsymbol{z}_{>1})\cdot\boldsymbol{M}(\boldsymbol{z}_{>1})=1$.

We will assume that the prior variations are much slower than the likelihood
variations, so the optimal model $\boldsymbol{m}(\widehat{\boldsymbol{z}})$ is
the point on the model manifold, which is the closest to the data. It has to
satisfy:

$(\boldsymbol{x}-\boldsymbol{m})\cdot\partial_{i}\boldsymbol{m}=0\qquad
i=1,\,2,\,...\,d.$ (7)

The log-likelihood ratio is related to the $\chi^{2}$ improvement of the fit:

$\displaystyle E$
$\displaystyle\equiv\log\frac{p(\hbox{\boldmath{$x$}}|\widehat{\boldsymbol{z}})}{p(\hbox{\boldmath{$x$}}|\mathcal{H}_{0})}=\frac{1}{2}\widehat{z}_{1}^{2},$
(8)

where we have used that
$\boldsymbol{x}\cdot\boldsymbol{m}(\widehat{\boldsymbol{z}})=\boldsymbol{m}^{2}(\widehat{\boldsymbol{z}})$
by Equation (7). We can identify $\widehat{z}_{1}$ as the signal-to-noise-
ratio.

The posterior is in general not Gaussian and the Laplace approximation may be
a poor approximation (e.g. Section V.3). Here we adopt it anyways to gain some
intuition. Equation (4) gives the posterior covariance matrix:

$\Sigma^{-1}_{ij}=\partial_{i}\boldsymbol{m}\cdot\partial_{j}\boldsymbol{m}-(\boldsymbol{x}-\boldsymbol{m})\cdot\partial_{i}\partial_{j}\boldsymbol{m}.$
(9)

The first term is independent of the data and thus equals the Fisher
information matrix (which is defined as an average over data realizations)

$I_{ij}=\partial_{i}\boldsymbol{m}\cdot\partial_{j}\boldsymbol{m},$ (10)

while we drop the second term since $\boldsymbol{x}-\boldsymbol{m}$ is rapidly
oscillating around zero if there are no systematic residuals and the Hessian
of the model is typically varying more slowly (Gauss-Newton approximation).
Although the Bayes factor is in general a function of the data
$\boldsymbol{x}$, we have here approximated it solely as a function of the
best fit parameters $\widehat{\boldsymbol{z}}$, inferred from the data:

$B_{10}(\widehat{\boldsymbol{z}})=\frac{(2\pi)^{d/2}p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})}{[\det
I(\widehat{\boldsymbol{z}})]^{1/2}}\exp\bigg{\\{}\frac{1}{2}\widehat{z_{1}}^{2}\bigg{\\}}.$
(11)

Note that the Fisher determinant can be restricted to the
$\boldsymbol{z}_{>1}$ components only, since there are no correlations between
$z_{1}$ and the other parameters: $I_{1i}=\boldsymbol{M}\cdot
z_{1}\partial_{i}\boldsymbol{M}=z_{1}\partial_{i}\boldsymbol{M}^{2}=0$ and
$I_{11}=\boldsymbol{M}^{2}=1$. This property is exact and does not rely on the
Laplace approximation, as can be seen by using the Gaussian likelihood (5) in
the Bayes factor computation (1). We also note that the errors on the non-
linear parameters scale as the inverse signal-to-noise-ratio, since
$I_{ij}=z_{1}^{2}\,\partial_{i}\boldsymbol{M}\cdot\partial_{j}\boldsymbol{M}$,
but this is only true under the Fisher approximation to the posterior.

Many cases of practical interest fall under the assumptions in the present
section, often with some necessary preprocessing. For example, the Kepler
space telescope data are not Gaussian distributed because of the outliers, but
they can be Gaussianized GMF . We will show multiple applications of this
setup throughout the paper: linear models in Subsection III.1 and V.1,
periodogram in Subsection III.2 and V.2 and exoplanet transit search in
Section V.3.

##  III Energy-Entropy decomposition

In a continuous parameter space, the nearby models are not independent. We can
think of two models as indistinguishable if we cannot favor one or the other
after seeing the posterior. This suggests that all models within a posterior
volume can be considered one unit, which we call a state222Equivalently one
can define a state as a collection of models whose relative entropy is small,
see jeffreypriorvolumes .. This is analogous to the shift from classical
statistical mechanics to quantum statistical mechanics, where the discrete
states are counted in units of their uncertainty volume.

Let’s first consider the Jeffrey’s prior choice, which is non-informative in
the sense that it assigns an equal probability to each state
jeffreypriorvolumes , so to each indistinguishable model. The logarithm of the
Bayes factor

$F\equiv\log
B_{10}=E-\log\frac{1}{p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})V_{\rm
post}(\widehat{\boldsymbol{z}})}=E-\log N\equiv E-S$ (12)

is reminiscent of the thermodynamics relation for the free energy if we
identify $E$ with energy, since the entropy $S=\log N$ is independent of
$\widehat{\boldsymbol{z}}$ and equals the logarithm of the number of states in
the prior volume. For a more general prior, the thermodynamics relation has to
be generalized to $F=E+U-S$, where the potential energy
$e^{U}=\frac{p(\boldsymbol{z})}{p_{Jeff}(\boldsymbol{z})}$ is given by the
prior relative to the Jeffrey’s prior and the total energy is $E+U$.

Our definition of energy and entropy should be viewed from the hypothesis
testing point of view, where the energy $E+U$ is the only ”macroscopic”
parameter that influences the outcome of the test. The other parameters are
”microscopic” in the sense that we do not care about their values in the test.
Entropy is the logarithm of the number of microstates, given the macro state,
as usual. To be precise, the entropy should only count the states which give
rise to the same macro state, so the states with the same energy. Such count
corresponds to the Bayes factor which ignores the look-elsewhere effect
associated with the amplitude parameter; we will call it $B_{>1}$. It is not
justified from the Bayesian perspective as the prior would have to be fixed
after seeing the data. Still, it makes more sense from the statistical
mechanics point of view, and we will return to this quantity when it appears
in the frequentist analysis in the next section.

Note that the entropy is always positive and therefore always favors
$\mathcal{H}_{0}$, since the posterior is narrower than the prior. Energy has
to surpass entropy for the alternative hypothesis to prevail. This is the
Occam’s razor penalty, which is built into the Bayes factor. The main
deficiency of the Bayes factor is in situations where the prior of the
alternative hypothesis is poorly known, as is invariably the case when we
encounter a new phenomenon we have never seen before. In this case Bayes
factor is sensitive to the assumption on the prior of alternative hypothesis,
which can be arbitrary, as we will now illustrate with two examples.

###  III.1 Example: linear model

First consider the setup in Subsection II.1 and further assume
$\mathcal{H}_{1}$ models the data as a linear superposition of $d$ features
$\boldsymbol{m}_{i}\in\mathbb{R}^{n}$:

$\boldsymbol{m}_{\text{linear}}(\boldsymbol{w})=\sum_{i=1}^{d}w_{i}\boldsymbol{m}_{i}$
(13)

Without loss of generality, we may assume the features are orthonormal
$\boldsymbol{m}_{i}\cdot\boldsymbol{m}_{j}=\delta_{ij}$ by applying the Gram-
Schmidt algorithm if they are not.

The optimal model is a projection of the data on the model plane
$\widehat{w}_{i}=\boldsymbol{x}\cdot\boldsymbol{m}_{i}$ and has the log-
likelihood ratio $E=\frac{1}{2}\sum_{i=1}^{d}\widehat{w}_{i}^{2}$. The
posterior is Gaussian, so the Laplace approximation is exact and gives the
posterior volume $(2\pi)^{d/2}$.

Often we do not have any prior information, and it makes sense to adopt the
Jeffrey’s prior, which is uniform for the linear parameters. For simplicity,
we assume the prior volume to be a ball with
$\sum_{i=1}^{d}w_{i}^{2}<w_{max}^{2}$, but we do not know how to choose
$w_{max}$. We cannot just send it to infinity since the Bayes factor

$B_{10}=\frac{(2\pi)^{d/2}e^{E}}{V(B^{d})w_{max}^{d}}$ (14)

would go to zero. Here, $V(B^{d})$ is the volume of the $d$ dimensional unit
ball. We will show another, more realistic, example with the same problem.

###  III.2 Example: periodogram

In this example, we are given $n$ time series measurements $x_{i}=x(t_{i})$
with uncertainties $\sigma$. We are searching for harmonic periodic signals
Lomb1976 ; Scargle1982 ; VanderPlas2018 :

$m_{i}(\hbox{\boldmath{$z$}})=A\sin(\omega t_{i}+\phi)$ (15)

Here $A$ is the signal’s amplitude, $\omega$ is the unknown frequency, and
$\phi$ is the phase.

To make the calculations tractable we assume dense data sampling over many
periods of the signal, so that we can replace the discrete sum in the scalar
products of Subsection II.1 with the integral
$\sum_{i=1}^{n}\xrightarrow[]{}\frac{n}{T}\int_{0}^{T}dt$, where $T$ is the
total time of the measurements. We will work with the amplitude $z_{1}$
relative to the normalized template, which has $A=\sigma\sqrt{2/n}$. We assume
a uniform prior on all three parameters, so

$p(z_{1},\omega,\phi)=\frac{1}{z_{max}\,\omega_{max}\,2\pi},$ (16)

where $z_{max}$ is some arbitrary large cutoff on the amplitude and
$\omega_{\rm max}$ is set by some physical properties of the signal or by the
experimental limitations, e.g. by the Nyquist frequency.

The inverse of the Fisher information matrix (10) is

$I^{-1}=\begin{bmatrix}\sigma_{\omega}^{2}&\sigma_{\omega\phi}\\\
\sigma_{\omega\phi}&\sigma_{\phi}^{2}\end{bmatrix}=E^{-1}\begin{bmatrix}6/T^{2}&-3/T\\\
-3/T&2\end{bmatrix}.$ (17)

Interestingly, the frequency-phase correlation is relatively high
$\sigma_{\omega\phi}/\sigma_{\omega}\sigma_{\phi}=-\sqrt{3}/2\approx-0.87$.
The Bayes factor

$B_{10}=\frac{\sqrt{6\pi}\,e^{E}}{z_{max}\,\omega_{max}T\,E},$ (18)

again suffers from the unknown amplitude cutoff. If we set it too high, we
might reject a true periodic signal.

To overcome the arbitrariness of the prior choice, one could try to learn the
prior from the data under some hyper-parametrization. This leads to a
hierarchical or empirical Bayesian analysis hierarchicalBayes ; empiricalBayes
, but this is not possible if we have no data that can inform the
hyperparameters.

Note that the constant energy Bayes factor $B_{>1}$ is independent of the
unknown cutoff in both examples, but we have not yet justified its use. We
will see its appearance when we phrase the test in terms of the frequentist
false positive rate, which we discuss now.

##  IV Frequentist hypothesis testing

In frequentist hypothesis testing, we first define a test statistic, which
should be chosen to maximize the contrast against the alternative hypothesis.
As argued in the previous section, the optimal statistic between two well-
specified hypotheses are the posterior odds. Any monotonically increasing
function of the posterior odds is an equally good test statistic in the
frequentist sense, so for convenience, we will work with $F=\log B_{10}$. If
we adopt the generalized Jeffrey’s prior, then
$p(\hbox{\boldmath{$z$}}|\mathcal{H}_{1})V_{\text{post}}(\hbox{\boldmath{$z$}})$
is a constant and the test statistic becomes the likelihood ratio. This is the
most common frequentist test statistics, but it is not prior independent.
Instead, it corresponds to a specific prior choice - the generalized Jeffrey’s
prior. Prior information about some parameters can reduce the Type II error at
a constant Type I error by the Neyman-Pearson argument, so in general, we will
work with $F$ as the optimal test statistic.

In frequentist methodology, we quantify a test statistic using its Type I
false positive rate (p-value) under the null hypothesis $\mathcal{H}_{0}$.
While $\mathcal{H}_{1}$ is used to define the test statistic, the distribution
of the test statistics under $\mathcal{H}_{0}$ is independent of whether
$\mathcal{H}_{1}$ distribution is true or not. If this number is sufficiently
small, the null hypothesis is rejected. Unlike the posterior odds, the p-value
is independent of the correctness of the prior of the alternative hypothesis,
such as the value of $w_{\rm max}$: it only depends on the null hypothesis
itself. This is a great advantage of p-value over the posterior odds, and is
the main reason why posterior odds have not been widely adopted for hypothesis
testing even in Bayesian textbooks BayesFactorOverview .

While the posterior odds can be derived entirely from the data using the
likelihood (the likelihood principle), the p-value is given by the frequency
distribution of the test statistic under the null hypothesis, which generally
requires simulating the null hypothesis many times to obtain it. Since this
requires evaluating the test statistic for each simulation, and since we argue
the test statistic itself should be the posterior odds, this can be
significantly more expensive in frequentist methodology than evaluating
posterior odds once on the data, as done in Bayesian methodology. For this
reason, we will study analytic techniques for evaluating the false-positive
rate.

Given the choice of the test statistic $F$, the frequentist analysis proceeds
to quantify the frequency of null events above its observed value
$F(\boldsymbol{x})$ in a frequentist sense, which means under a long series of
trials. This gives the p-value as a function of the observed test statistic
$F(\boldsymbol{x})$, which quantifies the false positive rate (FPR) as
$P(F>F(\boldsymbol{x}))$.

###  IV.1 p-value from frequentist property of Bayes factor

We will assume that the Bayes factor can be expressed as a function of the
optimal inferred parameters $\widehat{\boldsymbol{z}}$, as in Equation (11).
Then the distribution of $F$ can be inferred from the distribution of the
optimal parameters. For the null hypothesis, the majority of events will have
optimal parameters close to the null hypothesis values, but there will be
outliers whose frequency we would like to evaluate. If the prior is correct,
the Bayes factor predicts the relative frequency $\mathcal{P}$ of the outcomes
of the two hypotheses in a long series of trials:

$\mathcal{P}(\widehat{\boldsymbol{z}}|\mathcal{H}_{0})=\mathcal{P}(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})B_{01}(\widehat{\boldsymbol{z}})\approx
p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})e^{-F(\widehat{\boldsymbol{z}})}$
(19)

The frequency of $\mathcal{H}_{1}$ events follows the prior: each sample with
inferred MAP parameters $\widehat{\boldsymbol{z}}$ will have a true value of
parameters $z$ approximately within $V_{\rm post}$ of
$\widehat{\boldsymbol{z}}$. Thus, as long as the prior is sufficiently smooth
relative to the posterior, we have
$\mathcal{P}(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})\approx
p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})$.

The p-value is found by integrating over the parameters which yield the
desired test statistic:

$\displaystyle P(F>F(\boldsymbol{x}))$
$\displaystyle=\int_{F(\widehat{\boldsymbol{z}})>F(\boldsymbol{x})}p(\widehat{\boldsymbol{z}}|\mathcal{H}_{1})e^{-F(\widehat{\boldsymbol{z}})}d\widehat{\boldsymbol{z}}=$
(20) $\displaystyle=\int_{F(\boldsymbol{x})}^{\infty}p(F)e^{-F}dF,$ (21)

where $p(F)$ is the marginal prior density. This result is saying that the
probability density of finding some $F$ under the null hypothesis is given by
the probability density of finding it under the alternative hypothesis prior,
but exponentially damped with $F$. Note that Equation (20) is equivalent to

$P(F>F(\boldsymbol{x}))=\int_{F(\widehat{\boldsymbol{z}})>F(\boldsymbol{x})}\frac{e^{-E(\widehat{\boldsymbol{z}})}}{V_{\text{post}}(\widehat{\boldsymbol{z}})}d\widehat{\boldsymbol{z}},$
(22)

so the region of integration depends on the prior, but the integrand does not.
This means that the prior selects the parameter range where the false
positives can be generated, but the rate at which they are generated at those
parameters is prior-independent.

These expressions are useful when we have an analytic expression for the
posterior volume and can perform the integral analytically. We will show
examples of this situation in Section V. Here we will derive a useful
approximation which relates the p-value and the constant energy Bayes factor
$B_{>1}=e^{E}p_{>1}(\widehat{\boldsymbol{z}})V_{\text{post}}^{(>1)}(\widehat{\boldsymbol{z}})$.
This will enable the p-value calculation even when the Bayes factor is only
available numerically. Note that the variations of $V_{\text{post}}$ as a
function of $E$ are much slower than the suppression $e^{-E}$. One can
therefore consider the posterior volume to be a constant evaluated at
$z_{1}=\sqrt{2E(\boldsymbol{x})}$ when taking the $z_{1}$ integral. This gives

$P(F>F(\boldsymbol{x}))\approx\frac{e^{-E(\boldsymbol{x})}}{\sqrt{4\pi
E(\boldsymbol{x})}\,p_{>1}(\widehat{\boldsymbol{z}})V_{\text{post}}^{(>1)}(\widehat{\boldsymbol{z}})}=\frac{1}{B_{>1}\sqrt{4\pi
E}}$ (23)

which is approximately valid for small p-values. This expression implies that
the p-value can be directly inferred from the constant energy Bayes factor
$B_{>1}$ with no need to do any additional integrals.

###  IV.2 Frequentist derivation of p-value

Figure 1: Illustration of a two-dimensional data space ($n=2$). The concentric
circles are the likelihood contours of $\mathcal{H}_{0}$. The solid line is
the $\mathcal{H}_{1}$ manifold ($d=1$). If the data point $\boldsymbol{x}$
lies on the dashed line, the model $\boldsymbol{m}(\boldsymbol{z})$ satisfies
Equation (7) and is at least locally the optimal model. Integral of Equation
(25) runs over the dashed line.

We will here present an approach that is complementary to the previous section
in the sense that it does not use any of the Bayesian concepts but still
reproduces the same results. Since there is no concept of a prior, we will be
working with the likelihood ratio as a test statistic (which is the Bayes
factor with the Jeffrey’s prior). We will assume the setup as in Subsection
II.1.

First, let us determine how likely would we have seen some model
$\boldsymbol{m}(\boldsymbol{z})$ as the optimal model if the null hypothesis
was true. We add the probabilities over the suitable region of the data space:

$p(\boldsymbol{m}(\boldsymbol{z})|\mathcal{H}_{0})=\int_{\mathcal{N}}p(\boldsymbol{x}|\mathcal{H}_{0})d\boldsymbol{x},$
(24)

where $\mathcal{N}$ is the set of all points which have
$\boldsymbol{m}(\boldsymbol{z})$ as their optimal model. $\mathcal{N}$ is a
subset of the plane defined by Equation (7).This equation defines an affine
$n-d$ dimensional plane which is an orthogonal complement to the tangent plane
to the $\mathcal{H}_{1}$ manifold at the point
$\boldsymbol{m}(\boldsymbol{z})$ (for illustration, see Figure 1). However,
$\mathcal{N}$ is not the entire plane since Equation (7) is a necessary but
not sufficient condition. In other words, the posterior may be multi-modal. To
make progress, we will have to neglect this and assume $\mathcal{N}$ to be the
entire plane. By doing so, our final result for the p-value is in fact an
upper bound, which becomes more and more accurate the lower the p-value.

The model $\boldsymbol{m}(\boldsymbol{z})$ is both a point on the plane
$\mathcal{N}$ and a normal to the plane, as can be seen from Equation (7) with
$i=1$. Writing
$\boldsymbol{x}=\boldsymbol{m}(\boldsymbol{z})+\boldsymbol{x}_{\bot}$ for
points on the plane and using orthogonality, we can easily evaluate the
Gaussian integrals (24):

$p(\boldsymbol{m}(\boldsymbol{z})|\mathcal{H}_{0})=\int
p(\boldsymbol{x}|\mathcal{H}_{0})d\boldsymbol{x}_{\bot}=\frac{e^{-\boldsymbol{m}^{2}/2}}{(2\pi)^{d/2}}.$
(25)

This result is saying that the probability density in the data space of a
model occurring as an optimal model under $\mathcal{H}_{0}$ depends only on
$\boldsymbol{m}^{2}/2=E$. Therefore, the probability of observing $E$ under
the null hypothesis equals the probability of finding any $\boldsymbol{m}$
with $\boldsymbol{m}^{2}/2=E$ as the optimal model. This gives

$p(E|\mathcal{H}_{0})=\frac{1}{\sqrt{2E}}\frac{e^{-E}}{(2\pi)^{d/2}}V_{\text{shell}},$
(26)

where $V_{\text{shell}}$ is the data space volume of models with
$\boldsymbol{m}^{2}/2=E$, in other words, the constant likelihood surface
volume. We have picked up a factor $(2E)^{-1/2}$ when transforming the density
with respect to $\widehat{z_{1}}$ to the density with respect to $E$. The
volume of the constant likelihood shell in the data space is an integral over
the $\mathcal{H}_{1}$ manifold at fixed $z_{1}=\sqrt{2E}$. It can be evaluated
in the $\boldsymbol{z}_{>1}$ parameter space by integrating the square root of
the determinant of the metric over the coordinate range:

$V_{\text{shell}}=\int\,d\boldsymbol{m}(\sqrt{2E},\boldsymbol{z}_{>1})=\int\det
I_{ij}(\sqrt{2E},\boldsymbol{z}_{>1})^{1/2}d\boldsymbol{z}_{>1},$ (27)

where the metric on the $\mathcal{H}_{1}$ manifold coincides with the Fisher
information matrix (10). Combining Equations (26) and (27) we find

$p(E|\mathcal{H}_{0})=\frac{e^{-E}}{\sqrt{2E}}\int\frac{d\boldsymbol{z}_{>1}}{V_{\text{post}}(\sqrt{2E},\boldsymbol{z}_{>1})},$
(28)

where $V_{\text{post}}=(2\pi)^{d/2}\det I_{ij}^{-1/2}$.

Even though we use frequentist formalism, the posterior volume enters as the
determinant of the transformation from the data space to the parameter space.
We reproduce the result of Equation (22).

There are some exciting applications where the derived expressions are exact
(see Section V.1) since the posterior has a single mode. However, in general,
the frequentist derivation made it manifest that in the presence of multiple
modes, the p-value results are only an upper bound, which becomes useless if
the p-value is high. We now extend our results to all p-values.

###  IV.3 Non-asymptotic p-value

We have found an analytic expression for the false positive rate, both from
the frequentist and Bayesian perspectives. However, it is only valid when the
false positive rate is low. We will show that this result can be easily
extended to all p-values if the posterior volume is much smaller than the
prior volume, as was done in LEE1 .

Let’s partition the alternative hypothesis manifold in $K$ smaller manifolds
and consider searches over the smaller manifolds. Let’s choose the partition
in a way that FPR in all the small searches is the same and call it
$P_{\text{small}}(F>F(\boldsymbol{x}))$.

Suppose the posterior volume is much smaller than the prior volume. In that
case, most data realizations will have their posterior volume well within the
prior range of some smaller manifold, even when $K$ is relatively large.
Therefore, the probability of not finding a false positive in the original
search equals the probability of not finding it in any of the small searches:

$1-P(F>F(\boldsymbol{x}))=\big{\\{}1-P_{\text{small}}(F>F(\boldsymbol{x}))\big{\\}}^{K},$
(29)

If $K$ is large, $P_{\text{small}}(F>F(\boldsymbol{x}))$ becomes small, and we
can compute it using the asymptotic expression of Equation (20) or (23):
$P_{\text{small}}(F>F(\boldsymbol{x}))=P_{\text{asym}}(F>F(\boldsymbol{x}))/K$.
Taking the large $K$ limit we find

$P(F>F(\boldsymbol{x}))=1-\exp\\{-P_{\text{asym}}(F>F(\boldsymbol{x}))\\},$
(30)

This is a continuous parameters generalization of Sidák correction, which
itself is a generalization of Bonferroni correction, which are commonly used
for discrete states where the trials factor is a well-defined concept and
referred to as multiple test comparison.

###  IV.4 Statistical mechanics interpretation

Figure 2: Illustration of the statistical mechanics interpretation of the
p-value. We emphasize the difference between the amplitude and the other
parameters. Left panel: parameter space of a linear model (13) with two
degrees of freedom. This is equivalent to the periodogram model (15) at a
fixed frequency $\omega$. Constant energy (likelihood) shells are circles of
radius $z_{1}=(w_{1}^{2}+w_{2}^{2})^{1/2}$. We assume a uniform prior with
equal uncertainty on both variables and show posterior volume as the area of
each mesh cell. The expected ratio of the null hypothesis to the alternative
hypothesis events in each cell is given by the Boltzmann factor in Equation
(31). It dies off exponentially as a function of expected energy. The color
intensity in the plot is proportional to this Boltzmann factor, which is shown
along one radial direction as a blue line. False positive rate counts all the
states in the region exterior to the circle of observed $F(\boldsymbol{x})$
and weights them with the Boltzmann factor (31). Only the region close to the
observed test statistic circle will contribute due to the exponential
suppression. Increasing the prior on $w_{1}$ and $w_{2}$ increases the total
number of states and reduces the Bayes factor but has no impact on the p-value
(32). Right panel: now, we also vary the frequency, searching for sinusoidal
signal (15) over all frequencies within some prior range. The observed
$F(\boldsymbol{x})$ surface is a cylinder and the false positive rate is
proportional to its surface area. Increasing the prior range of $\omega$ over
which we search for sinusoidal signal reduces the Bayes factor because
additional states were introduced (14). Contrary to the first example, the
false positive rate has now also increased (33) because some new states are
close to the observed $F(\boldsymbol{x})$ shell.

Returning to the discrete picture of the continuous parameter space presented
in Section III we recognize the integral
$\int\frac{d\boldsymbol{z}}{V_{\text{post}}(\boldsymbol{z})}$ in Equation (22)
as the the continuous version of the sum over states. The p-value is therefore
the sum over all states which generate false positives, each weighted with the
Boltzmann factor $e^{-E}$. The p-value is therefore the partition function of
the canonical ensemble with $k_{B}T=1$:

$P(F>F(\boldsymbol{x}))=\sum_{F(\text{state})>F(\boldsymbol{x})}e^{-E(\text{state})}=e^{-F_{f}},$
(31)

where $F_{f}$ is the frequentist free energy. It is evident why the p-value is
not sensitive to the high amplitude cutoff: high energy states are suppressed
by the Boltzmann factor and do not contribute to the p-value. An illustration
of this interpretation is presented in Figure 2.

Note that the p-value approximation (23) implies that the Bayesian free energy
$\log B_{>1}$ and the frequentist free energy $-\log P(F>F(\boldsymbol{x}))$
are one and the same thing, up to logarithmic corrections in energy. Note that
in physics, the thermodynamic and statistical mechanic free energies are also
the same in the thermodynamic limit.

##  V Results

Figure 3: We show an example of a typical noise simulation with the injected
planet with a period $P=100$ days and transit duration $\tau=0.4$ days. We
show the likelihood ratio in the neighborhood of the MAP parameters
$\widehat{\boldsymbol{z}}$ (upper panels). Laplace approximation expands the
likelihood ratio’s logarithm to a quadratic order and approximates the peak
with a Gaussian (lower panels). Note that the Laplace approximation is poor,
corresponding to an inaccurate p-value estimation (Figure 4).

In this section, we apply the developed formalism to three examples with
increasing complexity. We start with the linear model and periodogram and
compare our results with the literature. We then turn to the more realistic
analysis of the exoplanet transit search.

###  V.1 Linear model

We start with the linear model as in Section III. As seen from Equation (14)
$F$ and $E$ only differ by a constant, so they are equally good test
statistics, and we might as well work with $E$. The posterior volume is a
constant, so the integral of Equation (22) just picks up the volume of the
constant energy surface, which is a $d-1$ dimensional sphere. We get

$P(E>E(\boldsymbol{x}))=\int_{E(\boldsymbol{x})}^{\infty}\frac{E^{d/2-1}\,e^{-E}}{\Gamma(d/2)}dE,$
(32)

where we have used that the volume of a $d-1$ dimensional sphere is
$2\pi^{d/2}/\Gamma(d/2)$, with $\Gamma$ the Gamma function.

The resulting cumulative distribution function is a $\chi^{2}$-distribution
with $d$ degrees of freedom, and we reproduce the well-known result that
$\chi^{2}$ of a linear model with $d$ features is distributed as a $\chi^{2}$
distribution with $d$ degrees of freedom wilksTheorem . The p-value is
increasing with $d$ at a constant $\chi^{2}=2E$, which is a reflection of the
entropy versus energy competition: there are more states on the shell of a
constant energy $E$ if $d$ is higher.

Note that the p-value is independent of the unknown amplitude cutoff
$w_{max}$.

###  V.2 Periodogram

We now turn to the periodogram example, which is defined in Section III.2. We
have seen in Equation (18) that the Bayes factor and the likelihood ratio are
a simple function of one another, so their null distributions are related by
the change of variable formula
$p(F|\mathcal{H}_{0})=|dF/dE|\,p(E|\mathcal{H}_{0})=|1-1/E|\,p(E|\mathcal{H}_{0})$.

Once again, the posterior volume is independent of the $\boldsymbol{z}_{>1}$
parameters, and the integral (22) picks up the volume of the constant energy
surface, which is, in this case, a cylinder, so $2\pi\omega_{max}$. We get

$p(E|\mathcal{H}_{0})=\frac{1}{B_{>1}\,\sqrt{4\pi E}},$ (33)

so the false positive rate is closely related to the Bayes factor, as in the
approximation (23). Note that while the amplitude cutoff cancels in the
p-value, the frequency cutoff does not. This is the look-elsewhere effect: the
larger the frequency space search, the larger the false positive rate at a
fixed likelihood-ratio.

These results agree with expressions of Baluev_2008 , which use the formalism
of Davies ; Davies2 . This mathematical formalism is based on extremes of
random processes of various distributions such as gamma ($\chi^{2}$)
distribution, and needs to be derived separately for each distribution. In
addition, this formalism formally only gives an analytic lower limit to the
corresponding extreme value distributions, while in practice, equality is
assumed without proper justification. Our formalism provides a different
derivation that results in the same expressions in the periodogram case.

###  V.3 Exoplanet search

Figure 4: We simulate 300 noise series and search for the planetary signal. In
each realization we find a candidate with highest Bayes factor and evaluate
it’s p-value using Equations (23) and (30). In the left panel, we then compare
it with a fraction of noise realizations where a more significant planet
candidate was found. A good agreement between the empirical and analytical
p-value is found, regardless of the prior choice. Dotted lines are the
corresponding results using the asymptotic limit of Laplace approximation,
which is inadequate in this case. In the right panel, we show the likelihood-
ratio test statistic distribution for the null hypothesis using the same
simulations. Note a large trials factor which causes a noise-only simulation
to produce a ${\rm SNR}=\sqrt{2E}>6$ in most realizations. There is an order
of magnitude difference in the trials factor for the circular prior versus the
wide or realistic prior.

As a more complex non-trivial example of the formalism we developed, we
consider exoplanet detections in the transit data, where the planet orbiting
the star dims the star when it transits across its surface. We have a time
series of star’s flux measurements $x_{i}=x(t_{i})$. In the absence of a
transiting planet, the data is described by a stationary correlated Gaussian
noise modeling stellar variability333Long-term trends in the Kepler data are
removed by the preprocessing module kepler_process . Outliers are a source of
non-Gaussianities, but they can be efficiently Gaussianized without affecting
the planet transits GMF . Here we ignore other defects in the data, such as
binary stars, sudden-pixel-sensitivity-dropouts, etc.. Such data are for
example collected by the Kepler Space Telescope KeplerSpaceTelescope  and
Transiting Exoplanet Survey Satellite (TESS) TESS .

One would like to compare the hypothesis $\mathcal{H}_{1}$ that we have a
planet in the data to the null hypothesis that there is only noise. As argued
in this paper, we will use the Bayes factor as a test statistic to incorporate
informative prior information and non-Gaussian posteriors. The significance of
a discovery is then reported as an associated false positive probability of
exceeding the observed value of the Bayes factor. We will first outline the
procedure for calculating the Bayes factor given the prior for the transit
model parameters. Then, using simulations of the null hypothesis, we will
demonstrate that Equation 30 gives reliable results for the p-value of the
Bayes factor (Subsection V.3.2). In Subsection V.3.3, we will discuss several
prominent prior choices and demonstrate how a realistic prior choice can
reduce Type II error at a constant Type I error.

####  V.3.1 Bayes factor

Planet transit model $m(t)$ can be parametrized by $d=4$ parameters: transit
amplitude, period, phase, and transit duration,
$\hbox{\boldmath{$z$}}=(A,P,\phi,\tau)$. The transit model is of the form

$m_{i}(\hbox{\boldmath{$z$}})=A\sum_{m=1}^{M}U\bigg{(}\frac{t_{i}-(m+\phi)P}{\tau}\bigg{)}.$
(34)

It is a periodic train of $M$ transits with period $P$, phase $\phi$,
amplitude $A$, and transit duration $\tau$. $U(x)$ is a U-shaped transit
template that is nonzero in the region (-1/2, 1/2) and depends on the limb
darkening of the stellar surface simple_limb_darkening .

The likelihood is Gaussian and stationarity ensures that the Fourier
transformation $\mathcal{F}$ diagonalizes the covariance matrix with the power
spectrum on the diagonal. The power spectrum
$\mathcal{P}(\omega)=\langle|\mathcal{F}\\{n\\}|^{2}\rangle$ can be learned
from the data FourierGP .

In the language of Subsection II.1 we introduce

$\boldsymbol{x}\xrightarrow[]{}\frac{\mathcal{F}\\{\boldsymbol{x}\\}}{\mathcal{P}^{1/2}}\qquad\boldsymbol{m}\xrightarrow[]{}\frac{\mathcal{F}\\{\boldsymbol{m}\\}}{\mathcal{P}^{1/2}}$
(35)

to make the likelihood a standard Gaussian. We also rescale the amplitude $A$
to make the model normalized.

It will be convenient to first optimize over the linear parameter $z_{1}$ and
then over the remaining non-linear parameters. Amplitude’s optimal value at
fixed $\boldsymbol{z}_{>1}$ is by the Equation (7)

$\widehat{z}_{1}(\boldsymbol{z}_{>1})=\frac{\int\mathcal{F}\\{x\\}^{*}\mathcal{F}\\{m\\}d\omega/\mathcal{P}}{\sqrt{\int|\mathcal{F}\\{m\\}|^{2}d\omega/\mathcal{P}}}.$
(36)

where we have restored a model normalization factor. This is the matched
filtering expression for the signal-to-noise ratio GMF  and can be computed
efficiently using Fourier transforms GMF  for all phases $\phi$ on a fine grid
at once. Our strategy for finding the MAP solution is therefore to first find
a maximum likelihood estimator (MLE) by scanning over the entire parameter
space (for details see GMF ) and then use it as an initial guess in a
nonlinear MAP optimization.

We calculate the posterior volume associated with the most promising peaks by
marginalizing over the planet’s parameters in the vicinity of the peak:

$V_{\text{post}}^{>1}=\int\frac{p(\boldsymbol{z}_{>1})\,e^{\widehat{z}_{1}(\boldsymbol{z}_{>1})^{2}/2}}{p(\widehat{\boldsymbol{z}})\,e^{\widehat{z}_{1}^{2}/2}}\,d\boldsymbol{z}_{>1},$
(37)

One would be tempted to employ the Laplace approximation. Figures 3 and 4 show
this does not give satisfactory results: we are not in the asymptotic limit.
While a frequentist approach using Wilks’ theorem would become invalid, here,
we can compute the full marginal integral of Bayesian evidence. Integral (37)
is only three-dimensional, so we take the Hessian at the peak to define the
Gaussian quadrature integration scheme gaussquadratureformulas  of degree 7,
implemented in quadpy , which requires 24 integrand evaluations.

####  V.3.2 p-value

To obtain the p-value, we use Equation 23 and extend it to all p-values using
(30).

We test our analytical expression for the p-value of the Bayes factor by
simulating the null hypothesis and comparing the computed p-value with the
empirically determined value. We evaluate 300 realizations of the null
hypothesis. We take a realistic power spectrum extracted from Kepler’s data
for the star Kepler 90. The power spectrum and example realizations are shown
in Figures 2 and 4 of GMF . A realization is a uniformly distributed choice of
the phases of the Fourier components and normally distributed amplitudes of
the Fourier components with zero mean and variance given by the power spectrum
(an example realization is shown in Figure 4 of GMF ). In each of the
resulting time-series, we then determine the Bayes factor of the planet
hypothesis against the null hypothesis (Subsection V.3.1) and its p-value.

Analytic and empirical p-value are compared in Figure 4, achieving an
excellent agreement. This shows that our formalism enables the evaluation of
p-value beyond the asymptotic limit, thus generalizing Wilks’ theorem.

####  V.3.3 The choice of prior

Figure 5: We show the potential $U=\log
p(\boldsymbol{z})/p_{Jeff}(\boldsymbol{z})$ for different choices of the
$\tau$ prior. This quantity measures to what extent is the prior informative
with the Jeffrey’s prior being non-informative. This quantity directly impacts
the hypothesis test’s outcome, see Section III. Fixing $\tau$ to the Kepler
value is equivalent to a constant potential on the red line and $-\infty$
everywhere else. A broad Jeffrey’s prior is equivalent to a constant potential
over the entire search domain. The informative prior is shown with the color
scale. It peaks close to the Kepler’s value, but is broadened to account for
variation in planet eccentricity orbit orientation and stellar density
uncertainty. Figure 6: We show marginal probability density distributions over
the period and transit duration of the five highest $F_{f}$ candidates from
each of the 300 noise-only simulations. The distributions follow Jeffrey’s
prior, which is approximated reasonably well with the scaling estimates from
Equations (41) and (42). This confirms that regions with smaller $V_{\rm
post}$ have more false positives. The blue plot is showing a situation where
the transit duration is a free parameter in the search, yielding $p(P)\sim{\rm
const}$ and $p(\tau)\sim\tau^{-2}$. A simple analytical estimate does not give
an accurate Jeffrey’s prior for the transit duration because of the non-white
power spectrum. Therefore, we complement this analysis with a numeric
evaluation of the expected posterior volume (solid line) and show that this
gives a better match. The orange plot shows a situation where the transit
duration is fixed to the Kepler’s value. It is well described by the scaling
estimate $p(P)\propto P^{-1/3}$.

Here we discuss a prior choice for the period, phase, and transit duration. We
will see that using informative prior when available can significantly improve
detection efficiency.

We scan over the period range from 3 to 300 days and at each period over all
phases from 0 to 1. As we argued in this paper, when we do not know the prior,
it is simplest to adopt Jeffrey’s choice. We will assume this for period and
phase parameters. Note that for the phase, Jeffrey’s prior is flat from 0 to
1, which we can also view as a known (informative) prior, since the phase
cannot affect the exoplanet detectability. There is a natural value for the
transit duration parameter $\tau$, given the planet’s period and assuming the
circular, non-inclined orbit. We illustrate the impact of the prior choice by
considering three scenarios:
1) A fixed value of $\tau$, defined by assuming a circular, perfectly aligned
orbit. The transit duration is fixed by the Kepler’s third law:
$\tau_{K}(P)=q\,P^{1/3}$, with the proportionality constant given by the
star’s density, $q=(3P/(\pi^{2}G\rho_{*}))^{1/3}$. This choice can be too
restrictive for non-circular or inclined orbits and can penalize real planets,
as shown in Figure 7. Under this prior choice, the Bayes factor is
proportional to the likelihood ratio.
2) Jeffrey’s prior on a broad domain
$\tau\in[\tau_{K}(P_{min}),\,2\,\tau_{K}(P_{max})]$. This choice may include
physically implausible transit times leading to large multiplicity penalty.
3) realistic prior distribution, taking into account inclined and eccentric
orbits. Orbits are assumed to be isotropic, with eccentricities drawn from a
beta function with parameters that match the observed planets in the Kepler’s
data EccentricityPriorKipping . In addition, the star’s density is measured
with an uncertainty on the order of 15 %, causing uncertainty in $\tau_{K}$.
The transit duration prior is obtained by marginalizing over the orbit
inclination, eccentricity and star’s parameters. The distribution is a
broadened version of the delta function at $\tau_{K}$.

All three choices are visualized in Figure 5. In appendix A we derive the
approximate Jeffrey’s prior by deriving the scaling of the Fisher information
matrix. In the first two scenarios, we respectively get $p\propto\tau^{-2}$
and $p\propto P^{-1/3}$. One can take the template and calculate the Fisher
information matrix numerically for more accurate results. In Figure 6 we show
that peaks of the null hypothesis are distributed by Jeffrey’s prior, in
agreement with our predictions.

Taking $\tau$ as a completely unconstrained parameter will in general lead to
detecting planets that are not physically plausible and will therefore
increase the false positive rate and force us to reject more real planets than
necessary. In the present example shown in figure 4 this effect is not very
large if the prior is chosen such that it still introduces a reasonable cutoff
on the transit duration. But choosing such a cutoff is a choice that must be
based on physical arguments: the prior plays an essential role.

A prior that is too narrow (case 1) is also suboptimal because it will reject
some real planets. Using a template with duration $\tau_{K}$ when in fact a
planet has transit duration $\tau\neq\tau_{K}$ will reduce planet’s $SNR$ by

$\frac{{\rm SNR(detected)}}{{\rm
SNR(true)}}=\frac{\int\mathcal{F}\\{m_{\tau}\\}^{*}\mathcal{F}\\{m_{\tau_{K}}\\}d\omega/\mathcal{P}}{\int\mathcal{F}\\{m_{\tau_{K}}\\}^{*}\mathcal{F}\\{m_{\tau_{K}}\\}d\omega/\mathcal{P}},$
(38)

where this is an expected reduction over the noise realizations.

In Figure 7 we show the ROC curves (Receiver Operated Characteristic) for all
three prior choices, that is, a true positive rate as a function of the false
positive rate, both parametrized by the detection threshold. False positive
probabilities are taken from Figure 4. True positive probability is a
probability that the signal with some ${\rm SNR(true)}$ is detected above the
threshold. The detected ${\rm SNR}$ can be approximated as a Gaussian
distributed variable with unit variance. Its mean is $\mu={\rm SNR(true)}$
with the realistic and wide priors. It is additionally reduced with the
circular orbit prior because the search template differs from the true
template (Equation (38)). Using a realistic prior improves the true detection
probability at a fixed false positive probability relative to a prior that is
too narrow (circular orbit prior) because it improves the fit. It also
improves ROC relative to a prior that is too broad because it does not include
the templates that rarely happen in the data and would lead to a larger
multiplicity and hence a larger false positive rate. This figure thus
demonstrates that hypothesis testing with a realistic prior Bayes factor gives
optimal ROC (power versus p-value).

Figure 7: ROC curve (see Subsection V.3.3) defined as true positive
probability, i.e. power (1-Type II error) versus false positive probability,
i.e. Type I error (p-value), for various prior choices. We search for a
hypothetical small planet with injected radius $R=1.45R_{\bigoplus}$,
corresponding to ${\rm SNR}(true)=6.25$ (solid lines), and
$R=1.8R_{\bigoplus}$ and ${\rm SNR}(true)=7$ (dotted lines). Using the
realistic prior results in an improved ROC curve (higher true positive
probability at a given false positive probability) implying that the Bayes
factor with a realistic prior is the test statistic with the highest
statistical power. The true positive probability at a fixed p-value is
decreased for the circular orbit prior because the prior does not include all
the allowed signal template forms. On the other hand, a wide prior includes
too many signal template forms that include signal templates that occur rarely
or never, leading to a larger multiplicity and a larger false positive rate at
a constant true positive probability.

##  VI Discussion

This paper compares frequentist and Bayesian significance testing between
hypotheses of different dimensionality, where the null hypothesis is a well-
defined accepted model for the reality of the data, while the alternative
hypothesis tries to replace the null hypothesis. Both Bayesian and frequentist
methodologies have advantages and disadvantages in the setting where the
alternative hypothesis has not been observed with sufficient frequency to
develop a reliable prior. We argue that for optimal significance testing, the
Bayes factor between the two hypotheses should be used as the test statistic
with the highest power. However, the Bayes factor should not be used to
quantify the test significance when the prior of the alternative hypothesis is
poorly known. Instead, the frequentist false positive rate of null hypothesis
(Type I error or p-value) can be used, which only depends on the properties of
the null hypothesis, which is assumed to be well understood. The sensitivity
of the Bayes factor to the choice of prior is known in the context of
Lindley’s paradox Lindley57 . While there is no actual paradox, it highlights
the dependence of the Bayes factor on the choice of prior, which is
undesirable since we often do not know it. Our solution to this paradox is to
relate the Bayes factor to the p-value, which is independent of the
alternative hypothesis, as it only tests the distribution of the null
hypothesis. In this way one can use the p-value for hypothesis testing even
when using Bayesian methods, by using the Bayes factor as a test statistic.

While it is common to use the likelihood ratio as a test statistic, this is
not prior independent, but corresponds to the Jeffrey’s prior. If some prior
information is available, as in our exoplanet example of transit duration
being determined by the period $P$ via the Kepler law, one should use it to
reduce the Type II error. We note that Jeffrey’s prior can be unreasonable
even in unknown prior situations: if, in a given experiment, the posterior
volume is strongly varying across the parameter space, the Jeffrey’s prior is
very experiment specific. A prior that is smooth across the parameter space is
undoubtedly a better prior, even if we do not know what the specific form
should be. Nevertheless, when Jeffrey’s prior is reasonable, it simplifies the
analytic calculation of the p-value.

We show that both the p-value and the Bayes factor can be expressed as energy
versus entropy competition. We define energy as the logarithm of the
likelihood ratio and Bayesian entropy as the logarithm of the number of
posterior volumes that fit in the prior volume. The constant energy Bayes
factor is analogous to the thermodynamical free energy. Conversely, the
p-value in the asymptotic regime corresponds to the canonical partition
function, which is a Boltzmann factor weighted sum over the posterior states
with the test statistic above the observed one. In the low p-value regime,
only the states close to the observed test statistic contribute to the sum.
Therefore, the constant energy Bayes factor and the asymptotic p-value are
related. This also happens in physics, where the statistical and
thermodynamical definitions of the free energy coincide in the thermodynamical
limit. As an example, we show that the p-value of the standard $\chi^{2}$
distribution of $d$ degrees of freedom can be interpreted as energy versus
entropy competition, with the latter defined as the logarithm of the number of
states on the constant energy shell. Entropy grows as the log of the area of a
sphere in $d$ dimensions with a radius proportional to the square root of
energy.

The formalism developed here extends the Wilks’ theorem wilksTheorem  in
several ways. First, the connection to the Bayes factor allows us to define
the posterior volume beyond the Gaussian approximation inherent in the
asymptotic limit assumed for the Wilks’ theorem. Second, the Wilks’ theorem
assumes the parameter values are inside the boundaries. We show that a
generalization counting the states as a function of energy provides a proper
generalization that gives better results. Third, Wilks’ theorem does not
account for the multiplicity of the Look Elsewhere Effect and thus fails for
coordinate parameters, which do not change the expected energy. Our method
correctly handles these situations.

As an example application, we apply the formalism to the exoplanet transit
search in the stellar variability polluted data. We search for exoplanets by
scanning over the period, phase, and transit duration and we show that the
multiplicity from the Look Elsewhere effect is of order $10^{7}$. We find that
the Laplace approximation for uncertainty volume $V_{\rm post}$ is inaccurate,
while the numerical integration of the Bayesian evidence gives very accurate
results when compared to simulations. We emphasize the role of informative
priors such as planet transit time prior, which reduce the Type II error,
leading to a higher fraction of true planets discovered at a fixed Type I
error (p-value) threshold. The method enables fast evaluation of the false
positive rate for every exoplanet candidate without running expensive
simulations.

There are other practical implications that follow from our analysis. For
example, the multiplicity depends not only on the prior range but also on the
posterior error on the scanning parameters. If this error is small in one part
of the parameter space but large in others, then the likelihood ratio test
statistic leads to a large multiplicity that will increase the p-value for all
events. We show that analytic predictions reproduce the distribution of false
positives as a function of period and transit duration. An informed choice of
the prior guided by what we know about the problem and what our goals are may
change this balance and reduce the multiplicity penalty, thus reducing the
Type II error: Bayes factor can be a better test statistic for Type II error
than the likelihood ratio. One is of course not allowed to pick and choose the
prior a posteriori: we must choose it prior to the data analysis.

In many situations, it is possible to analytically obtain the false positive
rate as a function of Bayes factor test statistic from Equation (23), which
gives a p-value estimate that is more reliable for hypothesis testing than the
corresponding Bayes factor in situations of a new discovery where the prior is
not yet known. As a general recommendation, we thus advocate that Bayesian
analyses report the frequentist p-value using the Bayes factor test statistic
against the null hypothesis as a way to quantify the significance of a new
discovery, and that frequentist analyses use the Bayes factor as the optimal
test statistic for hypothesis testing while using frequentist methods to
quantify its significance.

###### Acknowledgements.

We thank Adrian Bayer and Natan Dominko-Kobilica for many valuable discussions
and comments. This material is based upon work supported in part by the
Heising-Simons Foundation grant 2021-3282 and by the U.S. Department of
Energy, Office of Science, Office of Advanced Scientific Computing Research
under Contract No. DE-AC02-05CH11231 at Lawrence Berkeley National Laboratory
to enable research for Data-intensive Machine Learning and Analysis.

##  Appendix A Jeffrey’s prior in the exoplanet example

Here, we derive the approximate Jeffrey’s prior in the exoplanet example by
analyzing the Fisher information matrix (10):

$I_{ij}({\bf z})=Re\int\frac{\partial\mathcal{F}\\{m\\}}{\partial
z_{i}}\frac{\partial\mathcal{F}\\{m\\}^{*}}{\partial
z_{j}}\frac{d\omega}{P(\omega)},$ (39)

Let us first determine the scaling of the $I_{\phi\phi}$. Fourier transforming
the signal from Eq. 34, differentiating it with respect to $\phi$, inserting
into equation 39 and expressing the amplitude $A=z_{1}$ with $E$ we obtain

$I_{\phi\phi}=2E\frac{\int\frac{\omega^{2}P^{2}|\mathcal{F}\\{m\\}|^{2}d\omega}{P(\omega)}}{\int\frac{|\mathcal{F}\\{m\\}|^{2}d\omega}{P(\omega)}}\propto\frac{P^{2}}{\tau^{2}}\,,$
(40)

where in the second step we change the variable of integration to the
dimensionless quantity $\omega\tau$ and neglect a residual mild dependence of
the dimensionless integral on $P$ and $\tau$.

In a scenario where $\tau$ is an independent quantity we similarly get
$I_{PP}\propto P^{-2}\tau^{-1}$ and $I_{\tau\tau}\propto\tau^{-1}$, therefore
the Jeffrey’s prior is:

$p({\bf z})\propto\det I^{1/2}\propto\tau^{-2}\,.$ (41)

In the other scenario, where $\tau$ is fixed by the period we have
$I_{PP}\propto P^{-2}$ and Jeffrey’s prior is:

$p({\bf z})\propto\tau_{K}^{-1}\propto P^{-1/3}\,.$ (42)

An estimate of the Jeffrey’s prior is obtained by normalizing in the given
range of the parameters.

## References

  * [1] Harold Jeffreys.  The theory of probability.  OUP Oxford, 1998.
  * [2] David MacKay.  Information theory, inference and learning algorithms.  Cambridge university press, 2003.
  * [3] Adrian E. Bayer and Uroš Seljak.  The look-elsewhere effect from a unified bayesian and frequentist perspective.  Journal of Cosmology and Astroparticle Physics, 2020(10):009–009, oct 2020.
  * [4] Rupert G. Miller.  Simultaneous Statistical Inference.  Springer New York, 1981.
  * [5] J P Shaffer.  Multiple hypothesis testing.  Annual Review of Psychology, 46(1):561–584, 1995.
  * [6] Vijay Balasubramanian.  Statistical inference, occam’s razor, and statistical mechanics on the space of probability distributions.  Neural computation, 9(2):349–368, 1997.
  * [7] Jakob Robnik and Uroš Seljak.  Matched filtering with non-Gaussian noise for planet transit detections.  Monthly Notices of the RAS, 504(4):5829–5839, July 2021.
  * [8] N. R. Lomb.  Least-Squares Frequency Analysis of Unequally Spaced Data.  Astrophysics and Space Science, 39(2):447–462, February 1976.
  * [9] J. D. Scargle.  Studies in astronomical time series analysis. II. Statistical aspects of spectral analysis of unevenly spaced data.  The Astrophysical Journal, 263:835–853, December 1982.
  * [10] Jacob T. VanderPlas.  Understanding the lomb–scargle periodogram.  The Astrophysical Journal Supplement Series, 236(1):16, May 2018\.
  * [11] Yee Whye Teh and Michael I Jordan.  Hierarchical bayesian nonparametric models with applications.  Bayesian nonparametrics, 1:158–207, 2010.
  * [12] George Casella.  An introduction to empirical bayes data analysis.  The American Statistician, 39(2):83–87, 1985.
  * [13] Robert E Kass and Adrian E Raftery.  Bayes factors.  Journal of the american statistical association, 90(430):773–795, 1995.
  * [14] Samuel S Wilks.  The large-sample distribution of the likelihood ratio for testing composite hypotheses.  The annals of mathematical statistics, 9(1):60–62, 1938.
  * [15] R. V. Baluev.  Assessing the statistical significance of periodogram peaks.  Monthly Notices of the Royal Astronomical Society, 385(3):1279–1285, Apr 2008.
  * [16] R. B. Davies.  Hypothesis testing when a nuisance parameter is present only under the alternative.  Biometrika, 64(2):247–254, 1977.
  * [17] Robert B. Davies.  Hypothesis testing when a nuisance parameter is present only under the alternative.  Biometrika, 74:33–43, 1987.
  * [18] Jon M. Jenkins, Peter Tenenbaum, Shawn Seader, Christopher J. Burke, Sean D. McCauliff, Jeffrey C. Smith, Joseph D. Twicken, and Hema Chandrasekaran.  Kepler Data Processing Handbook: Transiting Planet Search.  Kepler Science Document, January 2017.
  * [19] David G Koch, William J Borucki, Gibor Basri, Natalie M Batalha, Timothy M Brown, Douglas Caldwell, Jørgen Christensen-Dalsgaard, William D Cochran, Edna DeVore, Edward W Dunham, et al.  Kepler mission design, realized photometric performance, and early science.  The Astrophysical Journal Letters, 713(2):L79, 2010.
  * [20] George R. et al Ricker.  Transiting exoplanet survey satellite (tess).  Journal of Astronomical Telescopes, Instruments, and Systems, 1:014003, jan 2015.
  * [21] David M Kipping.  Efficient, uninformative sampling of limb darkening coefficients for two-parameter laws.  Monthly Notices of the Royal Astronomical Society, 435(3):2152–2160, 2013.
  * [22] Jakob Robnik and Uroš Seljak.  Kepler data analysis: Non-gaussian noise and fourier gaussian process analysis of stellar variability.  The Astronomical Journal, 159(5):224, 2020.
  * [23] AH Stroud and Don Secrest.  Approximate integration formulas for certain spherically symmetric regions.  Mathematics of Computation, 17(82):105–135, 1963.
  * [24] Nico Schlomer.  quadpy.
  * [25] David M. Kipping.  Bayesian priors for the eccentricity of transiting planets.  Monthly Notices of the Royal Astronomical Society, 444(3):2263–2269, 09 2014.
  * [26] D. V. Lindley.  A statistical paradox.  Biometrika, 44:187–192, 1957.