# A UNIVERSAL ROBUST BOUND FOR THE INTRINSIC BAYES FACTOR.

Richard Clare

Abstract

In this work, we undertake a comprehensive reformulation, modification, and
extension of Smith & Spiegelhalter’s [14] and [15] Bayes Factor work within
the evolving subject of Objective Bayes Factors. Our primary focus centers on
defining and computing empirical and theoretical bounds for the Intrinsic
Bayes Factor (IBF) across various models, including normal, exponential,
Poisson, geometric, linear, and ANOVA. We show that our new bounds are useful,
feasible, and change with the amount of information. We also propose a
methodology to construct the least favorable (for the null model) intrinsic
priors that result in the lower and upper bounds of the Intrinsic Bayes
Factors under certain conditions. Notably, our lower bounds exhibit superior
performance compared to the well-known $-ep\log(p)$ bound proposed by Sellke
et al. (2001) [13] based on p-values.

###### Contents

  1. 1 Introduction
  2. 2 LITERATURE REVIEW
    1. 2.1 Bayes factors
    2. 2.2 Prior and Posterior Model Probabilities
    3. 2.3 Jeffreys’s objective priors
    4. 2.4 Objective Bayes Factors
    5. 2.5 Minimal training samples
      1. 2.5.1 A formal definition for the set of all minimal training samples
      2. 2.5.2 An algorithm to generate the set of all minimal training samples
    6. 2.6 The Intrinsic Bayes Factor
    7. 2.7 The Expected Posterior Bayes Factors
  3. 3 A UNIVERSAL ROUBUST BOUND
    1. 3.1 The Smith and Spiegelhalter’s method
    2. 3.2 The universal robust bound for the Intrinsic Bayes Factors
    3. 3.3 A bridge between the universal bound and the IBF
      1. 3.3.1 A comparison with the p-value based lower bound
  4. 4 The SP Bayes Factors
    1. 4.1 The SP Prior
    2. 4.2 The SP Bayes Factor
    3. 4.3 Results
  5. 5 Hypothesis Testing with the Normal Distribution
    1. 5.1 The Normal Scale Hypothesis Test
      1. 5.1.1 The universal lower bound
      2. 5.1.2 Remarks
      3. 5.1.3 The SP and Intrinsic Priors
    2. 5.2 The Normal Mean Hypothesis Test
      1. 5.2.1 The Expected Intrinsic Bayes Factor for the Normal Mean Hypothesis Test
      2. 5.2.2 The Intrinsic Bayes Factor Bounds for the Normal Mean Hypothesis Test
      3. 5.2.3 Normal mean hypothesis test with unknown variance.
  6. 6 The Exponential Model Hypothesis Test
    1. 6.1 The Intrinsic Bayes Factor Dynamic and Universal Lower Bounds
    2. 6.2 The SP Bayes Factors
    3. 6.3 The EP Bayes Factors
    4. 6.4 The -eplogp lower bound
    5. 6.5 Conclusions
  7. 7 Linear regression and ANOVA Hypothesis Testing with Bayes Factors
    1. 7.1 Linear regression
    2. 7.2 Two nested normal-linear models
    3. 7.3 The One-way Layout (ANOVA)
      1. 7.3.1 The Full Jeffrey’s (no independence) $q_{k}=p_{k},$, for $k=i,j$ for ANOVA
      2. 7.3.2 The Reference prior ($q_{k}=0,k=i,j$) for ANOVA
      3. 7.3.3 The Modified Jeffrey’s Prior ($q_{j}=p_{j}-p_{i},q_{i}=0$) for ANOVA
    4. 7.4 Conclusions
  8. 8 Separated Hypothesis Testing
    1. 8.1 The Poisson vs. Geometric Separated Hypothesis Test
    2. 8.2 The Negative Binomial vs. Poisson Separated Hypothesis Test
    3. 8.3 Conclusions
  9. 9 Ghosh and Samanta Alternative Approach
    1. 9.1 The Normal Mean Hypothesis Test Example
    2. 9.2 The Exponential Hypothesis Test Example
      1. 9.2.1 Conclusions
  10. 10 Research Contributions

##  Chapter 1 Introduction

The misuse of classical hypothesis testing methods, specifically p-values, has
been the subject of substantial critique in the statistical literature
(Wasserstein and Lazar, 2016 [17]; Wasserstein et al., 2019 [16]). Despite
these warnings, the practice of p-value null hypothesis significance testing
(NHST) continues to be a predominant technique for model selection. However,
alternatives do exist. Among them, the Bayesian testing methodology developed
by Sir Harold Jeffreys (1961) [7] stands as one of the most widely accepted
and used. Jeffreys’s methodology, premised on the concept that the null
hypothesis carries a positive point mass, quantifies the support for or
against a model based on the data. Additionally, it offers practical
applications through derived approximations for complex or analytically
infeasible solutions. A significant contribution of Jeffreys’s framework is
the Bayes factor, a measure of support for one model over another,
irrespective of the absolute correctness of these models. Its effectiveness
spans various inference desiderata, including interpretability, adherence to
Occam’s razor, and consistency in large samples. Kass and Raftery (1995) [8]
offer an extensive introduction to Bayes factors. Our research primarily
concentrates on the pivotal role of training samples within statistical
methodologies. These methodologies include classification and discrimination,
robustness, model selection, and cross-validation. Recent advancements in
Bayesian model selection, namely the intrinsic Bayes factor by Berger and
Pericchi (1996a) [3] and the expected posterior prior by Perez and Berger
(2002) [10], have employed training samples to transition improper objective
priors into more suitable distributions for model selection.

A challenge often faced is the indeterminacy of classical Bayes Factors, due
to the dependence on an undefined ratio of constants derived from the improper
priors. To circumnavigate this, Jeffrey suggested conventional proper priors
for the extra-parameters under the larger hypothesis, albeit improper for
common parameters. Various techniques have been proposed to address this
issue, including those by Smith and Spiegelhalter (1980) [14] and
Spiegelhalter and Smith (1982) [15]. Though approximate, these techniques
facilitate sensible scaling for Bayes Factors and bear a direct relationship
with more recent approaches like the Intrinsic Bayes Factor, the Intrinsic
Priors, and EP-priors. To elucidate these connections, we will re-
conceptualize Smith and Spiegelhalter’s approach and the associated Bayes
Factor, hereby referred to as SSBF’s ($B^{SS}$), drawing parallels to other
methodologies.

Chapter three specifically focuses on the quasi-Bayes factor, which emerges as
a universal lower bound for the Intrinsic Bayes factor. In Chapter four we
introduce the concept of SP Priors (Superior Posterior priors) that either
generate the Intrinsic Bayes Factors bound under specific assumptions or serve
as an asymptotic approximation of the bounds. These priors are the least
favorable priors for the null hypothesis and through chapter three to chapter
eight we performed numerical simulations and verified the consistency of the
SP, the empirical and the theoretical Intrinsic Bayes Factor bounds across
various essential distributions such as the normal, exponential, Poisson,
geometric, and negative binomial. Furthermore, we calculated the bounds for
two nested normal-linear models and ANOVA. We discussed model selection
problems across different distributions and hypothesis tests such as nested
and separated hypotheses, and for those cases we provided closed forms for the
Intrinsic Bayes factors Bounds, the SP Priors, and the SP Bayes Factors.

##  Chapter 2 LITERATURE REVIEW

###  2.1 Bayes factors

Traditional point estimation methods, such as the maximum likelihood estimator
(MLE) or the method of moments (MoM), typically involve the utilization of
specific estimates for each model parameter based on the available data. In
contrast, Bayesian inference takes a different approach by assigning a
probability distribution to each model parameter.

Let $\textbf{y}=(y_{1},...,y_{n})$ represent a set of $n$ random samples
obtained from a random variable $Y$ with a density $f_{i}(y|\theta_{i})$,
where $\theta_{i}$ denotes a vector of parameters associated with model
$M_{i}$. Employing the Bayes’ Theorem, the posterior distribution for model
$M_{i}$ can be expressed as:

$\pi(\boldsymbol{\theta}_{i}|\textbf{y})=\displaystyle\frac{f_{i}(y|\theta))\pi_{i}(\boldsymbol{\theta}_{i})}{m_{i}(\textbf{y})}=\displaystyle\frac{f_{i}(\textbf{y}|\theta_{i})\pi(\boldsymbol{\theta}_{i})}{\displaystyle\int
f_{i}(\textbf{y}|\theta_{i})\pi(\boldsymbol{\theta}_{i})d\boldsymbol{\theta}_{i}}$

where $m_{i}(\textbf{y})$ is called the marginal, evidence, or global
likelihood, and $\pi_{i}(\boldsymbol{\theta})$ is the prior distribution of
the parameters for $M_{i}$. Suppose that we are comparing two models for the
data y

$M_{k}:\textbf{Y}\mbox{ has a density
}f_{k}(\textbf{y}|\boldsymbol{\theta}_{k})\hskip 5.0ptk=i,j$

where $\boldsymbol{\theta}_{k}$ are unknown model parameters associated with
model $M_{k}$. Suppose we have available prior distributions
$\pi_{k}(\boldsymbol{\theta}_{k})$ for the unknown parameters. Define the
marginal or predictive densities of y as

$m_{i}(\textbf{y})=\int
f_{i}(\textbf{y}|\boldsymbol{\theta}_{i})\pi_{i}(\boldsymbol{\theta}_{i})d\boldsymbol{\theta}_{i}.$

The _Bayes factor_ of $M_{i}$ to $M_{j}$ is given by

$B_{ij}=\displaystyle\frac{m_{i}(\textbf{y})}{m_{j}(\textbf{y})}=\displaystyle\frac{\displaystyle\int
f_{i}(\textbf{y}|\boldsymbol{\theta}_{i})\pi_{i}(\boldsymbol{\theta}_{i})d\boldsymbol{\theta}_{i}}{\displaystyle\int
f_{j}(\textbf{y}|\boldsymbol{\theta}_{j})\pi_{j}(\boldsymbol{\theta}_{j})d\boldsymbol{\theta}_{j}}.$

The _Bayes factor_ is often interpreted as the evidence provided by the data
for the model $M_{i}$ vs. the alternative model $M_{j}$, but the _Bayes
factor_ depends also on the _priors_. Alternatively, $B_{ij}$ can be
interpreted as the weighted likelihood ratio of $M_{i}$ to $M_{j}$. with the
priors being the "weighting functions."

###  2.2 Prior and Posterior Model Probabilities

Suppose that prior probabilities $\pi(M_{i}),i=1,...,n$, of the $n$ models,
are available, then we can obtain for each model $M_{i}$ the posterior
probability from the Bayes factor given the data y by calculating

$P(M_{i}|\textbf{y})=\frac{\pi(M_{i})m_{i}(\textbf{y})}{\displaystyle\sum_{k=1}^{n}\pi(M_{k})m_{k}(\textbf{y})}$

This equation can also be written as

$P(M_{i}|\textbf{y})=\Big{[}\sum_{k=1}^{n}\frac{\pi(M_{k})}{\pi(M_{i})}B_{ki}(\textbf{y})\Big{]}^{-1}$

where $O_{ki}=\pi(M_{k})/\pi(M_{i})$ is known as the prior odds of the model
$M_{k}$ over the model $M_{i}$, and the posterior probabilities odds are given
by

$\frac{P(M_{i}|\textbf{y})}{P(M_{j}|\textbf{y})}=\frac{\pi(M_{i})m_{i}(\textbf{y})}{\pi(M_{j})m_{j}(\textbf{y})}=O_{ij}B_{ij}$

Hence, the posterior odds are given the product of the prior odds and the
Bayes factor. A common unbiased choice of the prior model probabilities is
given by $\pi(M_{i})=1/n$ where each model has the same prior probability.
Note that using those uniform priors we can obtain the re-normalized marginal
distributions from the posterior model probabilities as

$\overline{m}_{i}(\textbf{y})=\frac{m_{i}(\textbf{y})}{\sum_{k=1}^{n}m_{k}(\textbf{y})}$

and also the Bayes Factor can be in terms of the re-normalized marginals as

$B_{ij}(\textbf{y})=\frac{\overline{m}_{i}(\textbf{y})}{\overline{m}_{j}(\textbf{y})}$

###  2.3 Jeffreys’s objective priors

Consider the Fisher information matrix when there are n parameters, so that
$\boldsymbol{\theta}$ is a $n\times 1$ vector
$\boldsymbol{\theta}=\begin{bmatrix}\theta_{1},\theta_{2},\dots,\theta_{n}\end{bmatrix}^{\mathrm{T}}$,
then the Fisher information takes the form of an $n\times n$ matrix with
typical element

${I_{ij}\left(\theta\right)=\operatorname{E}\left[\left.\left(\frac{\partial}{\partial\theta_{i}}\log
f(\textbf{y}|\boldsymbol{\theta})\right)\left(\frac{\partial}{\partial\theta_{j}}\log
f(\textbf{y}|\boldsymbol{\theta})\right)\right|\boldsymbol{\theta}\right]}.$

Jeffreys prior is a non-informative (objective) prior distribution for a
parameter space; it is proportional to the square root of the determinant of
the Fisher information matrix

$\pi^{N}\left(\boldsymbol{\theta}\right)\propto\sqrt{\det
I\left(\boldsymbol{\theta}\right)}\,$

and this prior is invariant under reparameterization of the parameter vector
$\boldsymbol{\theta}$ and this property makes it of special interest for use
with location and scale parameters.

###  2.4 Objective Bayes Factors

The unscaled Bayes factors are obtained using the non-informative objective
priors $\pi_{i}^{N}(\theta_{i})$ and $\pi_{j}^{N}(\theta_{j})$ for the models
$M_{i}$ and $M_{j}$ respectively and can be calculated as

$B^{N}_{ij}(\textbf{y})=\frac{\displaystyle\int
f_{i}(\textbf{y}|\theta_{i})\pi_{i}^{N}(\theta_{i})d\theta_{i}}{\displaystyle\int
f_{j}(\textbf{y}|\theta_{j})\pi_{j}^{N}(\theta_{j})d\theta_{j}}$

Non-informative priors can be used for one-sided hypothesis testing but
unfortunately, when we are using improper priors for point null hypotheses the
Bayes Factor depends on an undetermined ratio of constants $C_{i}/C_{j}$,
which comes from the improper priors. Some techniques are widely used to
overcome this issue which can approximate the desired ratio of constants and
those will be studied and extended in the next sections.

###  2.5 Minimal training samples

There are different approaches to calculate Bayes Factors based on the concept
of a proper ’minimal training sample’ which is a subset of the entire data y.
Definition: A training sample $\textbf{y}(\ell)\subset\textbf{y}$, is called
proper if $0<m(\textbf{y}(\ell))<\infty$, and minimal if it is proper and no
subset is proper.

####  2.5.1 A formal definition for the set of all minimal training samples

The formal definition and formula for the set of all possible minimal training
samples can be obtained using the definition above and extending the concept.
Let
$y(\ell_{i})=y(i_{1},i_{2},...,i_{k})=\\{y_{i_{1}},y_{i_{2}},...,y_{i_{k}}\\}$
be any minimal training sample with $k\leq n$. The set
$\\{y(\ell_{1}),...,y(\ell_{L})\\}$ of all possible minimal training samples
(non-repeated subsets of size $k$) can be defined as the set

$D_{n,k}=\Bigg{\\{}\Big{\\{}\big{\\{}\cdot\cdot\cdot\\{y(i_{1},i_{2},...,i_{k})\\}_{i_{k}=i_{(k-1)}}^{i_{n}}\cdot\cdot\cdot\big{\\}}^{i_{n}}_{i_{2}=i_{1}}\Big{\\}}_{i_{1}=1}^{i_{n}}:1\leq
i_{1}<i_{2}<i_{3}<\dots<i_{k}\leq n\Bigg{\\}}$

Hence,

$|D_{n,k}|={n\choose k}=L\mbox{ \hskip 10.0pt and \hskip 10.0pt
}D_{\infty,k}=\lim_{n\to\infty}D_{n,k}$

Therefore,

$|D_{\infty,k}|=|\lim_{n\to\infty}D_{n,k}|=\lim_{n\to\infty}|D_{n,k}|=\lim_{n\to\infty}\frac{(n)(n-1)(n-(k+1))}{k!}=\infty$

For ease of notation, we will denote the set $D_{\infty,k}$ as $D$ throughout
this thesis.

The determination of the set $D_{n,k}$, encompassing all minimal training
samples of size $k$ from a sample of size $n$, will be facilitated by an
original, exhaustive, albeit sub-optimal algorithm. This algorithm yields the
computation of Bayes factors involving these terms?such as the Arithmetic
Intrinsic Bayes Factor, Median Intrinsic Bayes Factor, EP-Bayes factors, and
SS-Bayes factors?which will be expounded upon in subsequent chapters.

Similarly, considering $L=$ ${n}\choose{k}$,, we can derive the following
computationally intuitive formula applicable to the Arithmetic Intrinsic Bayes
Factor (AIBF)

$\boxed{\frac{1}{L}\displaystyle\sum_{i_{1}=1}^{i_{n}}\sum_{i_{2}=i_{1}}^{i_{n}}\cdot\cdot\cdot\sum_{i_{k}=i_{(k-1)}}^{i_{n}}\mathds{1}(1\leq
i_{1}<i_{2}<\dots<i_{k}\leq n)B_{01}(y(i_{1},...,i_{k}))=\frac{1}{L}\sum_{d\in
D_{n,k}}B_{01}(d).}$

In the next section, we will present the general definition algorithm that can
receive a vector of observations, and will generate the set of all possible
minimal training samples.

####  2.5.2 An algorithm to generate the set of all minimal training samples

We present here an original algorithm devised for the purpose of generating
the set of all minimal training samples. It is important to note that the
algorithm proposed in this study may not be suitable for processing very large
sample sizes due to its computational complexity, which amounts to
approximately $O(n^{k})$ operations. For future investigations, it is
recommended to explore more efficient algorithms for generating the matrix M
based on the algorithm discussed herein.

Alternatively, a generalized version of this algorithm can be constructed
using a recurrence relation instead of nested for loops. Such an approach
would enable the development of a versatile algorithm applicable to any given
values of $n$ and $k$. In cases where the sample sizes are substantial, it is
advisable to employ randomization techniques to generate a representative
subset of minimal training samples, which can yield satisfactory
approximations of the Intrinsic Bayes Factors (IBFs).

Algorithm 1 An algorithm to produce the set of all possible minimal training
samples

$n,k\in\mathbb{N},n\geq k$

A matrix with k columns and ${n}\choose{k}$ rows with the indexes of all of
the possible minimal training samples of size k from a sample of size n.

function MTS($(n,k)$)

Initialize an empty matrix M

for $i_{1}\mbox{ in }1:n$ do

for $i_{2}\mbox{ in }i_{1}:n$ do

⋮

for $i_{k}\mbox{ in }i_{k-1}:n$ do

if $1\leq i_{1}<i_{2}<...<i_{k}\leq n$ then

Add the row with values $(i_{1},...,i_{k})$ to the bottom of the matrix M.

end if

end for

⋮

end for

end for

Return M

end function

###  2.6 The Intrinsic Bayes Factor

Suppose that non-informative (usually improper) priors
$\pi_{i}^{N}(\boldsymbol{\theta}_{i}),i=1,...,n$ are available for the models
$M_{1},...M_{n}$ respectively. The general recommendation is that we choose
those priors as "reference priors" (Berger and Bernardo 1992 [2]), but there
are many other choices and approximations with good results and advantages
over the reference priors in terms of computational and mathematical
complexity.
The _Intrinsic Bayes Factor_ (IBF) is calculated using a minimal training
sample $\textbf{y}(\ell)$ which converts the improper prior
$\pi_{i}^{N}(\boldsymbol{\theta}_{i})$ to proper posteriors
$\pi_{i}^{N}(\boldsymbol{\theta}_{i}|\textbf{y}(\ell))$. The Bayes factor is
calculated using this proper posterior as the prior distribution of the
remaining of the data $\textbf{y}(-\ell)=\textbf{y}\setminus\textbf{y}(\ell)$
as follows

$\boxed{B_{ij}^{I}=B_{ij}(\ell)=\frac{\displaystyle\int
f_{i}(\textbf{y}(-\ell)|\boldsymbol{\theta}_{i},\textbf{y}(\ell))\pi_{i}(\boldsymbol{\theta}_{i}|\textbf{y}(\ell))d\boldsymbol{\theta}_{i}}{\displaystyle\int
f_{j}(\textbf{y}(-\ell)|\boldsymbol{\theta}_{j},\textbf{y}(\ell))\pi_{j}(\boldsymbol{\theta}_{i}|\textbf{y}(\ell))d\boldsymbol{\theta}_{j}}}$

where $B_{ij}(\ell)=B_{ij}(\textbf{y}(-\ell)|\textbf{y}(\ell))$. There are
$L={n\choose k}$ different minimal training samples for $n$ i.i.d samples,
where $k$ is the number of parameters of the density under consideration.
Denote by $B^{N}_{ij}(\mathbf{y})$, the Bayes Factor with the whole sample
$\mathbf{y}$ and un-scaled objective (Non-Informative) priors
$\pi^{N}_{k}(\boldsymbol{\theta}_{k}),k=i,j$, and,
$B_{ij}(\mathbf{y(-\ell)}|\mathbf{y(\ell)})$, the "trained" Bayes Factor after
the training sample has been employed to get a proper scaling of the Bayes
Factor. A fundamental Lemma with Bayes Factors is the following (Berger and
Pericchi(1996) [3]

######  Lemma 2.6.1.

Let $\textbf{y}=\\{y_{1},...,y_{n}\\}$ be independent and identically
distributed random samples, and assume that $\textbf{y}(\ell)\in D$ with
$\textbf{y}(-\ell)\cup\textbf{y}(\ell)=\textbf{y}$, then

$\mbox{
}B_{ij}(\ell)=B^{N}_{ij}(\mathbf{y})B^{N}_{ji}(\mathbf{y(\ell)}),\forall\ell.$
(2.1)

###### Proof.

The Bayes Factor of the remaining samples given the training samples can be
written as

$B_{ij}(\mathbf{y(-\ell)}|\mathbf{y(\ell)})=\frac{m_{i}(\mathbf{y(-\ell)}|\mathbf{y(\ell)})}{m_{j}(\mathbf{y(-\ell)}|\mathbf{y(\ell)})}=\frac{\displaystyle\int
f_{i}(\mathbf{y(-\ell)}|\theta_{i})\pi(\theta_{i}|\mathbf{y(\ell)})d\theta_{i}}{\displaystyle\int
f_{i}(\mathbf{y(-\ell)}|\theta_{j})\pi(\theta_{j}|\mathbf{y(\ell)})d\theta_{j}}$

The proper posterior distribution of the parameters given the MTS can be
written as

$\pi(\theta_{k}|\mathbf{y(\ell)})=\frac{f_{k}(\theta_{k}|\mathbf{y(\ell)})\pi_{k}^{N}(\theta_{k})}{m_{k}^{N}(\mathbf{y(\ell)})},\mbox{
for }k=i,j$

Hence,

$\frac{\displaystyle\int
f_{i}(\mathbf{y(-\ell)}|\theta_{i})\pi(\theta_{i}|\mathbf{y(\ell)})d\theta_{i}}{\displaystyle\int
f_{j}(\mathbf{y(-\ell)}|\theta_{j})\pi(\theta_{j}|\mathbf{y(\ell)})d\theta_{j}}=\frac{\displaystyle\int
f_{i}(\mathbf{y(-\ell)}|\theta_{i})\frac{f_{i}(\theta_{i}|\mathbf{y(\ell)})\pi_{i}^{N}(\theta_{i})}{m_{i}^{N}(\mathbf{y(\ell)})}d\theta_{i}}{\displaystyle\int
f_{j}(\mathbf{y(-\ell)}|\theta_{j})\frac{f_{j}(\theta_{j}|\mathbf{y(\ell)})\pi_{j}^{N}(\theta_{j})}{m_{j}^{N}(\mathbf{y(\ell)})}d\theta_{j}}$
$=\frac{m_{j}^{N}(\mathbf{y(\ell)})}{m_{i}^{N}(\mathbf{y(\ell)})}\frac{\displaystyle\int
f_{i}(\mathbf{y(-\ell)}|\theta_{i})f_{i}(\theta_{i}|\mathbf{y(\ell)})\pi_{i}^{N}(\theta_{i})d\theta_{i}}{\displaystyle\int
f_{j}(\mathbf{y(-\ell)}|\theta_{j})f_{j}(\theta_{j}|\mathbf{y(\ell)})\pi_{j}^{N}(\theta_{j})d\theta_{j}}$
(2.2)

Since
$f_{k}(\mathbf{y}(\ell)|\boldsymbol{\theta}_{k})f_{k}(\mathbf{y}(-\ell)|\boldsymbol{\theta}_{k})=f_{k}(\mathbf{y}|\boldsymbol{\theta}_{k})\mbox{
for }k=i,j$, the ratio of the marginals is

$\frac{m_{j}^{N}(\mathbf{y(\ell)})}{m_{i}^{N}(\mathbf{y(\ell)})}=B_{ji}^{N}(\mathbf{y(\ell)})$

Hence, equation (2) becomes

$B_{ji}^{N}(\mathbf{y(\ell)})\frac{\displaystyle\int
f_{i}(\mathbf{y}|\theta_{i})\pi_{i}^{N}(\theta_{i})d\theta_{i}}{\displaystyle\int
f_{j}(\mathbf{y}|\theta_{j})\pi_{j}^{N}(\theta_{j})d\theta_{j}}=B_{ji}^{N}(\mathbf{y(\ell)})B_{ij}^{N}(\mathbf{y}).$

∎

By employing the methodology elucidated in the aforementioned lemmas, we have
successfully mitigated the reliance of $B_{ij}$ on the scales of the priors,
thereby resolving the issue of undetermined constants associated with the
Bayes factor. However, a new dependence arises, contingent upon the arbitrary
selection of the minimal training sample denoted as $\textbf{y}(\ell)$. To
address this dependency and enhance stability, we propose to perform an
"average" of $B_{ij}(l)$ over all feasible training samples
$\textbf{y}(\ell)$, where $\ell=1,...,L$.

###  2.7 The Expected Posterior Bayes Factors

Suppose that non-informative (typically improper) priors
$\pi_{i}^{N}(\boldsymbol{\theta}_{i})$, where $i=1,\ldots,n$, associated with
the models $M_{1},\ldots,M_{n}$, respectively are available. The concept of
Expected Posterior (EP) priors, introduced by Perez and Berger in 2002 [10],
allows for the calculation of proper posteriors
$\pi_{i}^{N}(\boldsymbol{\theta}_{i}|\textbf{y}(\ell))$ using a minimal
training sample $\textbf{y}(\ell)$. The EP prior for model $M_{i}$ is then
defined as:

$\boxed{\pi_{i}^{EP}(\theta_{i})=\frac{1}{L}\displaystyle\sum_{\ell=1}^{L}\pi_{i}^{N}(\boldsymbol{\theta}_{i}|\textbf{y}(\ell))}$

The _EP Bayes factor_ , which enables the comparison between models $M_{i}$
and $M_{j}$, is computed as follows:

$\boxed{B_{ij}^{EP}(\textbf{y})=\frac{\displaystyle\int
f_{i}(\textbf{y}|\boldsymbol{\theta}_{i})\pi_{i}^{EP}(\boldsymbol{\theta}_{i})d\boldsymbol{\theta}_{i}}{\displaystyle\int
f_{j}(\textbf{y}|\boldsymbol{\theta}_{j})\pi^{EP}_{j}(\boldsymbol{\theta}_{j})d\boldsymbol{\theta}_{j}}}$

In the above equations, $f_{i}(\textbf{y}|\boldsymbol{\theta}_{i})$ and
$f_{j}(\textbf{y}|\boldsymbol{\theta}_{j})$ represent the likelihood functions
associated with models $M_{i}$ and $M_{j}$, respectively. The EP prior,
derived from improper priors using the minimal training sample
$\textbf{y}(\ell)$, allows for the calculation of the EP Bayes factor,
facilitating model comparison and selection.

##  Chapter 3 A UNIVERSAL ROUBUST BOUND

###  3.1 The Smith and Spiegelhalter’s method

This chapter presents an innovative methodology for computing Bayes Factors,
building upon the seminal work of Spiegelhalter and Smith (1982) [15]. Their
approach, which encompassed examining an alternative prior specification,
could be interpreted as a genuine subjective Bayesian analysis based on
specific prior beliefs, or alternatively, as a formal analysis serving as a
theoretical tool for comparing a model $M_{0}$ with a localized subset of
models within another model $M_{1}$. Spiegelhalter and Smith employed the
concept of "imaginary training samples" to maximize support for the simpler
model $M_{0}$. While their paper primarily focused on assigning the constant
$c_{0}/c_{1}$ in equation (5), Lempers (1971) [9] and Atkinson (1978) [1]
emphasized that assigning vague prior information directly to $c_{0}$ and
$c_{1}$ is inherently arbitrary without incorporating a training sample to
determine appropriate variable scaling. Consequently, they proposed
partitioning a portion of the data for a pre-comparison inference phase,
enabling the formation of suitable priors for the parameters within each
model.

The method suggested by Spiegelhalter and Smith involved determining
$\textbf{y}_{0}$ as the argument that maximizes
${m_{0}^{N}(\textbf{y}(\ell))}/{m_{1}^{N}(\textbf{y}(\ell))}$ and subsequently
solving the equation

$cB_{01}^{N}(\textbf{y}_{0})=1.$

This approach yielded the Bayes factor expressed as

$\boxed{B^{SS}_{10}(\textbf{y})=cB^{N}_{10}(\textbf{y})=B^{N}_{01}(\textbf{y}_{0})B_{10}^{N}(\textbf{y}).}$

While Spiegelhalter and Smith’s technique focused exclusively on nested point
hypotheses, the upcoming sections will expand upon this method to encompass
more general cases and examples. Furthermore, a novel lower bound will be
introduced, leveraging the concept of "imaginary training samples" to maximize
support for $M_{0}$. However, there is a distinction between this proposed
method and the aforementioned approach by Smith and Spiegelhalter. Our method
does not concern itself with the constant $c$; instead, we aim to establish a
lower bound for various forms of Intrinsic Bayes Factors, including the
arithmetic, geometric, expected and harmonic. This bound encompasses both
theoretical and empirical versions, with the empirical version providing the
tightest bound. The subsequent section will outline the procedures for
obtaining both theoretical and empirical bounds.

###  3.2 The universal robust bound for the Intrinsic Bayes Factors

In this study, we propose a novel method that generates a lower bound for the
Bayes factor and corresponding priors, employing an approach that is closely
related to the concept of Intrinsic Bayes Factors and the technique employed
by Spiegelhalter and Smith. Let $\textbf{y}(\ell)$ denote a minimal training
sample, and $\pi(\boldsymbol{\theta_{i}}|\textbf{y}(\ell))$ represent the
posterior distribution of parameter $\boldsymbol{\theta}_{i}$ given the
minimal training sample for two models, $M_{i}$ and $M_{j}$. Our aim is to
compare models $M_{i}$ and $M_{j}$. To achieve this comparison, we define the
set $D$ as the aggregate of all possible minimal training samples with size
$k$, where $k$ is determined as the maximum of
$\\{\dim(\boldsymbol{\theta_{i}}),\dim(\boldsymbol{\theta_{j}})\\}$.

######  Definition 1.

The Intrinsic Bayes Factor Lower Bound ($\underline{B}^{I}_{01}$), is the
theoretical lower bound of all possible Intrinsic Bayes Factors obtained by
computing the infimum over the set of all possible imaginary training samples
D.

######  Definition 2.

The Empirical Intrinsic Bayes Factor Lower Bound
($\underline{B}^{I^{*}}_{01}$) is the empirical lower bound obtained by
computing the minimum over the set of all possible empirical training samples
$D_{n,k}$.

###  3.3 A bridge between the universal bound and the IBF

In this section, we will explore various approaches to construct lower bounds
for the Intrinsic Bayes factor, ultimately contributing to the field of model
comparison and selection. There is a bridge between IBF and the Least
Favorable Intrinsic Bayes Factor approach that we want to explore by
considering the AIBF, defined as,

$\boxed{B^{AI}_{10}(\textbf{y})=B^{N}_{10}(\mathbf{y})\times\frac{1}{L}\sum_{\ell=1}^{L}B^{N}_{01}(\mathbf{y(\ell)})}$
(3.1)

where $L={n\choose m}$, where m is the minimal training sample size. The AIBF
is asymmetric the more complex model has to be placed above to guarantee the
convergence.

$\boxed{\overline{B}^{I}_{10}:=B^{N}_{10}\sup_{\textbf{y}(\ell)\in
D}B_{01}(\textbf{y}(\ell))\geq B^{N}_{10}\frac{1}{L}\sum
B_{01}(\textbf{y}(\ell))=B^{A}_{10}.}$ (3.2)

By forming the reciprocals, we have the important consequence for a bound on
the probability of a null hypothesis:

$\boxed{\underline{B}^{I}_{01}:=B^{N}_{01}\inf_{\ell}B_{10}(\textbf{y}(\ell))=B^{N}_{01}\frac{1}{\sup_{\textbf{y}(\ell)\in
D}B_{01}(\textbf{y}(\ell))}\leq B^{N}_{01}\frac{1}{\frac{1}{L}\sum
B_{01}(\textbf{y}(\ell))}=B^{AI}_{01}.}$ (3.3)

where $\textbf{y}(\ell)$ is a theoretical training sample from the set of all
possible training samples. By construction, all Bayes Factors obey:
$B_{01}=1/B_{10}$ and $\overline{B}^{I}_{10}\geq B^{AI}_{10}$, since the
supremum is necessarily bigger or equal than the arithmetic average. On the
other hand, for generalizing the bound we may consider another bound based on
empirical training samples,

$B^{AI}_{01}\geq B^{N}_{01}\inf_{\textbf{y}(\ell)\in
D_{n,k}}B_{10}(\textbf{y}(\ell))=\underline{B}^{I^{*}}_{01},$ (3.4)

where $\textbf{y}(\ell)$ are empirical real training samples, and
$\underline{B}_{01}^{I^{*}}$ is formed over all the real training samples, of
a minimal size $m$. Clearly, $\overline{B}^{I*}_{10}\leq\overline{B}^{I}_{10}$
and $B^{AIBF}_{01}\geq\underline{B}^{I*}_{01}\geq\underline{B}^{I}_{01}.$
Thus, $\underline{B}^{I^{*}}_{01}$ and $\underline{B}^{I}_{01}$ are the least
favorable lower bounds on the Bayes Factor. In addition, the results obtained
with the AIBF are valid for other related types of Intrinsic Bayes Factors.

####  3.3.1 A comparison with the p-value based lower bound

In this section, we aim to illustrate the superiority of our bound in
comparison to the well-known p-value-based robust lower bound introduced by
Sellke et al. (2001) [13]. To commence our discussion, we will first establish
the definition of a p-value.

######  Definition 3.

A p-value $p(\textbf{Y})$ is a statistic satisfying $0\leq p(\textbf{y})\leq
1$ for every sample point y. Small values of p(Y) give evidence that
$H_{1}:\theta\in\Theta_{0}^{c}$ is true, where $\Theta_{0}$ is some subset of
the parameter space and $\Theta_{0}^{c}$ is its complement. A p-value is valid
if, for every $\theta\in\Theta_{0}$ and every $0\leq\alpha\leq 1$,

$P_{\theta}(p(\textbf{Y})\leq\alpha)\leq\alpha.$

The asymptotic behavior of the p-value as the sample size tends toward
infinity is contingent upon both the null hypothesis and the fundamental
distribution of the test statistic. Generally speaking, in the context of a
two-sided hypothesis test, if the null hypothesis holds true, the p-value
follows a uniform distribution within the range [0,1]. This implies that,
irrespective of the sample size, the p-value will be uniformly distributed
provided that the null hypothesis is valid.

The _Robust Lower Bound_ offers a framework for calibrating p-values and is
defined as follows:

$B_{01}^{L}(p)=\begin{cases}-ep\log(p)&\text{if }p<\frac{1}{e}\\\ 1&\text{if
}p\geq\frac{1}{e}\end{cases}$ (3.5)

We postulate that given any sample size $n$ and for a fixed $p$, our IBF Lower
Bound as a function of the p-value, $\underline{B}^{I}_{01}(p)$, is superior
to the robust lower bound $B_{01}^{L}(p)$. It can be mathematically expressed
as

$\underline{B}^{I}_{01}(p)\geq\underline{B}^{I}_{01}(p)\geq B_{01}^{L}(p).$

The last part of this inequality was proven by Sellke et. al (2001) [13] and
it should be noted that $\underline{B}^{I}_{01}(p)$ is dynamic, changing with
the sample size, akin to a real Bayes Factor, whereas $B_{01}^{L}(p)$ remains
static with $n$. In the next theorem, we use the Laplace approximation to show
that the robust lower bound stays constant under $H_{0}$ while our lower bound
converges to the correct decision as the sample size increases. The Laplace
approximation is a method used to approximate a probability distribution with
a normal distribution and it’s valid under the following regularity
conditions:

  1. 1.

Smoothness of the Logarithm of the Density: For a density function $f(y)$, if
$\log(f(y))$ is twice continuously differentiable in a neighborhood around the
mode of the distribution, then the Laplace approximation tends to work well.

  2. 2.

Existence of a Unique Mode: The distribution should have a single, well-
defined mode within the region of interest.

  3. 3.

Curvature of the Log Density: The curvature (second derivative) of
$\log(f(y))$ at the mode should not be zero. A non-zero curvature ensures that
the distribution can be well-approximated by a quadratic form, which is the
characteristic shape of a normal distribution.

  4. 4.

Convergence to Normality: As the sample size tends to infinity, the
distribution should converge to a normal distribution in a suitable sense
(like the Central Limit Theorem).

These conditions, particularly the smoothness and curvature of the log
density, are crucial for the Laplace approximation to accurately represent the
shape of the distribution around its mode.

However, it’s important to note that while the Laplace approximation is often
effective for many distributions that meet these conditions, it might not work
well for highly skewed or heavy-tailed distributions where the normality
assumption might not hold. In such cases, other approximation methods or
numerical techniques might be more appropriate.

######  Theorem 3.3.1.

Let $f_{0}(\textbf{y}|\boldsymbol{\theta_{0}})$ and
$f_{1}(\textbf{y}|\boldsymbol{\theta_{1}})$ be the densities that satisfies
regularity conditions (1)-(4). Let $H_{0}\subset H_{1}$,
$k_{0}=dim(\mathbf{\theta}_{0})$ and $k_{1}=dim(\mathbf{\theta}_{1})$ where
$k_{0}<k_{1}$, then under $H_{0}$

$\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}\to 0\mbox{ as }n\to\infty$

###### Proof.

For a proper minimal training sample $y^{*}(\ell),\exists b,c$ such that
$0<b\leq c<\infty$ and

$0<b\leq B_{10}^{N}(\textbf{y}^{*}(\ell))\leq c$

The Laplace approximation for $B_{01}^{N}(\textbf{y})$ (Berger & Pericchi
(2001) [4]) yields

$\Big{(}\frac{f_{0}(\textbf{y}|\hat{\theta}_{0})|\hat{I}_{0}|^{-1/2}}{f_{1}(\textbf{y}|\hat{\theta}_{1})|\hat{I}_{1}|^{-1/2}}\cdot\frac{\pi^{N}_{0}(\hat{\theta}_{0})}{\pi^{N}_{1}(\hat{\theta}_{1})}\Big{)}\cdot
b\leq
B_{01}^{N}(\textbf{y})B_{10}^{N}(\textbf{y}^{*}(\ell))\leq\Big{(}\frac{f_{0}(\textbf{y}|\hat{\theta}_{0})|\hat{I}_{0}|^{-1/2}}{f_{1}(\textbf{y}|\hat{\theta}_{1})|\hat{I}_{1}|^{-1/2}}\cdot\frac{\pi^{N}_{0}(\hat{\theta}_{0})}{\pi^{N}_{1}(\hat{\theta}_{1})}\Big{)}\cdot
c$

where $\hat{I}_{i}$ and $\hat{\theta}_{i}$ are the observed information matrix
and m.l.e, respectively under model $M_{i}$. The BIC approximation for
$\underline{B}_{01}^{I}$ gives us

$b\cdot\frac{f_{0}(\textbf{y}|\hat{\theta}_{0})}{f_{1}(\textbf{y}|\hat{\theta}_{1})}\cdot\big{(}n\big{)}^{(k_{1}-k_{0})/2}\leq\underline{B}_{01}^{I}\leq
c\cdot\frac{f_{0}(\textbf{y}|\hat{\theta}_{0})}{f_{1}(\textbf{y}|\hat{\theta}_{1})}\cdot\big{(}n\big{)}^{(k_{1}-k_{0})/2}$

In the nested hypothesis case under the null hypothesis, the likelihood ratio
in the equation above is bounded because $-2\cdot\log(LR_{01})$ is bounded (it
converges to a central Chi-square distribution by Wilks theorem). Hence,
$\underline{B}_{01}^{I}(\textbf{y})\to\infty$ and
$\overline{B}_{10}^{I}(\textbf{y})\to 0\mbox{ as }n\to\infty$ and by
definition, $0\leq B_{01}^{L}(p)\leq 1$. Hence,

$0\leq\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}\leq\frac{1}{\underline{B}_{01}^{I}}=\overline{B}_{10}^{I}\longrightarrow\lim_{n\to\infty}0=0\leq\lim_{n\to\infty}\frac{\underline{B}_{01}^{L}(p)}{\underline{B}_{01}^{I}}\leq\lim_{n\to\infty}\overline{B}_{10}^{I}=0$

By the Squeezing Theorem,

$\lim_{n\to\infty}\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}=0.\hskip
30.0pt$

∎

Although this proof is for nested hypotheses, we believe that this is true in
general and we support this claim with examples for separated hypotheses (see
Geometric vs. Poisson, Negative Binomial vs Poisson examples).

A study of the p-value and the RBL vs. the SS Bayes factor under $H_{1}$
We want to study the following limit under $H_{1}$

$\lim_{n\to\infty}\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}$

Using an approach similar to the one used in the previous theorem, we can see
that

$\frac{f_{1}(\textbf{y}|\hat{\theta_{1}})}{c\cdot
n^{k_{1}-k_{0}}f_{0}(\textbf{y}|\hat{\theta_{0}})}\leq\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}\leq\frac{f_{1}(\textbf{y}|\hat{\theta_{1}})}{b\cdot
n^{k_{1}-k_{0}}f_{0}(\textbf{y}|\hat{\theta_{0}})}$

To explore this limit we need the asymptotic distribution of the likelihood
ratio under the alternative hypothesis and related results by Self, S. G., &
Liang, K. Y. (1987) [12] known as the generalized likelihood ratio test (GLRT)
can be used to handle the case where the alternative hypothesis is true. The
GLRT is a modification of the likelihood ratio test that allows for testing of
both null and alternative hypotheses, and it is based on the same asymptotic
theory as Wilks’ theorem. Specifically, as the sample size goes to infinity,
the distribution of the GLRT under the alternative hypothesis approaches a
non-central chi-squared distribution with a non-centrality parameter
proportional to the squared distance between the true parameter value and the
null parameter value, and degrees of freedom equal to the difference in the
number of parameters between the null and alternative hypotheses.
Mathematically, we can write this as:

$LR\sim\chi^{2}(p,\lambda)\mbox{ as }n\to\infty\mbox{ under }H_{1}$

where $\chi^{2}(p,\lambda)$ denotes the non-central chi-squared distribution
with p degrees of freedom, the non-centrality parameter $\lambda$ and
$LR=2*\log(L(\hat{\theta}_{1})/L(\hat{\theta}_{0}))$. The GLRT can be used to
test hypotheses in a wide range of statistical models, including linear
regression, logistic regression, and ANOVA, among others. So, while Wilks’
theorem applies only to the null hypothesis, the related GLRT provides a way
to handle testing of the alternative hypothesis using the same underlying
asymptotic theory.

In some cases, it may be possible to write the asymptotic distribution of the
likelihood ratio test statistic as a function of n.

However, note that this formula applies only to specific cases, and the
asymptotic distribution of the likelihood ratio test statistic may be
different for other models and hypotheses. Determining the non-centrality
parameter $\lambda$ for the asymptotic distribution of the likelihood ratio
test statistic can be challenging and requires knowledge of the true parameter
values under the alternative hypothesis. In general, there is no simple
formula that can be used to calculate $\lambda$ for any hypotheses test. In
some cases, simulation methods may be used to estimate the non-centrality
parameter $\lambda$. This involves generating simulated data sets from the
assumed distribution under the alternative hypothesis, calculating the
likelihood ratio test statistic for each data set, and then estimating the
distribution of the test statistic based on the simulated values. The
estimated distribution can then be compared to the theoretical distribution to
estimate the non-centrality parameter.

In general, determining the non-centrality parameter and its rate of
convergence for the asymptotic distribution of the likelihood ratio test
statistic is an active area of research, and the methods used to estimate
$\lambda$ depend on the specific model and hypothesis being tested. We
conjecture that if there exist $r>0$ such that
$\lim_{n\to\infty}\frac{LR}{n^{r}}$ is constant under $H_{1}$, then

$\lim_{n\to\infty}\frac{B_{01}^{L}(p)}{\underline{B}_{01}^{I}}=0\mbox{ if
}k_{1}-k_{0}>r$ (3.6)

because if this limit is zero, then similarly to the case when $H_{0}$ is
true, we can show that the limit of the ratio is also zero. Further research
is needed to determine if Theorem 3.1 is true under the alternative hypothesis
in general, or if this theorem only holds under certain conditions.

##  Chapter 4 The SP Bayes Factors

In this chapter, we introduce a novel prior known as the SP Prior and leverage
it to establish what we term the SP Bayes Factor. The designation ’SP’ denotes
the ’superior posterior’ approach, which underpins the creation of these
priors. The underlying concept involves transforming an improper prior into a
proper posterior distribution by identifying either the theoretical or
empirical minimal training sample that maximizes evidence for the null
hypothesis (or yields the least favorable prior for the alternative).

Both the prior and the Bayes factor are presented in two versions-empirical
and theoretical-similar to the universal lower bound discussed in the
preceding chapter. These priors exhibit a close relationship with IBF bounds,
which we elucidate within this chapter.

The empirical prior proves valuable in cases where theoretical maximization
poses challenges. Remarkably, the empirical SP Bayes Factor closely
approximates the actual Bayes factor value compared to the theoretical
counterpart. Nevertheless, in scenarios where theoretical maximization remains
viable, the theoretical prior offers computational efficiency over the
empirical one. Employing both priors, the SP Bayes factors converge towards
accurate decisions with increasing sample sizes-a point we substantiate both
theoretically and through numerical examples in subsequent sections.

###  4.1 The SP Prior

######  Definition 4 (The Theoretical SP Prior).

The Theoretical SP prior for $B_{10}$ is defined as

$\pi^{SP}_{k}(\theta_{k})=\pi_{k}(\theta_{1}|\textbf{y}^{*}_{T}(\ell))=\frac{f_{k}(\textbf{y}^{*}_{T}(\ell)|\theta_{k})\pi_{k}^{N}(\theta_{k})}{m_{k}(\textbf{y}^{*}_{T}(\ell))},k=0,1$

$\mbox{where }\textbf{y}^{*}_{T}(\ell)=\arg[\sup_{\textbf{y}(\ell)\in
D}B_{01}(\textbf{y}(\ell))].$

######  Definition 5 (The Empirical SP Prior).

The Empirical SP prior for $B_{10}$ is defined as

$\pi_{k}^{SP*}(\theta_{k})=\pi_{k}(\theta_{k}|\textbf{y}_{E}(\ell))=\frac{f_{k}(\textbf{y}^{*}_{E}(\ell)|\theta_{k})\pi_{k}^{N}(\theta_{k})}{m_{k}(\textbf{y}^{*}_{E}(\ell))},k=0,1$

$\mbox{where }\textbf{y}^{*}_{E}(\ell)=\arg[\sup_{\textbf{y}(\ell)\in
D_{n,k}}\frac{m_{0}(\textbf{y}(\ell))}{m_{1}(\textbf{y}(\ell))}].$

###  4.2 The SP Bayes Factor

######  Definition 6 (The Theoretical SP Bayes Factor).

If $\displaystyle\arg\sup_{\textbf{y}(\ell)\in
D}B_{01}(\textbf{y}(\ell))<\infty$, then

$B^{SP}_{10}=\frac{m_{1}^{SP}(\textbf{y})}{m_{0}^{SP}(\textbf{y})}=\frac{\displaystyle\int
f_{1}(\textbf{y}|\theta_{1})\pi_{1}^{SP}(\theta_{1})d\theta_{1}}{\displaystyle\int
f_{0}(\textbf{y}|\theta_{0})\pi_{0}^{SP}(\theta_{0})d\theta_{0}}.$

######  Definition 7 (The Empirical SP Bayes Factor).

Let $\textbf{y}^{*}_{T}(\ell)=\displaystyle\arg\sup_{\textbf{y}(\ell)\in
D_{n,k}}B_{01}(\textbf{y}(\ell))$, and
$\textbf{y}(-\ell)=\textbf{y}\setminus\textbf{y}^{*}_{T}(\ell)$, then

$B^{SP*}_{10}=\frac{m_{1}^{SP*}(\textbf{y})}{m_{0}^{SP*}(\textbf{y})}=\frac{\displaystyle\int
f_{1}(\textbf{y}(-\ell)|\theta_{1})\pi_{1}^{SP*}(\theta_{1})d\theta_{1}}{\displaystyle\int
f_{0}(\textbf{y}(-\ell)|\theta_{0})\pi_{0}^{SP*}(\theta_{0})d\theta_{0}}.$

###  4.3 Results

######  Theorem 4.3.1.

The Empirical IBF upper and lower bounds can be obtained using the Empirical
SP Prior, that is

$B^{SP*}_{10}=B_{10}^{N}[\textbf{y}]\sup_{\ell=1,...,L}B_{01}{[\textbf{y}(\ell)]}=\overline{B}_{10}^{I*}(\textbf{y}).$

###### Proof.

Let $\textbf{y}(-\ell)=\textbf{y}\setminus\\{\textbf{y}^{*}_{E}(\ell)\\}$,
without loss of generality, we will prove it for the upper-bound with the
supremum since the proof for the infimum lower bound is identical. The
Empirical SP Bayes factor is

$B^{SP*}_{10}=\frac{m_{1}^{SP*}(\textbf{y})}{m_{0}^{SP*}(\textbf{y})}=\frac{\displaystyle\int
f_{1}(\textbf{y}(-\ell)|\boldsymbol{\theta}_{1})\frac{f(\textbf{y}^{*}_{E}(\ell)|\boldsymbol{\theta}_{1})\pi_{i}^{N}(\boldsymbol{\theta}_{1})}{m_{1}(\textbf{y}^{*}_{E}(\ell))}d\boldsymbol{\boldsymbol{\theta}}_{1}}{\displaystyle\int
f_{0}(\textbf{y}(-\ell)|\boldsymbol{\theta}_{0})\frac{f(\textbf{y}^{*}_{E}(\ell)|\boldsymbol{\theta}_{0})\pi_{j}^{N}(\boldsymbol{\theta}_{0})}{m_{0}(\textbf{y}^{*}_{E}(\ell))}d\boldsymbol{\theta}_{0}}$

The marginal distributions $m_{0}$ and $m_{1}$ do not depend on
$\boldsymbol{\theta}_{0}$ and $\boldsymbol{\theta}_{1}$, hence we can take
those terms out of the integrals to obtain

$\frac{m_{0}(\textbf{y}^{*}_{E}(\ell))}{m_{1}(\textbf{y}^{*}_{E}(\ell))}\frac{\displaystyle\int
f_{1}(\textbf{y}(-\ell)|\theta_{1})f(\textbf{y}^{*}_{E}(\ell)|\theta_{1})\pi_{i}^{N}(\theta_{1})d\theta_{1}}{\displaystyle\int
f_{0}(\textbf{y}(-\ell)|\theta_{0})f(\textbf{y}^{*}_{E}(\ell)|\theta_{0})\pi_{0}^{N}(\theta_{0})d\theta_{0}}=\frac{m_{0}(\textbf{y}^{*}_{E}(\ell))}{m_{1}(\textbf{y}^{*}_{E}(\ell))}\frac{m_{1}(\textbf{y}(-\ell),\textbf{y}^{*}_{E}(\ell))}{m_{0}(\textbf{y}(-\ell),\textbf{y}^{*}_{E}(\ell))}$
$=B_{01}[\textbf{y}^{*}_{E}(\ell)]B^{N}_{10}[\textbf{y}(-\ell),\textbf{y}^{*}_{E}(\ell)]$

By definition, we have that

$B_{01}[\textbf{y}^{*}_{E}(\ell)]=B_{01}[\arg(\sup_{\ell=1,...,L}B_{01}(\textbf{y}(\ell)))]=\sup_{\ell=1,...,L}B_{01}[\textbf{y}(\ell)]$

and $\\{\textbf{y}^{*}_{E}(\ell)\\}\cup\\{\textbf{y}(-\ell)\\}=\textbf{y}$ and
$m_{k}(\textbf{y}(-\ell),\textbf{y}^{*}_{E}(\ell))=m_{k}(\textbf{y}),$ for
$k=0,1$

Therefore,

$B^{SP*}_{10}=B^{N}_{10}[\textbf{y}]\sup_{\ell=1,...,L}B^{N}_{01}[\textbf{y}(\ell)]=\overline{B}_{10}^{I*}(\textbf{y})$

∎

######  Theorem 4.3.2.

Let $f(\textbf{y}|\boldsymbol{\theta})$ be a p.d.f or p.m.f a random variable
$Y$ with support $S=\\{y:f(y|\theta)>0\\}$. Suppose we want to calculate the
Bayes factor for the following nested hypothesis test
$H_{0}:\theta=\theta_{0}\mbox{ vs }H_{1}:\theta\neq\theta_{0}$, then the
empirical SP-Prior and SP-BF converges to the Theoretical SP-Prior and SP-BF.

###### Proof.

Without loss of generality we can assume that $dim(\theta)=1$, so that the
minimal training sample size is one. In this case, the set
$D_{n,k}=D_{n,1}=\textbf{y}_{n}=(y_{1},...,y_{n})\subseteq S$,
$\displaystyle\lim_{n\to\infty}D_{n,1}=D$ and $L={n\choose k}={n\choose 1}=n$.
The proof will be divided in the discrete and continuous support cases.

Case 1: The support $S$ is discrete

$\lim_{n\to\infty}D_{n,1}=\lim_{L\to\infty}\\{\textbf{y}(\ell)\\}_{l=1}^{L}=\\{\textbf{y}(\ell)\\}_{l=1}^{\infty}=D=S$

Therefore,

$\displaystyle\lim_{n\to\infty}B_{01}^{SP*}(\textbf{y})=B^{N}_{01}[\textbf{y}]\displaystyle\lim_{L\to\infty}\sup_{\\{\textbf{y}(\ell)\\}_{\ell=1}^{L}}B_{10}[\textbf{y}(\ell)]=B^{N}_{01}[\textbf{y}]\sup_{D}B_{10}[\textbf{y}(\ell)]$

Hence,

$\lim_{n\to\infty}B_{01}^{SP*}(\textbf{y})=B_{01}^{SP}(\textbf{y})$

Clearly, in the discrete support case, we can eventually (for some n) obtain
the minimal training sample $\textbf{y}_{n}^{*}(\ell)$ that generates the
theoretical supremum, and the empirical SS-BF will be exactly equal to the
theoretical SSBF.

Case 2: The support $S$ is continuous
In the case of a continuous support $S$, we can only find a sequence of
minimal training samples $\\{\textbf{y}_{n}^{*}(\ell)\\}$ that gets
arbitrarily close to a minimal imaginary sample $\textbf{y}^{*}(\ell)$ that
produces the theoretical supremum. We will construct such a sequence and prove
its existence. Assume that
$\textbf{y}^{*}(\ell)=\arg\sup_{D}B_{10}(\textbf{y}(\ell))$ and for each $n$,

$\textbf{y}_{n}^{*}(\ell)=\textbf{y}(\ell_{i_{n}}),\mbox{ where
}i_{n}=\arg\min_{i\in\\{1,...,{n\choose
1}\\}}|\textbf{y}^{*}(\ell)-\textbf{y}(\ell_{i})|$

We claim that $\textbf{y}_{n}^{*}(\ell)\to\textbf{y}^{*}(\ell)$ for any
$\textbf{y}^{*}(\ell)\in D$
Proof:
Assume that $\textbf{y}^{*}(\ell)\in D$ and $\exists\epsilon>0$ such that

$\forall
n\in\mathbb{N},|\textbf{y}^{*}(\ell)-\textbf{y}_{n}^{*}(\ell)|>\epsilon$

Then the interval
$(\textbf{y}^{*}(\ell)-\frac{\epsilon}{2},\textbf{y}^{*}(\ell)+\frac{\epsilon}{2})\cap
D=\emptyset$ which is a contradiction since $\textbf{y}^{*}(\ell)\in D$, thus
the intersection can not be empty. Similarly, if such a point does exist in
$D$, then

$\int^{\textbf{y}^{*}(\ell)+\frac{\epsilon}{2}}_{\textbf{y}^{*}(\ell)-\frac{\epsilon}{2}}f(\textbf{y}|\theta)d\textbf{y}=0\Longrightarrow\mbox{
f is not continuous in
}\textbf{y}^{*}(\ell)\Longrightarrow\textbf{y}^{*}(\ell)\notin D$

∎

##  Chapter 5 Hypothesis Testing with the Normal Distribution

This chapter explores hypothesis testing for mean and precision (the inverse
of scale) parameters within the Normal Distribution. It covers essential
concepts such as the universal lower bound, the SP/Intrinsic Priors, and the
SP/Intrinsic Bayes factor. Through numerical experiments and visual charts, we
compare these approaches and visualize the results. By blending theoretical
foundations with practical applications, this chapter aims to provide readers
with a robust understanding of conducting and interpreting these crucial
statistical tests, leveraging our novel theoretical framework..

###  5.1 The Normal Scale Hypothesis Test

Assume $\mathbf{y}=\\{y_{1},y_{2},...,y_{n}\\}$ are i.i.d samples where
$y_{i}\sim N(y_{i}|\mu,h)$ for every $i$, and that the mean $\mu$ is unknown.
We wish to decide if the precision parameter $h=1/\sigma^{2}$ is a specified
value. If we want to make a decision, we should perform the following
hypothesis test

$H_{0}:h=h_{0}\mbox{ vs }H_{1}:h\neq h_{0},\ $

The likelihood is
$f(\textbf{y}|\mu,h)=(\frac{h}{2\pi})^{n/2}\exp[-\frac{h}{2}(S^{2}+n(\bar{\textbf{y}}-\mu)^{2}$,
where $S^{2}=\sum_{i=1}^{n}(y_{i}-\bar{\textbf{y}})^{2}$. As a first step, use
the so-called "Jeffrey’s rule" equating the prior to the square root of the
determinant of the Fisher Information. In the motivating example this leads
to: $\pi^{J}(\mu,h)=C_{0}\cdot\frac{d\mu dh}{\sqrt{h}}$. For the so-called
"independence" Jeffrey’s we have instead

$\pi_{I}^{J}(\mu,h)=\pi^{J}(\mu)\pi^{J}(h)=C_{1}\cdot\frac{d\mu dh}{h}$

Putting together both priors $\pi^{J}(\mu,h)=C/h^{r}$ with $r=1/2,1$ for the
dependent or independent Jeffreys prior respectively. It turns out that the
Bayes Factor is,

$B^{N}_{01}(\mathbf{y})=\frac{h_{0}^{(n-1)/2}\exp(-h_{0}S^{2}/2)}{(2/S^{2})^{(n-r)/2}\Gamma((n-r)/2)}\times\frac{C_{0}}{C_{1}},$
(5.1)

where $r$ is equal to $1/2$ for the dependent prior, and for the independent
prior equal to one. In this thesis, we assume the independent prior, as put
forward by Jeffreys.It is apparent from this expression that the Bayes Factor
is undetermined since it depends on the undefined ratio $C_{0}/C_{1}$, the
undefined constants that come from the improper priors. It can be argued, from
several points of view, that the constants regarding the location $\mu$ cancel
out, leaving only the indeterminacy related to the hypothesis parameter $h$.
Some of these points of view are among others:

  1. 1.

Mean and Precision parameters are orthogonal in the Fisher Information Matrix,
and thus "cancel-out" in the ratio of marginal likelihoods

  2. 2.

It turns out that two location models are predictively matched (see Pericchi
(2005) [11]) and when the scale is integrated into the denominator, models
under both hypotheses become location models, and the corrections cancel out.

####  5.1.1 The universal lower bound

It is evident from the previous example that, unfortunately, the Bayes factor
remains defined only up to the constant that pertains to the precision
parameter under test. This is why Jeffrey’s suggested conventional proper
priors for the extra-parameters under the larger hypothesis (but improper for
common parameters like $\mu$ here). There have been suggested techniques
around this problem, as in Smith and Spiegelhalter (1980) [14] and (1982)
[15]. These techniques, although approximate, are useful in devising sensible
scaling for Bayes Factors. In fact, we proved there for the first time that
this approach has a very direct relationship with more recent approaches that
have been studied in detail, particularly the Intrinsic Bayes Factor, the
Intrinsic Priors, and EP-priors.

In the example of the section, the minimal training sample consists of two
observations $(y(\ell_{1}),y(\ell_{2}))$, and the expression (2) turns out to
be

$B^{N}_{01}(\mathbf{y})\times
B^{N}_{10}(\mathbf{\textbf{y}(\ell)})=\frac{h_{0}^{(n-1)/2}\exp(-h_{0}S^{2}/2)}{(2/S^{2})^{p}\frac{C_{0}}{C_{1}}\sqrt{h_{0}/\pi}\cdot\exp(-h_{0}\textbf{D}(\ell)^{2}/4)|\textbf{D}(\ell)|\cdot\frac{C_{1}}{C_{0}}},$

where $\textbf{D}(\ell)={y(\ell_{1})-y(\ell_{2})}$ and
$p=(n-1)/2\Gamma((n-1)/2)$. It is apparent from (3) that the undefined
constants cancel out. As a specific example, take $H_{0}:h_{0}=1$. The crucial
step is: what to do about the (theoretical or imaginary) training samples
summary statistics $\textbf{D}(\ell)$? SSBF’s practical and simplifying
approach is to take:

$\sup_{\textbf{y}(\ell)\in D}B^{N}_{10}(\textbf{D}(\ell)),$ (5.2)

which is attained at $\hat{D}=\sqrt{2}$ and the "correction" (4) is equal to
0.484. The $IBF$ bound is then
$\overline{B}_{01}^{I}=B^{N}_{01}\times\sup_{\textbf{y}(\ell)\in
D}B^{N}_{10}(\textbf{y}(\ell)).$ To see the relationship with more recent
approaches like Intrinsic Bayes Factors the $\sup$ above is replaced by the
arithmetic mean or by the theoretical expectation, on $\textbf{y}(\ell)$,
under $H_{1}\supset H_{0}$. The theoretical expressions in (4) and (5), which
take the formal maximum possible value for all possible training samples, can
be refined by the empirical observed maximum for the observed
$\textbf{y}(\ell)$’s:

$\sup_{\ell=1,...,L}B^{N}_{10}(\mathbf{y(\ell)})\leq\sup_{\textbf{y}(\ell)\in
D}B^{N}_{10}(\textbf{y}(\ell)),$ (5.3)

where the distinction arises as $\mathbf{y(\ell)}$ is the set of empirically
observed training samples and D which is the set of all possible minimal
training samples from the whole support. The empirical bound is then sharper
than the theoretical bound, and we can refer to the theoretical bound and the
empirical bound of the Intrinsic Bayes Factors. In the Expected IBF case, for
any pair of i.i.d. training samples $\\{y(\ell_{1}),y(\ell_{2})\\}$, assumed
(when taking expectations) to be drawn from $H_{1}\supset H_{0}$, results in
$\textbf{D}(\ell)\sim N(0,2/h)$. Then under $H_{1}$:

$\frac{B^{EIBF}_{10}(\textbf{y})}{B^{N}_{10}(\mathbf{y})}=\int\frac{|\textbf{D}|}{\sqrt{\pi}}exp(-\frac{\textbf{D}^{2}}{4})\cdot\frac{\sqrt{h}\cdot
exp(-h\cdot\textbf{D}^{2}/4)}{2\sqrt{\pi}}d\textbf{D}=\frac{2\sqrt{h}}{\pi(h+1)}.$

In an applied problem with real data, the parameter $h$ will be estimated by
its MLE. In this theoretical exercise, we study it as a function of $h$, and
computation yields that its maximum is attained at the null hypothesis
$h=h_{0}=1$ and its value at its maximum is $1/\pi=0.318$.

####  5.1.2 Remarks

  1. 1.

Our ratio of constants is 0.484 and the Maximum of the Expected IBF of the
Training Sample is 0.318. This may seem like a big difference (the ratio is
about .66). These values are attained close to the null hypothesis, the
difference far from the null, may be much bigger but there, the dominant
factor $B^{N}_{01}(\mathbf{y})$, will overwhelm that difference for a moderate
sample, and even small but larger than the training sample sizes.

  2. 2.

The lower bound is robust, it does not depend on the fact that all the
training samples belong to the alternative hypothesis $H_{1}$ and that all
belong to the same population. This bound achieves a very important robust
effect. The bound has Universal validity: It is still valid under violations
of distributional assumptions, for example, Normality or homogeneity of the
population of training samples. It can also be used to check the computations
of IBFs and with Intrinsic Priors.

  3. 3.

The SS’s bound on Intrinsic Bayes Factors is typically easy to compute and
generalize to realistic scenarios.

  4. 4.

In terms of p-values, if it can be found a "tight" bound $g(p)$ so that
$g(p)\leq B^{N}_{01}$ then

$g(p)\times\frac{1}{\sup
B^{N}_{01}(\textbf{y}(\ell))}\leq\underline{B}^{I}_{01}\leq B_{01}$ (5.4)

which is a Universal Bound that adapts to the sample size n.

####  5.1.3 The SP and Intrinsic Priors

In this section, we turn our attention to the calculation of Superior
Posterior (SP) and Intrinsic priors within the context of the Normal Scale
Hypothesis testing scenario. An in-depth exploration will be conducted to
illuminate the nuanced attributes and similarities shared between these two
types of priors. To enhance the comprehension and tangibility of the
underlying principles, we shall include illustrative representations in the
form of charts. These graphical elucidations will serve as effective tools for
visualizing the properties and similarities of both SP and Intrinsic priors.

Moreover, to ensure results and facilitate the application of the presented
methods, we will accompany our discussion with executable R code snippets.
These will illustrate the computational procedures required to carry out the
prior calculations and generate the corresponding charts. By making this code
available, we aim to provide a practical and accessible means for readers to
engage with our methodology. Consider the Bayes factor for the minimal
training sample of the normal scale example with an unknown mean

$B_{10}^{N}(\textbf{y}(\ell))\propto\frac{1}{\exp\\{-h_{0}\textbf{D}(\ell)^{2}/4\\}|\textbf{D}(\ell)|}=\frac{\exp\\{h_{0}\textbf{D}(\ell)^{2}/4\\}}{|\textbf{D}(\ell)|}$

To find the value of $\textbf{D}(\ell)$ that generates the $\sup$ of
$B_{10}^{N}(\textbf{y}(\ell))$, we must take logarithms on both sides, and
differentiate the equation $B_{10}^{N}(\textbf{y}(\ell))=0$ with respect to
$\textbf{D}(\ell)$,

$\frac{d}{d\textbf{D}(\ell)}\log(\frac{\exp\\{h_{0}\textbf{D}(\ell)^{2}/4\\}}{|\textbf{D}(\ell)|})=2h_{0}\frac{\hat{D}}{4}-\frac{1}{\hat{D}}=0\Longrightarrow
h_{0}\frac{\hat{D}^{2}}{2}=1\Longrightarrow\boxed{\hat{D}=\sqrt{\frac{2}{h_{0}}}}$

After calculating this value, we must consider the prior of $\mu$ and $h$
given $\textbf{y}(\ell)$ which is

$\pi(\mu,h|\textbf{y}(\ell))=\frac{f(\textbf{y}(\ell)|\mu,h)\pi(\mu)\pi(h)}{\int
f(\textbf{y}(\ell)|\mu,h)\pi(\mu)\pi(h)d\mu
dh}=\frac{\frac{1}{h}\frac{h}{2\pi}\exp\\{-\frac{h}{2}[(y(1)-\mu)^{2}+(y(2)-\mu)^{2}]\\}}{\frac{1}{2\pi}\displaystyle\int\exp\\{-\frac{h}{2}(\frac{\textbf{D}(\ell)^{2}}{2}+2(\mu-\bar{\textbf{y}}(\ell))^{2}\\}d\mu
dh}$ (5.5)

Simplifying this equation, we obtain

$\frac{\exp\\{-\frac{h}{2}[\frac{\textbf{D}(\ell)^{2}}{2}+2(\mu-\bar{\textbf{y}}(\ell))^{2}]\\}}{\displaystyle\int^{\infty}_{0}\exp\\{-\frac{h}{4}\textbf{D}(\ell)^{2}\\}\frac{\sqrt{2\pi}}{\sqrt{2h}}dh}$

To solve the last integral, we should consider the following change of
variable

$u=\frac{h\textbf{D}(\ell)^{2}}{4},\hskip
10.0pth=\frac{4u}{\textbf{D}(\ell)^{2}},\hskip
10.0ptdh=\frac{4}{\textbf{D}(\ell)^{2}}du$
$\frac{4\sqrt{\pi}}{\textbf{D}(\ell)^{2}}\displaystyle\int_{0}^{\infty}(\frac{4u}{\textbf{D}(\ell)^{2}})^{\frac{1}{2}-1}\exp\\{-u\\}du=\Gamma(\frac{1}{2})\frac{4\sqrt{\pi}}{\textbf{D}^{2}}\frac{D}{2}=\frac{2\pi}{\textbf{D}}$

Hence,

$\pi(\mu,h|D=\hat{D})=\frac{\exp\\{-\frac{h}{2}[\frac{\hat{D}^{2}}{2}+2(\mu-\bar{y}(\ell))^{2}]\\}\hat{D}}{2\pi}$
(5.6)

Integrating with respect to $\mu$, we obtain the prior for $h$

$\displaystyle\pi(h|\textbf{D})$
$\displaystyle=\int\pi(\mu,h|\textbf{D})\,d\mu=\int\frac{\textbf{D}}{2\pi}\exp\left\\{-\frac{h}{2}\left(\frac{\textbf{D}^{2}}{2}+2(\mu-\bar{\textbf{y}}(\ell))^{2}\right)\right\\}\,d\mu$
$\displaystyle=\frac{\textbf{D}}{2\pi}\exp\left\\{-\frac{h}{2}\frac{\textbf{D}^{2}}{2}\right\\}\frac{\sqrt{2\pi}}{\sqrt{2h}}=\frac{\textbf{D}\exp\left\\{-\frac{h\textbf{D}^{2}}{4}\right\\}}{\sqrt{2\pi}\sqrt{2h}}$

and by substituting D by $\hat{D}$, we can obtain the SP prior for $h$

$\boxed{\pi^{SP}(h|\hat{D}=\sqrt{\frac{2}{h_{0}}})=\frac{1}{\sqrt{2\pi
hh_{0}}}\exp\\{-\frac{h}{2h_{0}}\\}}$

Evidently, the prior $\pi^{SP}(h)$ is $Gamma(\alpha=1/2,\beta=2h_{0})$. To
find the prior for $\mu$, we must integrate with respect to $h$ as follows

$\pi^{SP}(\mu|\textbf{D})=\displaystyle\int_{0}^{\infty}\exp\\{-\frac{h}{2}[\frac{\textbf{D}^{2}}{2}+2(\mu-\bar{\textbf{y}}(\ell))^{2}]\\}\frac{\textbf{D}}{2\pi}dh=\frac{8\pi}{\textbf{D}^{2}+4(\mu-\bar{\textbf{y}}(\ell))}$

$\pi^{SP}(\mu|\hat{D}=\sqrt{\frac{2}{h_{0}}})=\frac{8\pi}{\frac{2}{h_{0}}+4(\mu-\bar{\textbf{y}}(\ell))^{2}}=\frac{2\pi}{\frac{1}{2h_{0}}+(\mu-\bar{y}(\ell))^{2}}=\frac{4h_{0}\pi}{1+(\frac{\mu-\bar{\textbf{y}}(\ell)}{\sqrt{2h_{0}}})^{2}}$

The distribution of $\mu$ for $B_{01}$ is proportional to a
$Cauchy(\bar{\textbf{y}},\gamma)$ where $\gamma=\sqrt{2h_{0}}$.

$\boxed{\pi^{SP}(\mu)\propto\frac{1}{\pi\gamma\Big{(}1+(\frac{\mu-\bar{\textbf{y}}(\ell)}{\gamma})^{2}\Big{)}}}$

Next, we will compute the intrinsic priors for the same hypothesis test above
to compare those with the SP priors. The improper intrinsic prior for $\mu$ is
proportional to a constant $\pi^{I}(\mu)\propto c$. In the other hand, the
intrinsic prior for the scale $h$ is

$\pi^{I}(\mu,h)=\pi^{N}(\mu,h)\displaystyle\int\frac{m_{0}^{N}(\textbf{y}(\ell))}{m_{1}^{N}(\textbf{y}(\ell))}f(\textbf{y}(\ell)|\mu,h)d\textbf{y}(\ell)=\displaystyle\int\pi(\mu,h|\textbf{y}(\ell))m_{0}(\textbf{y}(\ell))d\textbf{y}(\ell)$
$\boxed{\pi^{I}(h)\propto\frac{1}{\pi}\frac{\sqrt{\frac{h}{h_{0}}}}{\frac{h+1}{h_{0}}}\frac{1}{h}=\frac{1}{\pi
h_{0}}\frac{(\frac{h}{h_{0}})^{1/2-1}}{(\frac{h}{h_{0}}+1)}}$

Evidently, the prior $\pi^{I}(h)$, is a $SBeta2(h|p=1/2,q=1/2,b=h_{0})$ (Scale
beta2 distribution). The Scale beta two density is given by

$f(y|p,q,b)=\frac{\Gamma(p+q)}{\Gamma(p)\Gamma(q)}\frac{(y/b)^{p-1}}{(\frac{y}{b}+1)^{(p+q)}},\hskip
5.0pty,p,q,b>0$

The subsequent script is written in the R programming language that
facilitates the generation of a comparative graphical representation of both
priors. These priors are superimposed on the same set of axes to foster a
direct comparison. This visual comparison is executed across distinct values
of the null hypothesis ($H_{0}=1$ and $H_{0}=10$).
R program

⬇

1 sp_prior = function(h,h0) { (1/sqrt(2*pi*h*h0))*exp(-h/(2*h0)) }

2 i_prior = function(h,h0) { (1/(pi*h0))*( (h/h0)^(-0.5) )/(h/h0+1) }

3 par(mfrow=c(1,2)); h = seq(0.1,5,0.01);h0=1; iprior = i_prior(h,h0); spprior
= sp_prior(h,h0)

4 ymin = min(spprior,iprior); ymax = max(spprior,iprior)

5 plot(h,spprior,type="l",main="Comparison of the SP and Intrinsic Priors of
h",ylab="Priors")

6 mtext(text="Hypothesis test: H0 = 1 vs H0 != 1",side=3);
lines(h,iprior,type="l",col="red")

7 legend(3.4, 1.2, legend=c("SP", "Intrinsic"),col=c("black", "red"),
lty=c(1,1), cex=0.8,title="Priors", text.font=4, bg=’lightblue’)

8 h = seq(1,15,0.01);h0=10; iprior = i_prior(h,h0); spprior = sp_prior(h,h0)

9 ymin = min(spprior,iprior); ymax = max(spprior,iprior)

10 plot(h,spprior,type="l",main="Comparison of the SP and Intrinsic Priors of
h",ylim=c(ymin,ymax),ylab="Priors")

11 mtext(text="Hypothesis test: H0 = 10 vs H0 != 10",side=3);
lines(h,iprior,type="l",col="red")

12 legend(10.5, 0.12, legend=c("SP", "Intrinsic"),col=c("black", "red"),
lty=c(1,1), cex=0.8,title="Priors", text.font=4, bg=’lightblue’)

Figure 5.1: A comparison between the SP Prior and the Intrinsic Prior for the
Normal Scale Hypothesis test for different null hypothesis values.

As illustrated in the preceding figures, a noteworthy characteristic of the
Superior Posterior (SP) Prior emerges. This prior assigns a heightened
probability to the null hypothesis in comparison to the Intrinsic Prior,
specifically when the values of $h$ are less than that of the null hypothesis.
It is furthermore observed that the probability ascribed by the SP is
invariably equal to or exceeds that of the Intrinsic Prior under these
conditions. This pattern elucidates the differential weights that these priors
place on the null hypothesis, offering valuable insights for our
investigation.

###  5.2 The Normal Mean Hypothesis Test

This section is dedicated to exploring the domain of Normal Mean Hypothesis
Testing. Here, we uncover the core principles and methodologies that drive
inference regarding the mean parameter within the Normal Distribution. Assume
that $y_{i}\sim N(\mu,\sigma_{0}^{2})$ and consider the following hypothesis
test

$H_{0}:\mu=\mu_{0}\mbox{ vs }H_{1}:\mu\neq\mu_{0},\sigma_{0},\mbox{ known }$

Hence, the non-informative prior is $\pi^{N}(\mu)=1$ and a minimal training
sample is only one observation $y(\ell)$, for any $\ell=1,...,n$. We will use
that information to calculate different Bayes factors and their bounds in the
next sections.

####  5.2.1 The Expected Intrinsic Bayes Factor for the Normal Mean Hypothesis
Test

Consider the general formula for the Expected Intrinsic Bayes Factor (EIBF)

$\boxed{B_{10}^{EI}(\textbf{y})=B_{10}^{N}(\textbf{y})E^{\hat{\theta}}_{\textbf{y}(\ell)|H_{0}}(\frac{m_{0}^{N}(\textbf{y}(\ell))}{m_{1}^{N}(\textbf{y}(\ell))})}$

For this specific hypothesis test, the EIBF can be obtained by solving the
following integral

$E^{\hat{\theta}}_{\textbf{y}(\ell)|H_{0}}(\frac{m_{0}^{N}(\textbf{y}(\ell))}{m_{1}^{N}(\textbf{y}(\ell))})=\int_{-\infty}^{\infty}\Big{(}\frac{1}{\sqrt{2\pi}\sigma_{0}}\Big{)}^{2}\exp\\{-\frac{1}{2\sigma_{0}^{2}}(\textbf{y}(\ell)-\mu_{0})^{2}\\}\exp\\{-\frac{1}{2\sigma_{0}^{2}}(\textbf{y}(\ell)-\mu)^{2}\\}d\textbf{y}(\ell)$

This integral can be easily solved to obtain the following equality

$E^{\hat{\theta}}_{\textbf{y}(\ell)|H_{0}}(\frac{m_{0}^{N}(\textbf{y}(\ell))}{m_{1}^{N}(\textbf{y}(\ell))})=N(\mu|\mu_{0},2\sigma_{0}^{2})=\frac{1}{\sqrt{2\pi}\sqrt{2\sigma_{0}^{2}}}\exp\\{\frac{-1}{4\sigma_{0}^{2}}(\mu-\mu_{0})^{2}\\}$

In this case, the parameter $\mu$ is unknown and we should estimate it using
the maximum likelihood estimator $\bar{\textbf{y}}$ to obtain

$E^{\hat{\theta}}_{\textbf{y}(\ell)|H_{0}}(\frac{m_{0}^{N}(\textbf{y}(\ell))}{m_{1}^{N}(\textbf{y}(\ell))})=\frac{1}{\sqrt{2\pi}\sqrt{2\sigma_{0}^{2}}}\exp\\{\frac{-1}{4\sigma_{0}^{2}}(\bar{\textbf{y}}-\mu_{0})^{2}\\}$

After finding this term, we can proceed to calculate the Expected Intrinsic
Bayes Factor

$B_{10}^{EI}=\frac{1}{\sqrt{2\pi}\sqrt{2\sigma_{0}^{2}}}\exp\\{\frac{-1}{4\sigma_{0}^{2}}(\bar{\textbf{y}}-\mu_{0})^{2}\\}\frac{\sqrt{2\pi}\frac{\sigma_{0}}{\sqrt{n}}}{\exp\\{\frac{-1}{2\sigma_{0}^{2}}n(\bar{\textbf{y}}-\mu_{0})^{2}\\}}$

Simplifying the equation above yields

$\boxed{B_{10}^{EI}=\frac{1}{\sqrt{2n}}\exp\\{\frac{1}{2\sigma_{0}^{2}}(n-\frac{1}{2})(\bar{\textbf{y}}-\mu_{0})^{2}\\}}$

####  5.2.2 The Intrinsic Bayes Factor Bounds for the Normal Mean Hypothesis
Test

Let us consider a scenario where the random vector
$\mathbf{Y}=(Y_{1},\ldots,Y_{n})$ is composed of independent and identically
distributed (i.i.d) random variables $Y_{i}$, following a normal distribution
with mean $\mu_{0}$ and standard deviation $\sigma_{0}$ (known) under the null
hypothesis $M_{0}$, and a normal distribution with mean $\mu$ and standard
deviation $\sigma_{0}$ (known) under the alternative hypothesis $M_{1}$. We
adopt the non-informative prior, specifically $\pi_{1}^{N}(\mu)={1}$.

$B_{01}^{N}(\textbf{y}(\ell))=\frac{\frac{1}{\sqrt{2\pi}\sigma_{0}}\exp\\{\frac{-1}{2\sigma_{0}^{2}}(\textbf{y}(\ell)-\mu_{0})^{2}\\}}{\displaystyle\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi}\sigma_{0}}\exp\\{\frac{-1}{2\sigma_{0}^{2}}(\textbf{y}(\ell)-\mu)^{2}\\}d\mu}$

From the equation above, it is evident that the maximum value of the Bayes
factor is obtained when $\textbf{y}(\ell)=\mu_{0}$, and the maximum value is

$\boxed{\sup_{\textbf{y}(\ell)\in
D}B_{01}^{N}(\textbf{y}(\ell))=\frac{1}{\sqrt{2\pi}\sigma_{0}}}.$

By using the result above, we can compute the Bayes factor bound

$\displaystyle\overline{B}_{10}^{I}$
$\displaystyle=B_{10}^{N}(\textbf{y})\times\sup B_{01}^{N}(\textbf{y}(\ell))$
$\displaystyle=\frac{\displaystyle\int_{0}^{\infty}(\frac{1}{\sqrt{2\pi}\sigma_{0}})^{n}\exp\\{-\frac{1}{2\sigma_{0}^{2}}\sum[(\textbf{y}(\ell)-\bar{\textbf{y}})^{2}+n(\bar{\textbf{y}}-\mu)^{2}]\\}d\mu}{\displaystyle(\frac{1}{\sqrt{2\pi}\sigma_{0}})^{n}\exp\\{-\frac{1}{2\sigma_{0}^{2}}\sum(\textbf{y}(\ell)-\mu_{0})^{2}\\}}\times\sup
B_{01}^{N}(\textbf{y}(\ell))$
$\displaystyle=\frac{\sqrt{2\pi}\frac{\sigma_{0}}{\sqrt{n}}}{\exp\\{\frac{-1}{2\sigma_{0}^{2}}n(\bar{\textbf{y}}-\mu_{0})^{2}\\}}\sup
B_{01}^{N}(y(\ell))=\frac{\sqrt{2\pi}\frac{\sigma_{0}}{\sqrt{n}}}{\exp\\{\frac{-1}{2\sigma_{0}^{2}}n(\bar{\textbf{y}}-\mu_{0})^{2}\\}}\frac{1}{\sqrt{2\pi}\sigma_{0}}$

Simplifying the equation above yields

$\boxed{\overline{B}_{10}^{I}=\frac{\exp\\{\frac{1}{2\sigma_{0}^{2}}n(\bar{\textbf{y}}-\mu_{0})^{2}\\}}{\sqrt{n}}}$

Similarly, we can obtain the IBF lower bound

$\boxed{\underline{B}_{01}^{I}={\sqrt{n}}{\Big{(}\exp-\\{\frac{1}{2\sigma_{0}^{2}}n(\bar{\textbf{y}}-\mu_{0})^{2}\\}}\Big{)}}.$

####  5.2.3 Normal mean hypothesis test with unknown variance.

Let us consider a scenario where the random vector
$\mathbf{Y}=(Y_{1},\ldots,Y_{n})$ is composed of independent and identically
distributed (i.i.d) random variables $Y_{i}$, following a normal distribution
with mean $\mu_{0}$ and standard deviation $\sigma_{0}$ under the null
hypothesis $M_{0}$, and a normal distribution with mean $\mu$ and standard
deviation $\sigma_{1}$ under the alternative hypothesis $M_{1}$.

$H_{0}:\mu=\mu_{0}\mbox{ vs }H_{1}:\mu\neq\mu_{0},\sigma\mbox{ unknown}$

We adopt non-informative priors, specifically
$\pi_{0}^{N}(\sigma_{0})={1}/{\sigma_{0}}$ and
$\pi_{1}^{N}(\mu,\sigma_{1})={1}/{\sigma_{1}^{2}}$. It is important to note
that the choice of $\pi_{1}^{N}$ differs from that of $\pi_{0}^{N}$ due to a
seminal work by Berger, Pericchi, and Varshavsky in 1998 [5], wherein they
calculated the marginal of $\textbf{y}(\ell)$ using
$\pi_{1}^{N}(\mu,\sigma_{1})={1}/{\sigma_{1}}$. In this context, the
expression for the density function under the alternative hypothesis takes the
form:

$m_{1}(\textbf{y}(\ell))=\frac{1}{2|y_{1}-y_{2}|}$

where $\mathbf{y}(\ell)$ represents the observed data. However, we face a
challenge in evaluating the Bayes Factor due to the nature of the above
density function. Specifically, the infimum and supremum of this function can
potentially become either zero or infinity, leading to indeterminate bounds
for the Bayes Factor. To circumvent this issue and ensure well-defined bounds,
we opt for the aforementioned choice of priors. This strategic choice of
priors enables us to derive meaningful marginals and subsequently evaluate the
Bayes Factor.

$m_{0}^{N}(\textbf{y}(\ell))=\frac{1}{2\pi(y_{1}^{2}+y_{2}^{2})},\hskip
10.0ptm_{1}^{N}(\textbf{y}(\ell))=\frac{1}{\sqrt{\pi}(y_{1}-y_{2})^{2}}$

and the Bayes factor is

$B_{01}(\textbf{y}(\ell))=\frac{m_{0}(\textbf{y}(\ell))}{m_{1}(\textbf{y}(\ell))}=\frac{(y_{1}-y_{2})^{2}}{2\sqrt{\pi}(y_{1}^{2}+y_{2}^{2})}$

To find the IBF bounds, we found the supremum of the Bayes factor above

$\sup_{\textbf{y}(\ell)\in
D}B_{01}(\textbf{y}(\ell))=\frac{m_{0}(\textbf{y}(\ell))}{m_{1}(\textbf{y}(\ell))}=\frac{(y_{1}-y_{2})^{2}}{2\sqrt{\pi}(y_{1}^{2}+y_{2}^{2})}=\frac{1}{\sqrt{\pi}}$

The IBF upper bound can be obtained as

$\overline{B}_{10}=B_{10}(\textbf{y})\sup
B_{01}(\textbf{y}(\ell))=\frac{m_{1}(\textbf{y})}{m_{0}(\textbf{y})}\frac{1}{\sqrt{\pi}}$

To find the upper bound, we must calculate the marginals $m_{0}$ and $m_{1}$.
To obtain $m_{0}$, we must solve the following integral

$m_{0}(\textbf{y}))=\int_{0}^{\infty}(\frac{1}{\sqrt{2\pi}})^{n}\exp\\{\frac{-1}{2\sigma_{0}^{2}}[\sum(y_{i}-\mu_{0})^{2}]\\}\frac{1}{\sigma_{0}^{(n+1)}}d\sigma_{0}$

Consider the change of variable
$\beta=\frac{1}{2}[\sum(y_{i}-\mu_{0})^{2}],\frac{1}{\sigma^{2}}=\eta$, then

$\frac{-2}{\sigma^{3}}d\sigma=d\eta\Longrightarrow
d\sigma=-\frac{\sigma^{3}}{2}d\eta$
$\frac{1}{\eta^{1/2}}=\sigma\Longrightarrow d\sigma=-\eta^{3/2}/2d\eta$
$\boxed{m_{0}(\textbf{y})=\frac{1}{2}(\frac{1}{\sqrt{2\pi}})^{n}\int_{0}^{\infty}\exp\\{-\beta\eta\\}{\eta^{n/2-1}}d\eta=\frac{1}{2}(\frac{1}{\sqrt{2\pi}})^{n}\frac{\Gamma(n/2)}{\beta^{n/2}}}$

Similarly, we can calculate $m_{1}$

$m_{1}(\textbf{y})=(\frac{1}{\sqrt{2\pi}})^{n}\int_{0}^{\infty}\frac{1}{\sigma_{1}^{n+2}}\int_{-\infty}^{\infty}\exp\\{-\frac{1}{2\sigma_{1}^{2}}(\sum(y_{i}-\mu)^{2}))\\}d\mu
d\sigma_{1}$

Adding and subtracting $\bar{\textbf{y}}$ inside the exponential function
yields

$\sum(y_{i}-\mu)^{2}=\sum((y_{i}+\bar{\textbf{y}})-(\mu+\bar{\textbf{y}}))^{2}=\sum(y_{i}-\bar{\textbf{y}})^{2}+n(\mu-\bar{\textbf{y}})^{2}$

Inserting this result inside the integral with respect to $\mu$ results in

$\displaystyle\int_{-\infty}^{\infty}\exp\\{-\frac{1}{2\sigma_{1}^{2}}(\sum(y_{i}-\mu)^{2}))\\}d\mu$
$\displaystyle=\exp\\{\frac{-1}{2\sigma_{1}^{2}}\sum(y_{i}-\bar{y})^{2}\\}\int_{-\infty}^{\infty}\exp\\{\frac{-1}{2(\sigma_{1}/\sqrt{n})^{2}}(\mu-\bar{\textbf{y}})^{2}\\}d\mu$
$\displaystyle=\sigma_{1}\sqrt{2\pi/n}\exp\\{-\frac{1}{2\sigma_{1}^{2}}(\sum(y_{i}-\bar{\textbf{y}})^{2})\\}$

To calculate $m_{1}$, we insert this result to solve the integral with respect
to $\sigma_{1}$

$\frac{1}{\sqrt{n}}(\frac{1}{\sqrt{2\pi}})^{n-1}\int_{0}^{\infty}\frac{1}{\sigma_{1}^{n+1}}\exp\\{-\frac{1}{2\sigma_{1}^{2}}(\sum(y_{i}-\bar{y})^{2})\\}d\sigma_{1}$

Let $\beta=\frac{\sum(y_{i}-\bar{y})^{2}}{2}$ and
$\eta=\frac{1}{\sigma_{1}^{2}}$, then

$\frac{1}{2\sqrt{n}}(\frac{1}{\sqrt{2\pi}})^{n-1}\int_{0}^{\infty}\eta^{n/2-1}\exp\\{-\beta\eta\\}d\eta$

This integral is the kernel of a Gamma distribution which results in

$\boxed{m_{1}(\textbf{y})=\frac{1}{2\sqrt{n}}(\frac{1}{\sqrt{2\pi}})^{n-1}\frac{\Gamma(n/2)}{(\frac{\sum(y_{i}-\bar{y})^{2}}{2})^{n/2}}}$

After obtaining $m_{0}$ and $m_{1}$, we can calculate the Bayes Factor upper
bound

$\boxed{\overline{B}_{10}=\sqrt{2/n}(\frac{\sum(y_{i}-\mu_{0})^{2}}{\sum(y_{i}-\bar{y})^{2}})^{n/2}}.$

Similarly, the Bayes Factor Lower bound is

$\boxed{\underline{B}_{01}=\sqrt{n/2}(\frac{\sum(y_{i}-\bar{y})^{2}}{\sum(y_{i}-\mu_{0})^{2}})^{n/2}.}$

##  Chapter 6 The Exponential Model Hypothesis Test

In the exploration of statistical hypothesis testing, this chapter embarks on
a comprehensive investigation into exponential hypothesis testing.
Specifically, we delve into the comparison of two hypotheses
$H_{0}:\lambda=\lambda_{0}$ and $H_{1}:\lambda\neq\lambda_{0}$ where $\lambda$
is the parameter for an exponential density $f(y|\lambda)=\lambda e^{-\lambda
y}$. We will employ several robust methodologies to asses this hypothesis test
from different but consistent approaches.

Our focus centers on evaluating and contrasting several pivotal methodologies:
the EP Approach, the $-ep\log(p)$ lower bound, the Universal Robust Bound for
the Intrinsic Bayes Factor, and the SP Approach. These methodologies serve as
our guiding compasses in the pursuit of discerning their effectiveness in
hypothesis validation. We rely on simulations to shed light on the performance
and consistency of these methodologies. By employing these simulations, we aim
to craft illustrative charts that will serve as visual aids, elucidating the
comparative efficacy and nuances of each approach. This chapter unfolds as a
meticulous exploration, aiming not only to identify the most effective
methodology but also to unravel the intricacies underlying their application.

The exponential density is defined as

$f(y|\lambda)={\lambda}e^{-\lambda y},\mbox{ where }y\geq 0\mbox{ and
}\lambda>0$

The likelihood, Jeffrey’s prior, and posterior distribution for an exponential
density with parameter $\lambda$ are given by

$L(\lambda|\textbf{y})={\lambda^{n}}\exp(\displaystyle-\lambda\sum_{i=1}^{n}y_{i}),\hskip
5.0pt\pi^{N}(\lambda)\propto\frac{1}{\lambda},\hskip
5.0pt\pi(\lambda|\textbf{y})=\frac{L(\lambda|\textbf{y})\pi_{J}(\lambda)}{m(\textbf{y})}$

After calculating the likelihood and Jeffrey’s prior for the exponential
density we can easily obtain the following marginal

$m(\textbf{y})\propto\int_{0}^{\infty}{\lambda^{(n-1)}}\exp(-\displaystyle{\lambda\sum_{i=1}^{n}y_{i}})d\lambda=\displaystyle\frac{(\sum_{i=1}^{n}y_{i})^{n}}{\Gamma(n)}$

Finally, we can obtain the posterior distribution for $\lambda$

$\boxed{\pi(\lambda|\textbf{y})\propto\frac{\displaystyle{\Gamma(n)}{\lambda^{n-1}}\exp\Big{(}-{\lambda}{\sum_{i=1}^{n}y_{i}}\Big{)}}{{(\sum_{i=1}^{n}y_{i})^{n}}}}$

In the following sections, we’ll tackle a hypothesis test:

$H_{0}:\lambda=\lambda_{0}\mbox{ vs. }H_{1}:\lambda\neq\lambda_{0}$

We’ll use different methods to calculate Bayes Factors, bounds, and priors.
This helps us compare traditional approaches with our new methods, giving us a
clearer understanding of their outcomes.

###  6.1 The Intrinsic Bayes Factor Dynamic and Universal Lower Bounds

In this section, we’re computing our new bound for the Intrinsic Bayes Factor.
To derive this bound, we start by examining the Bayes factor for the minimal
training sample:

$B_{01}^{N}(y(\ell))=\frac{\mathbf{y(\ell)}\lambda_{0}}{\exp^{\lambda_{0}\mathbf{y(\ell)}}}$

Additionally, we need the full Bayes factor for this test

$B_{01}^{N}(\textbf{y})=\frac{\lambda_{0}^{n}\exp(-\lambda_{0}\sum
y_{i})}{\int_{0}^{\infty}{\lambda^{n}\exp(-\lambda\sum
y_{i})}\frac{1}{\lambda}d\lambda}=\frac{\lambda_{0}^{n}\exp(-\lambda_{0}\sum
y_{i})(\sum y_{i})^{n}}{\Gamma(n)}$

By minimizing the first Bayes factor, we obtain the Empirical IBF lower bound

$\boxed{\underline{B}^{I^{*}}_{01}=\frac{(\sum
y_{i})^{n}\lambda_{0}^{n-1}}{\textbf{y}^{*}(\ell)\Gamma(n)}\exp(\lambda_{0}[y^{*}(\ell)-\sum
y_{i}]),\hskip
5.0pty^{*}(\ell)=\arg\displaystyle\min_{y_{i}\in\\{y_{1},...,y_{n}\\}}\Big{\\{}\frac{\mathbf{y_{i}}\lambda_{0}}{\exp^{\lambda_{0}\mathbf{y_{i}}}}\Big{\\}}}$

Our next goal is to obtain the theoretical IBF bound, and to achieve this we
need to find $\displaystyle\sup_{y(\ell)\in D}B^{N}_{01}(y(\ell))$. This
maximization can be attained by solving the following equation

$\displaystyle\frac{d}{dy(\ell)}[\log(B_{01}^{N}(y(\ell))]$
$\displaystyle=\frac{d}{dy(\ell)}[\log(y(\ell))+\log(\lambda_{0})-\lambda_{0}y(\ell)]$
$\displaystyle=0$ $\displaystyle=\frac{1}{y^{*}(\ell)}+\lambda_{0}$
$\displaystyle=0$ $\displaystyle y^{*}(\ell)=\frac{1}{\lambda_{0}}$

To verify that this is a maximum, we take the second derivative with respect
to $y(\ell)$

$\frac{d^{2}}{dy(\ell)}[\log(B_{01}^{N}(y(\ell))]=\frac{-1}{y(\ell)^{2}}<0,\forall\mbox{
}y(\ell)\in\mathbb{R}$

Therefore, $y^{*}(\ell)$ the supremum in $(0,\infty)$ and
$B_{01}^{N}(y^{*}(\ell))=1/e$. By using these results, we can obtain the
following Intrinsic Bayes Factor upper bound

$\boxed{\overline{B}^{I}_{10}=1/e\frac{\Gamma(n)\exp(\lambda_{0}\sum
y_{i})}{(\sum y_{i})^{n}\lambda_{0}^{n}}}.$

Similarly, the IBF lower bound is given by

$\boxed{\underline{B}_{01}^{I}=\frac{e(\sum
y_{i})^{n}\lambda_{0}^{n}}{\Gamma(n)\exp(\lambda_{0}\sum y_{i})}}$

###  6.2 The SP Bayes Factors

In this section, our focus lies in computing the SP priors and the SP Bayes
factor for the exponential hypothesis test. We aim to contrast the outcomes
from the SP approach with the intrinsic Bayes factor bounds derived earlier.
To bolster our findings, we’ll conduct simulations and generate illustrative
charts. These visual aids will fortify and elucidate our conclusions.

The theoretical SP Priors for the exponential model for $B_{10}$ are given by

$\boxed{\pi_{0}^{SP}(\lambda)=1,\hskip
5.0pt\pi_{1}^{SP}(\lambda)=\frac{1}{\lambda_{0}}\exp^{-\lambda/\lambda_{0}},\hskip
5.0pt\lambda>0}$

It’s noteworthy that $\pi_{1}^{SS}(\lambda)$ follows an exponential
distribution with the parameter $1/\lambda_{0}$. Let’s assume, without loss of
generality, that $y_{1}=1/{\lambda_{0}}$. The SP marginals for the hypothesis
test $H_{0}:\lambda=\lambda_{0}\text{ versus }H_{1}:\lambda\neq\lambda_{0}$
are

$\displaystyle m^{SP}_{0}(\textbf{y})$
$\displaystyle\propto\lambda_{0}^{n}e^{-\lambda_{0}\sum_{i=1}^{n}y_{i}}$
$\displaystyle=\lambda_{0}^{n}e^{-\lambda_{0}\sum_{i=1}^{n}y_{i}}$
$\displaystyle m^{SP}_{1}(\textbf{y})$
$\displaystyle\propto\frac{1}{\lambda_{0}}\int_{0}^{\infty}\lambda^{n}e^{-\lambda(\sum_{i=1}^{n}y_{i}+\frac{1}{\lambda_{0}})}$
$\displaystyle=\frac{1}{\lambda_{0}}[\Gamma(n+1)(\sum_{i=1}^{n}y_{i}+\frac{1}{\lambda_{0}})^{(-(n+1))}]$
$\displaystyle=\frac{1}{\lambda_{0}}[\Gamma(n+1)(\sum_{i=1}^{n}y_{i}+\frac{1}{\lambda_{0}})^{-(n+1)}]$

Utilizing these marginals leads to the results of the SP Bayes factor as
follows:

$\boxed{B^{SP}_{10}=\frac{\Gamma(n+1)e^{\lambda_{0}\sum_{i=1}^{n}y_{i}}}{(\sum_{i=1}^{n}y_{i}+\frac{1}{\lambda_{0}})^{n+1}\lambda_{0}^{n+1}}}$

In order verify the validity of our calculations, and to compare the numerical
values of $B^{SP}_{10}$ with $\overline{B}^{I}_{10}$, a simulation was
conducted and the corresponding chart was obtained. The simulation was
implemented using the following code snippet:

R Program

⬇

1 b10_theoretical = function(n,h0,lambda=lambda0) {

2 y = c(); b10_t = c(); b10_sp = c(); s = c(); e = exp(1)

3 for (i in seq(1,n,1)) {

4 y = c(y,rexp(1,rate=1/lambda))

5 s = sum(y)

6 val = -log(e,e) + - n*log(s,e) - n*log(h0,e) + log(gamma(n),e)+ h0*s

7 b10_t = c(b10_t,val)

8 val = log(gamma(n+1),e) + h0*s - (n+1)*log(s + 1/h0,e) - (n+1)*log(h0,e)

9 b10_sp = c(b10_sp,val)

10 }

11 plot(seq(1,n,1),b10_t,type="l",col="red",

12 xlab="Number of samples",

13 ylab="log(BF10)", ylim=c(min(b10_t,b10_sp),max(b10_t,b10_sp)),

14 main="Theorical SS Bayes Factor (Red) vs SP Bayes Factor (Blue) (H0 TRUE)")

15 lines(b10_sp,col="blue")

16 }

Figure 6.1: $\overline{B}^{SS}_{10}$ vs $B^{SP}_{10}$ under $H_{0}$ for
increasing sample size and exponential data.

This comparison serves to analyze the numerical relationship between
${B}^{SP}_{10}$ and $\overline{B}^{I}_{10}$, thereby contributing to the
evaluation and understanding of their respective characteristics. The obtained
chart provides a visual representation of the convergence of the
$\overline{B}^{I}_{10}$ to ${B}_{10}^{SP}$ as the number of samples increases.

###  6.3 The EP Bayes Factors

The exponential Expected posterior prior under the alternative (EP-Prior) for
the exponential hypothesis test $H_{0}:\lambda=\lambda_{0}\mbox{ vs
}H_{1}:\lambda\neq\lambda_{0}$ can be calculated as

$\pi_{1}^{EP}(\lambda)=\frac{\displaystyle\sum_{i=1}^{L}\pi(y(\ell)|\theta_{i})}{n}=\frac{\displaystyle\sum_{i=1}^{n}y_{i}e^{-\lambda
y_{i}}}{n}$

It can be verified that the prior in question integrates unity with respect to
$\lambda$ by performing integration overall positive values of $\lambda$.
Additionally, the SP prior, which is an exponential distribution with a
parameter of $1/{\lambda_{0}}$, also integrates one.

The objective of this study is to compare the SP prior with the EP prior under
different sample sizes, specifically for the cases of $n=10$ and $n=100$. To
achieve this, simulations were conducted utilizing the following R code:
R Program

⬇

1 prior = function(data,lambda0,lambda) {

2 data = c(); object = c(); lambda0 = 0; ep1 = c(); ss1 =
rep(0,length(lambda)); ss1 = (1/lambda0)*exp(-lambda/lambda0)

3 for (j in 1:length(lambda)) {

4 epnew = 0

5 for (i in 1:length(data)) { epnew = epnew + data_[i]*exp(-data[i]*lambda[j])
}

6 ep1 = c(ep1,epnew)

7 }

8 object$ss = ss1; object$ep = ep1/length(data)

9 return(object)

10 }

11 data = rexp(n=100,r=1);lambda0=1; lambda = seq(0.1,2,0.1)

12 object = prior(data = data,lambda0=lambda0,lambda=lambda);

13 ss = exp(-lambda); ep = object$ep

14
plot(lambda,ep,col="red",ylim=c(min(ss,ep),max(ss,ep)),type="b",ylab="Prior",

15 main="EP Prior under H1 (Red) / SP Prior under H1 (Blue)";

16 lines(lambda,ss,col="blue")

By employing these simulations, we aim to investigate and evaluate the
performance of the SP prior relative to the EP prior in the context of various
sample sizes. The obtained results will contribute to enhancing our
understanding and analysis of the behavior and characteristics of these
priors.

(a) Exponential simulation with 10 samples for $\lambda_{0}=1$

(b) Exponential simulation with 100 samples for $\lambda_{0}=1$

Figure 6.2: Simulations with different number of samples for the exponential
SP and EP Priors for $\lambda_{0}=1$

The exponential EP-Bayes factor for $\lambda=\lambda_{0}\mbox{ vs
}\lambda\neq\lambda_{0}$ can be obtained by calculating the following
expression

$B_{10}^{EP}(\textbf{y})=\frac{m_{1}^{EP}(\textbf{y})}{m_{0}^{EP}(\textbf{y})}=\frac{\frac{1}{n}\displaystyle\int_{0}^{\infty}\lambda^{n}\exp\\{-\lambda\sum_{i=1}^{n}y_{i}\\}\sum_{i=1}^{n}y_{i}\exp\\{-\lambda
y_{i}\\}d\lambda}{\lambda_{0}^{n}\exp\\{-\lambda_{0}\displaystyle\sum_{i=1}^{n}y_{i}\\}}$

Our objective is to resolve the integral mentioned earlier. To facilitate
this, let’s define $S=\displaystyle\sum_{i=1}^{n}y_{i}$. Subsequently, the
integral transforms into:

$\displaystyle\sum_{i=1}^{n}y_{i}\int_{0}^{\infty}\lambda^{n}\exp\\{-\lambda(S+y_{i})\\}d\lambda$

Consider the change of variable $\alpha=n+1$ and
$\beta_{i}=\frac{1}{S+y_{i}}$, which leads to the transformation:

$\displaystyle\sum_{i=1}^{n}y_{i}\int_{0}^{\infty}\lambda^{n}\exp\\{-\lambda(S+y_{i})\\}d\lambda=\displaystyle\sum_{i=1}^{n}y_{i}\int_{0}^{{\infty}}\lambda^{\alpha-1}\exp\\{-\lambda/\beta_{i}\\}d\lambda$

This new form represents the kernel of a Gamma density with parameters
$\alpha$ and $\beta_{i}$, resulting in the integral:

$\Gamma(n+1)\sum_{i=1}^{n}y_{i}\Big{(}\frac{1}{S+y_{i}}\Big{)}^{n+1}$

Hence, the EP Bayes Factor, in this case, becomes:

$\boxed{B_{10}^{EP}(\textbf{y})=\frac{\frac{1}{n}\Gamma(n+1)\displaystyle\sum_{i=1}^{n}y_{i}\Big{(}\frac{1}{S+y_{i}}\Big{)}^{n+1}}{\lambda_{0}^{n}\exp\\{-\lambda_{0}\displaystyle\sum_{i=1}^{n}y_{i}\\}}}$

The principal objective of this numerical experiment is to compare the EP and
SP Bayes factors, previously derived, by creating 100 simulations, each
consisting of 100 samples. The process entails averaging the Bayes factor for
each sample size, ranging from 1 to 100. Furthermore, we intend to plot the
logarithm of each Bayes factor on the same axis to discern their respective
behaviors under two distinct scenarios: when the null hypothesis is valid and
when it is not. The forthcoming R program has been devised to generate the
anticipated output.

R Program

⬇

1 #Initialize vectors and variables

2 finallogbf10 = c(); finallogbf10sp=c(); nsim = 100; n = 100;

3 data = c(); s = 0; sumdata = 0;logbf10=c();logbf10sp=c();

4 #H0 TRUE

5 true_lambda = 1; lambda0 = 1;

6 #Use this for H0 FALSE

7 #true_lambda = 1; lambda0 = 0.3;

8 for (k in 1:nsim) {

9 #Generate random samples

10 data =rexp(n=n,rate=true_lambda)

11 for (i in 1:n) {

12 s = sum(data[1:i])

13 for (j in 1:i) { sumdata = sumdata + data[j]*((1/(s+data[j]))^(i+1)) }

14 #EP marginals

15 m1_ep = (gamma(i+1)*sumdata)/i; m0_ep = (lambda0)^i*exp(-lambda0*s)

16 #Log EP Bayes factor

17 logbf10 = c(logbf10,log(m1_ep) - log(m0_ep))

18 #SP Marginals

19 m1_sp = gamma(i+1)*exp(lambda0*s); m0_sp = ((s + 1/lambda0)^(i+1)
)*(lambda0^(i+1))

20 #Log SP Bayes Factor

21 logbf10sp = c(logbf10sp,log(m1_sp) - log(m0_sp))

22 sumdata = 0

23 }

24 #Storing simulations

25 finallogbf10 = cbind(finallogbf10,logbf10); finallogbf10sp =
cbind(finallogbf10sp,logbf10sp)

26 }

27 logbf10_avg = rep(0,n); logbf10sp_avg = rep(0,n);

28 #Averaging Bayes factors for all the simulations

29 for (i in 1:n) {

30 logbf10_avg[i] = mean(finallogbf10[i,]); logbf10sp_avg[i] =
mean(finallogbf10sp[i,])

31 }

32 x = 1:n

33 y1 = logbf10_avg[x]; y2 = logbf10sp_avg[x]

34 if (lambda0 == true_lambda) { h0 =TRUE } else { h0 = FALSE }

35 plot(x,y1,type="l",col="red",ylim=c(min(y1,y2),max(y1,y2)), main="log(EP
Bayes Factor 10) (Red) vs log(SP Bayes Factor 10) (Black)",xlab="Sample
size",ylab="Log Bayes Factor 10")

36 lines(x,y2,type="l",col="black")

37 mtext(paste("H0 is ",h0," / Number of simulations: ",nsim," / Number of
samples: ",n," / lambda = ",true_lambda, " / lambda_0 = ",lambda0,
sep=""),side=3)

Figure 6.3: Logarithm of the SP and EP Bayes factor using one-hundred
exponential simulations (100 samples each) for $\lambda_{0}=1,H_{0}$ True

Figure 6.3 demonstrates the consistency of both the SP and EP Bayes factors.
The figure illustrates that both Bayes factors approach zero when calculating
the Bayes factor of the alternative hypothesis over the null hypothesis. This
observation is expected since the data used for the calculations are random
samples generated from the null hypothesis. The decreasing trend of the Bayes
factors indicates their convergence towards zero, indicating consistency in
the evidence provided by both methods.

Figure 6.4: Logarithm of the SP and EP Bayes factor using one-hundred
exponential simulations (100 samples each) for $\lambda_{0}=1,H_{0}$ False

Figure 6.4 demonstrates the consistency of both the SP and EP Bayes factors.
The figure illustrates that both Bayes factors increase in value with the
sample size when calculating the Bayes factor of the alternative hypothesis
over the null hypothesis. This observation is expected since the data used for
the calculations are random samples generated from the alternative hypothesis.
The increasing trend of the Bayes factors indicates consistency in the
evidence provided by both methods.

###  6.4 The -eplogp lower bound

In this section we revisit the Sellke, Bayarri, and Berger (2001) [13] lower
bound for the Bayes factor that employs the classical p-value of the
hypothesis test, represented as $-ep\log(p)\leq B_{01}$. We will perform a
simulation to compare this this p-value approximation of the Bayes factor to
contrast it with the intrinsic Bayes factor lower bound for the exponential
test. In an attempt to compute the $-ep\log(p)$ lower bound of $B_{01}$, we
determined the likelihood ratio for the hypotheses $H_{0}$ against $H_{1}$.
The approximation of the p-value was accomplished through the utilization of
Wilk’s theorem and the chi-squared distribution. The graphical representation
provided demonstrates the lower bound and the $-ep\log(p)$ in tandem as the
sample size increases.

The forthcoming illustration elucidates how both the empirical and theoretical
bounds converge toward the appropriate decision as the sample size escalates.
This is demonstrated by drawing samples from an exponential distribution
characterized by a known rate parameter, while simultaneously altering the
null hypothesis.

R Program

⬇

1

2 exponential_simulations=function(n_sim,n_data,h0,rate) {

3 lambda = h0; sim = c(); bfbound = c(); bfssbound = c(); pval_seq= c();
epbound = c(); ssbound = c();

4 for (j in 1:n_sim) {

5 for (n in 1:n_data) {

6 sim = c(sim,rexp(n=1,rate=rate))

7 lrt_null = n*log(lambda) - lambda*sum(sim)

8 lrt_alt = n*(log(n) - log(sum(sim)) - 1)

9 lrt = exp(lrt_null-lrt_alt)

10 ts = -2*log(lrt)

11 #p-value approximation wilks theorem

12 pval = 1-pchisq(ts,1)

13 pval_seq = c(pval_seq,pval)

14 #selke bound

15 if (pval > (1/exp(1))) { bfbound = c(bfbound,1) }

16 else { bfbound = c(bfbound,-exp(1)*pval*log(pval)) }

17 bfssbound =
c(bfssbound,(exp(1)*(sum(sim)*lambda)^n)/gamma(n)*exp(-lambda*sum(sim)))

18 }

19 epbound = cbind(epbound,bfbound); ssbound = cbind(ssbound,bfssbound)

20 }

21 epbound_avg = c(); ssbound_avg = c()

22 for (i in 1:n_data) {

23 epbound_avg = c(epbound_avg,mean(epbound[i,]))

24 ssbound_avg = c(ssbound_avg,mean(ssbound[i,]))

25 }

26 object = c(); object$epbound_avg = epbound_avg; object$ssbound_avg =
ssbound_avg

27 return(object)

28 }

29

30 #100 Simulations and 100 Samples; H0 True

31 n_sim = 100; n_data = 100; h0 = 1; rate = 1

32

33 simulations = exponential_simulations(n_sim=n_sim,n_data =
n_data,h0=h0,rate=rate)

34

35 par(mfrow=c(1,2))

36 plot(1:n_data,simulations$epbound_avg,type="l",col="red",xlab="Number of
samples",ylab="Average of the Lower
bounds",ylim=c(min(simulations$epbound_avg,simulations$ssbound_avg),max(simulations$ssbound_avg,simulations$epbound_avg)),

37 main="Comparison between the -eplogp and IBF lower bounds")

38 lines(simulations$ssbound_avg,col="blue")

39

40 mtext(paste("Simulations: lambda = ",rate, " / Null hypothesis: lambda =
",h0, sep=""),side=3)

41 legend(0, 10, legend=c("-eplogp", "SS"),col=c("red", "blue"), lty=c(1,1),
cex=0.8,title="Lower bounds", text.font=4, bg=’lightblue’)

By employing the aforementioned code, a simulation was carried out to
calculate both the empirical and theoretical lower bounds of $B^{I}_{01}$ for
the exponential distribution. This was executed to test the hypothesis $H_{0}$
against $H_{1}$ under two distinct circumstances: when $H_{0}$ holds true and
when it does not. The objective was to examine the convergence of the Bayes
factor lower bound toward the accurate decision as the sample size expands.
The outcomes of this simulation are presented in the following figures.

Figure 6.5: The $-ep\log(p)$ and $\underline{B}_{01}^{I}$ using 100
exponential simulations (100 samples each) for $\lambda_{0}=1$ (Null
hypothesis true), and 100 exponential simulations with 10 samples each for
$\lambda_{0}=1$ and the simulations rate parameter $r=5$ (false null
hypothesis)

In Figure 6.5 we illustrated the consistency of the IBF lower bound in both
scenarios. On the other hand, the second figure (right) illustrates that both
Bayes factor’s lower bounds approach zero when calculating the Bayes factor of
the null hypothesis over the alternative hypothesis. This observation is
expected since the data used for the calculations are random samples generated
from the alternative hypothesis. The decreasing trend of the lower bounds
indicates their convergence towards zero, indicating consistency in the
evidence provided by both methods under the alternative, however, the
$-ep\log(p)$ bound does not provide decisive evidence to accept the null
hypothesis when we draw random samples from it, because the evidence is
bounded by one.

###  6.5 Conclusions

Our research indicates that the Theoretical SP Bayes Factor, derived using the
SP-Priors, is congruent with the Intrinsic Bayes factor Bound
$\overline{B}^{I}_{10}$, albeit with the inclusion of an additional
observation being the supremum of $B^{N}_{01}(y(\ell))$. We demonstrated that
the SS-prior yields comparable results and that this methodology facilitates
the generation of a lower bound for the theoretical SS Bayes factor.

As the sample size increases, the discrepancy between these two methods
diminishes, given that we are more likely to encounter samples closer to the
supremum (under $H_{0}$). In the context of comparing the EP and SP Priors for
$H_{1}$, it becomes evident that the Theoretical SP prior remains unchanged
for any sample, while the EP-prior fluctuates in accordance with the samples.

Moreover, it can be inferred from the simulations and the chart that the SP
prior assigns a higher probability to the null hypothesis $\lambda_{0}=1$,
owing to its wider area from 0 to 1 compared to the EP-prior for both $n=10$
and $n=100$. Both the EP-Bayes factor and the SP-Bayes factors exhibit similar
patterns of convergence and values. However, the SP-Bayes factor presents
higher values under the null and lower values under the alternative.

The consistencies observed in the results validate the results of the SP-Bayes
factors and confirm the accuracy of the computations performed in this study.

##  Chapter 7 Linear regression and ANOVA Hypothesis Testing with Bayes
Factors

In the realm of statistical modeling, the avenues of linear regression and
analysis of variance (ANOVA) stand as cornerstones, offering valuable insights
into relationships between variables and differences among groups. Traditional
frequentist approaches have long governed these methodologies, relying on
p-values and significance testing. However, the burgeoning field of Bayesian
statistics introduces a paradigm shift, emphasizing coherent inference and
robust evidence assessment through Bayes factors.

This chapter embarks on an exploration of linear regression and ANOVA within
the Bayesian framework. Departing from traditional methods, we delve into the
application of Bayes factor approximations, building upon the discussions from
our previous chapters. Moreover, we extend the foundations laid by Smith and
Spiegelhalter, broadening our insights by calculating these approximations
through the utilization of a generalized Jeffrey’s prior.

ANOVA, a stalwart technique in comparing group means, undergoes a
transformation in this chapter as we explore its foundations through Bayesian
inference. The comparison of multiple groups finds new light as Bayes factors
unveil nuanced insights into the evidence supporting differing hypotheses.

###  7.1 Linear regression

Linear regression is a statistical model that is commonly utilized in several
disciplines, ranging from economics to engineering. At its core, linear
regression is a conditional model where the outcome variable is predicated on
a linear combination of the predictor variables, in conjunction with an
unobserved error term that introduces variability into the relationship
between the input and output variables.

The fundamental form of a simple linear regression model is as follows:

$Y=\theta_{0}+\theta_{1}X+\varepsilon$ (7.1)

where $Y$ is the dependent or outcome variable, $X$ is the independent or
predictor variable, $\theta_{0}$ is the Y-intercept, which represents the
expected value of $Y$ when all $X$ are 0, $\theta_{1}$ is the slope of the
regression line, indicating the degree to which $Y$ changes for each unit
change in $X$, and $\varepsilon$ is the error term, embodying the difference
between the observed and predicted values of $Y$.

The goal of linear regression analysis is to estimate the coefficients
$\theta_{0}$ and $\theta_{1}$ that minimize the sum of the squared residuals,
thus providing the "best fit" line for the observed data.

In a multiple linear regression scenario, the model is expanded to include
more than one independent variable:

$Y=\theta_{0}+\theta_{1}X_{1}+\theta_{2}X_{2}+...+\theta_{p}X_{p}+\varepsilon$
(7.2)

where $p$ represents the number of predictor variables. Here, each
$\theta_{i}$ (for $i=1,2,...,p$) signifies the change in the expected value of
$Y$ for each unit change in the corresponding predictor $X_{i}$, holding all
other predictors constant.

Linear regression models assume linearity, independence, homoscedasticity
(constant variance), and normality of residuals. Violations of these
assumptions may necessitate the application of more complex models or
transformations of the data. Linear regression serves as a fundamental
building block in understanding more intricate statistical models.

###  7.2 Two nested normal-linear models

Consider the comparison of two nested normal-linear models $M_{0}\subset
M_{1}$, defined as in Smith and Spiegelhalter (1980) [14] by

$M_{i}:\textbf{y}\sim
N(\textbf{A}_{i}\theta_{i},\sigma^{2}\textbf{I}_{n}),\hskip 10.0pti=0,1$

where $\textbf{A}_{i}$ is a full rank $p_{i}$ known matrix, y is a vector with
dimension $n$, $\boldsymbol{\theta_{i}}=(\theta_{i1},...,\theta_{ip_{i}})$ is
a vector of $p_{i}$ unknown parameters and $\sigma^{2}$ is unknown. This can
be written in matrix notation as

${\displaystyle\textbf{y}=\textbf{A}_{i}\theta_{i}+\epsilon_{i}},i=0,1$

where $\boldsymbol{\epsilon}_{i}\sim N(\textbf{0},\sigma^{2}I_{n})$ for
$i=0,1$. Let

$\hat{\theta}_{i}=(A_{i}^{T}A_{i})^{-1}A_{i}^{T}\textbf{y}\mbox{ and
}\textbf{R}_{i}=|\textbf{y}-A_{i}\hat{\theta}_{i}|^{2}$

denote the least squares estimator of $\theta_{i}$ and the residual sum of
squares, respectively. The Bayes factor in this case is given by

$B_{01}=\frac{\int\int
p(\textbf{y}|\textbf{A}_{0},\boldsymbol{\theta}_{0},\sigma)p(\boldsymbol{\theta}_{0},\sigma|\textbf{A}_{0})d\boldsymbol{\theta}_{0}d\sigma}{\int\int
p(\textbf{y}|\textbf{A}_{1},\boldsymbol{\theta}_{1},\sigma)p(\boldsymbol{\theta}_{1},\sigma|\textbf{A}_{1})d\boldsymbol{\theta}_{1}d\sigma}$

Smith and Spiegelhalter showed that for the case when $i=0,1$,
$p(\boldsymbol{\theta_{i}},\sigma|\textbf{A}_{i})$ has an improper limit form
representing vague prior information for all parameters inside each linear
model. The improper limiting version of the normal-inverse-$\chi^{2}$
conjugate prior can be written in the form

$p(\boldsymbol{\theta}_{i},\sigma|\textbf{A}_{i})=p(\boldsymbol{\theta}_{i}|\textbf{A}_{i},\sigma)p(\sigma)=c_{i}(2\pi\sigma^{2})^{-p_{i}/2}\sigma^{-1}$

and the likelihood is given by

$p(\textbf{y}|A_{i},\theta_{i},\sigma)=(\frac{1}{\sqrt{2\pi}\sigma})^{n}\exp\\{-\frac{1}{2\sigma^{2}}\Big{[}\textbf{R}_{i}(\textbf{y})+(\theta_{i}-\textbf{y})^{T}A_{i}^{T}A_{i}(\theta_{i}-\textbf{y})\Big{]}\\}.$

The Bayes factor in this case is given by:

$B_{01}=\frac{c_{0}}{c_{1}}[|\textbf{A}_{1}^{T}\textbf{A}_{1}|/|\textbf{A}_{0}^{T}\textbf{A}_{0}|]^{\frac{1}{2}}\big{[}1+\frac{(p_{1}-p_{0})}{(n-p_{1})}F\big{]}^{-(n/2)}$
(7.3)

Here, $F$ represents the $F$-test statistic for comparing models $M_{0}$ and
$M_{1}$. Due to the ratio of undefined constants, this Bayes factor is
inherently indeterminate. The formula for the $F$-test statistic is:

${\displaystyle
F={\frac{\left({\frac{\textbf{R}_{0}-\textbf{R}_{1}}{p_{1}-p_{0}}}\right)}{\left({\frac{{\textbf{R}}_{1}}{n-p_{1}}}\right)}},}$

Moreover, the $F$ value used in the Bayes Factor equation can be derived from
the following identity:

$\frac{\textbf{R}_{0}}{\textbf{R}_{1}}=1+\Big{(}\frac{p_{1}-p_{0}}{n-p_{1}}\Big{)}F$

Spiegelhalter and Smith proposed a satisfactory solution to the problem of
determining the ratio $c_{0}/c_{1}$ by introducing the concept of an imaginary
training sample. Let $\textbf{A}_{0}(\ell)$ and $\textbf{A}_{1}(\ell)$ be the
design matrix of $M_{0}$ and $M_{1}$ occurring in the "thought experiment"
generating the imaginary training sample, they obtained the following result

$\frac{c_{0}}{c_{1}}=[|\textbf{A}_{1}(\ell)^{T}\textbf{A}_{1}(\ell)|/|\textbf{A}_{0}(\ell)^{T}\textbf{A}_{0}(\ell)|]^{-\frac{1}{2}}$
(7.4)

Our goal is to calculate the SS Bayes Factor, in general, using the
generalized prior $\pi(\theta_{i},\sigma)\propto{\sigma}^{-(1+q_{i})}$ where
different values of q will result in well-known priors such as the reference,
complete and modified Jeffrey’s priors. Additionally, we will calculate this
Bayes factor for ANOVA, to compare it with the results obtained by
Spiegelhalter and Smith (1980)[14].

Suppose that we want to generalize the result above and compare $M_{0}$ vs
$M_{1}$ where $M_{0}$ has $p_{0}$ unknown parameters, and $M_{1}$ has $p_{1}$
unknown parameters. Assume that our prior has the general form

$\pi(\theta_{i},\sigma|A_{i})\propto\sigma^{-(1+q_{i})},i=0,1$

where there are many possible choices for $q_{i}$. According to Berger and
Pericchi (1996) [3], $q_{i}=0$ is the reference prior (example above), and
$q_{i}=p_{1}-p_{0}$ is the modified Jeffrey’s prior. In general, the marginal
distribution for $n$ samples and $i=0,1$ is given by

$m_{i}(\textbf{y})=\displaystyle\int_{-\infty}^{\infty}\int_{0}^{\infty}\frac{c_{i}}{\sigma^{1+q_{i}}}(\frac{1}{\sqrt{2\pi}\sigma})^{n}\exp\\{-\frac{1}{2\sigma^{2}}\Big{[}\textbf{R}_{i}(\textbf{y})+(\theta_{i}-\textbf{y})^{T}A_{i}^{T}A_{i}(\theta_{i}-\textbf{y})\Big{]}\\}d\theta_{i}d\sigma$

Where
$\textbf{R}_{i}(\textbf{y})=(\textbf{y}-A_{i}\hat{\theta_{i}})^{T}(\textbf{y}-A_{k}\hat{\theta_{i}})$.
Since the second term inside the exponential function is a multivariate normal
distribution, we can integrate it with respect to $\theta_{i}$. Let
$\dim(\theta_{i})=p_{i}$, then the marginal can be written as:

$m_{i}(\textbf{y})=c_{i}\displaystyle\Big{(}\frac{1}{\sqrt{2\pi}}\Big{)}^{n-p_{i}}|\textbf{A}_{i}^{T}\textbf{A}_{i}|^{-1/2}\int_{0}^{\infty}(\frac{1}{\sigma})^{n+q_{i}-p_{i}+1}\exp\\{-\frac{\textbf{R}_{i}(\textbf{y})}{2\sigma^{2}}\\}d\sigma.$

We can solve this integral by employing the following change of variables

$\eta=\frac{\textbf{R}_{i}(\textbf{y})}{2\sigma^{2}},\hskip
10.0pt\sigma=\frac{\textbf{R}^{1/2}_{i}(\textbf{y})}{\sqrt{2}\eta^{1/2}},\hskip
10.0ptd\sigma=-\frac{\textbf{R}^{1/2}_{i}(\textbf{y})}{\sqrt{2}}\frac{1}{2}\frac{1}{\eta^{3/2}}d\eta$
$\displaystyle m_{i}(\textbf{y})$
$\displaystyle=\Big{(}\frac{1}{\sqrt{2\pi}}\Big{)}^{n-p_{i}}c_{i}|\textbf{A}_{i}^{T}\textbf{A}_{i}|^{-1/2}\int_{0}^{\infty}\exp^{-\eta}\eta^{\frac{n+q_{i}+1-p_{i}}{2}}\Big{(}\frac{2}{\textbf{R}_{i}(\textbf{y})}\Big{)}^{(n+q_{i}+1-p_{i})/2}\Big{(}\frac{\textbf{R}_{i}(\textbf{y})}{2}\Big{)}^{1/2}\frac{1}{2}\frac{1}{\eta^{3/2}}d\eta$
$\displaystyle=2^{\frac{n+q-p_{i}-2}{2}}\Big{(}\frac{1}{\sqrt{2\pi}}\Big{)}^{n-p_{i}}\Big{(}\frac{1}{\textbf{R}_{i}(\textbf{y})}\Big{)}^{\frac{n+q_{i}-p_{i}+1}{2}}\Big{(}\textbf{R}_{i}(\textbf{y})\Big{)}^{1/2}|\textbf{A}_{i}^{T}\textbf{A}_{i}|^{-1/2}\frac{c_{i}}{2}\int_{0}^{\infty}\exp^{-\eta}\eta^{\frac{n+q_{i}-p_{i}+3}{2}-1}d\eta$

Solving the integral above yields

$\int_{0}^{\infty}\exp^{-\eta}\eta^{\frac{n+q_{i}-p_{i}+3}{2}-1}d\eta=\Gamma(\frac{n+q_{i}-p_{i}+3}{2})$

Therefore, the marginal for $M_{i}$ is

$\displaystyle m_{i}(\textbf{y})$ $\displaystyle=c_{i}\mbox{
}2^{\frac{n+q_{i}-p_{i}}{2}}\Big{(}\frac{1}{\sqrt{2\pi}}\Big{)}^{n-p_{i}}\Big{(}\frac{1}{\textbf{R}_{i}(\textbf{y})}\Big{)}^{\frac{n+q_{i}-p_{i}+1}{2}}\Big{(}\frac{\textbf{R}_{i}(\textbf{y})}{2}\Big{)}^{1/2}|\textbf{A}_{i}^{T}\textbf{A}_{i}|^{-1/2}\Gamma(\frac{n-p_{i}+q_{i}+3}{2})$
$\displaystyle\propto c_{i}\sqrt{|\textbf{A}_{i}^{T}\textbf{A}_{i}|^{-1}\hskip
3.0pt[4\pi]^{p_{i}}\hskip
3.0pt\textbf{R}_{i}(\textbf{y})^{-(n+q_{i}-p_{i})}}\hskip 3.0pt$

After we obtain the marginals, we can calculate the Bayes factor as

$B_{01}^{N}(\textbf{y})=\frac{m_{0}(\textbf{y})}{m_{1}(\textbf{y})}=\frac{c_{0}}{c_{1}}\sqrt{[4\pi]^{(p_{0}-p_{1})}\frac{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}\frac{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y})]^{n+q_{0}-p_{0}}}}$

The following step involves computing the Bayes factor for the minimal
training sample. For this situation, where at least one observation per
unknown parameter is necessary, the minimal training sample size can be found
as:

$n_{01}=\max\\{p_{0}+1,p_{1}+1\\}.$

In this case, the Bayes factor for the minimal training sample takes the form:

$B_{10}^{N}(\textbf{y}(\ell))=\frac{m_{1}(\textbf{y}(\ell))}{m_{0}(\textbf{y}(\ell))}=\frac{c_{1}}{c_{0}}\sqrt{[4\pi]^{(p_{1}-p_{0})}\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}}$

The previously discussed Lemma 1 establishes that:

$B_{01}(\textbf{y}(-\ell)|\textbf{y}(\ell))=B_{01}^{N}(\textbf{y})B_{10}^{N}(\textbf{y}(\ell))$

Consequently,

$B_{01}(\textbf{y}(-\ell)|\textbf{y}(\ell))=\sqrt{\frac{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\frac{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y})]^{n+q_{0}+p_{0}}}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}}$

By leveraging these results, we derive the Empirical Intrinsic Bayes Factor
(IBF) lower bound for General Linear Models:

$\underline{B}_{01}^{GLI^{*}}=\sqrt{\frac{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}\frac{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y})]^{n+q-p_{0}}}}\displaystyle\min_{\ell=1,...,L}\sqrt{\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}}$

Additionally, the Theoretical Intrinsic Bayes Factor (IBF) Lower bound for
General Linear Models can be expressed as:

$\underline{B}_{01}^{GLI}=\sqrt{\frac{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}\frac{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y})]^{n+q_{0}-p_{0}}}\inf_{y(\ell)\in
D}\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}}$
(7.5)

###  7.3 The One-way Layout (ANOVA)

The One-way Layout in analysis of variance (ANOVA) involves $m$ groups of
observations, each with $n_{i}$ observations ($j=1,...,n_{i}$) independently
drawn from a normal distribution $y_{ij}\sim N(\mu_{i},\sigma^{2})$, given
$\mu_{1},...,\mu_{m},\sigma^{2}$.

The models considered are:

$M_{0}:\mu_{1}=...=\mu_{m}\text{ versus }M_{1}:\mu_{i}\neq\mu_{j}\text{ for
some }i\neq j.$

In this scenario, $p_{1}=m$, $p_{0}=1$, and the general form of the matrices
$\textbf{A}_{0}$ and $\textbf{A}_{1}$ can be derived. Let
$n=\sum_{i=1}^{m}n_{i}$. Notably, the design matrices $A_{0}$ and $A_{1}$ for
the one-way ANOVA can generally be written as follows:

$A_{0}=\begin{bmatrix}1\\\ \cdot\\\ \cdot\\\ \cdot\\\ 1\end{bmatrix}$,
$A_{0}^{T}=\begin{bmatrix}1&\cdot\cdot\cdot&1\\\ \end{bmatrix},$
$A_{1}=\begin{bmatrix}1&0&\cdot\cdot\cdot&0&0\\\
\vdots&\vdots&\cdot\cdot\cdot&\vdots&\vdots\\\ 1&0&\cdot\cdot\cdot&0&0\\\
0&1&\cdot\cdot\cdot&0&0\\\ \vdots&\vdots&\cdot\cdot\cdot&\vdots&\vdots\\\
0&1&\cdot\cdot&0&0\\\ \vdots&\vdots&\cdot\cdot\cdot&\vdots&\vdots\\\
0&0&\cdot\cdot\cdot&0&1\\\ \vdots&\vdots&\cdot\cdot\cdot&\vdots&\vdots\\\
0&0&\cdot\cdot\cdot&0&1\end{bmatrix}$

The matrices $A_{0}$ and $A_{1}$ are ($n\times m$) matrices, and each block of
ones in $A_{1}$ corresponds to each of the m groups. In this case, every block
i has $n_{i}$ rows and $\det(\textbf{A}_{0}^{T}\textbf{A}_{0})=n$. The matrix
product of $\textbf{A}_{1}^{T}$ and $\textbf{A}_{1}$ can be expressed as:

$\textbf{A}_{1}^{T}\textbf{A}_{1}=\begin{bmatrix}n_{1}&0&\cdot\cdot\cdot&0\\\
\cdot&n_{2}&\cdot\cdot\cdot&\cdot\\\ \cdot&\cdot&\ddots&\cdot\\\
0&0&\cdot\cdot\cdot&n_{m}\end{bmatrix}$

The matrix resulting from the product of these two matrices is itself
diagonal, with its determinant expressed as:

$\det(\textbf{A}_{1}^{T}\textbf{A}_{1})=\prod_{i=1}^{m}n_{i}$ (7.6)

Thus, we can compute the term
$|\textbf{A}_{1}^{T}\textbf{A}_{1}|/|\textbf{A}_{0}^{T}\textbf{A}_{0}|$ in
equation (7.3) as:

$|\textbf{A}_{1}^{T}\textbf{A}_{1}|/|\textbf{A}_{0}^{T}\textbf{A}_{0}|=\frac{\prod_{i=1}^{n}n_{i}}{n}$
(7.7)

In this case, the minimal training sample requires at least one observation in
each group plus one extra observation in any of the groups to estimate
$\sigma^{2}$. Thus, we require that $n_{i}=1,i=1,...,j-1,j+1,...,n$, and
$n_{j}=2$, for some $j\in[1,m]$.

For the case mentioned above, we can use equation (7.4), and the result
obtained in (7.6) to get $c_{1}$ as:

$\det(\textbf{A}_{1}(\ell)^{T}\textbf{A}_{1}(\ell))=(1\cdot\cdot\cdot 1\cdot
2\cdot 1\cdot\cdot\cdot 1)=2$ (7.8)

Similarly, $c_{0}$ can be derived as:

$\det(\textbf{A}_{0}(\ell)^{T}\textbf{A}_{0}(\ell))=n=\sum_{j=1}^{m}n_{j}=(1+\cdot\cdot\cdot+1+2+1+\cdot\cdot\cdot+1)=m+1$
(7.9)

Thus, utilizing the outcomes from (7.8) and (7.9), we can evaluate the
expression in equation (7.4) as follows:

$\frac{c_{0}}{c_{1}}=\big{(}\frac{\det(\textbf{A}_{0}(\ell)^{T}\textbf{A}_{0}(\ell))}{\det(\textbf{A}_{1}(\ell)^{T}\textbf{A}_{1}(\ell))}\big{)}^{\frac{1}{2}}=\big{(}\frac{m+1}{2}\big{)}^{\frac{1}{2}}$
(7.10)

Now, by substituting (7.10) and (7.7) into (7.3), derive the following
expression for the Bayes factor:

$B^{SS}_{01}=\big{(}\frac{m+1}{2}\big{)}^{\frac{1}{2}}\big{[}\displaystyle\prod_{i=1}^{m}n_{i}\big{/}n\big{]}^{\frac{1}{2}}\big{[}1+\frac{(m-1)}{n-m}F\big{]}^{-n/2}$
(7.11)

The SS Bayes Factor can alternatively be expressed compactly as:

$\boxed{B^{SS}_{01}=\Bigg{(}\frac{\big{(}\frac{m+1}{2}\big{)}\big{[}\prod_{i=1}^{m}n_{i}\big{/}n\big{]}}{\big{[}1+\frac{(m-1)}{n-m}F\big{]}^{n}}\Bigg{)}^{\frac{1}{2}}}$
(7.12)

Furthermore, relating the $F$ statistics to the p-value can be achieved
through:

$F_{\nu_{1},\nu_{2}}=qf(1-p,\nu_{1},\nu_{2})$

Here, $qf$ denotes the quantile function for the $F$ distribution, while
$\nu_{1}$ and $\nu_{2}$ represent the degrees of freedom. For this case,
$\nu_{1}=p_{1}-p_{0}$ and $\nu_{2}=n-p_{1}$. Consequently, expressing the SS
Bayes factor in terms of the p-value yields:

$\boxed{B^{SS}_{01}(p)=\Bigg{(}\frac{\big{(}\frac{m+1}{2}\big{)}\big{[}\prod_{i=1}^{m}n_{i}\big{/}n\big{]}}{\big{[}1+\frac{(m-1)}{n-m}qf(1-p,p_{1}-p_{0},n-p_{1})\big{]}^{n}}\Bigg{)}^{\frac{1}{2}}}$

This formula offers a straightforward method to convert p-values into Bayes
factors. It’s akin to the Sellke et al. (2001) [13] bound, but distinguishes
itself by being dynamic, enhancing accuracy with increasing information. This
technique is applicable to any Bayes factor that relies on F-statistics.

Our subsequent aim involves computing $B^{GSS}_{01}$ as in equation
$\eqref{BFGSS}$ for ANOVA models and comparing it with their findings.
Considering that $A_{0}(\ell)$ and $A_{1}(\ell)$ remain independent of
observations in ANOVA models, and the square root is a monotonic function, we
can rewrite the expression as:

$\displaystyle\sup_{y(\ell)\in
D_{n}}\sqrt{\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{i}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q-p_{1}}}}=\sqrt{\frac{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}\displaystyle\sup_{y(\ell)\in
D_{n}}\frac{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{i}-p_{0}}}{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{i}-p_{1}}}}$

The General Theoretical and Empirical SS Bayes factors for ANOVA can be
expressed respectively as:

$\boxed{B_{10}^{GSSA}=\sqrt{\frac{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}\frac{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}\frac{[\textbf{R}_{0}(\textbf{y})]^{n+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}\sup_{y(\ell)\in
D}\frac{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}}}$
$\boxed{B_{10}^{GESSA}=\sqrt{\frac{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}\frac{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}\frac{[\textbf{R}_{0}(\textbf{y})]^{n+q_{0}-p_{0}}}{[\textbf{R}_{1}(\textbf{y})]^{n+q_{1}-p_{1}}}\displaystyle{\max_{\ell=1,...,L}}\frac{[\textbf{R}_{1}(\textbf{y}(\ell))]^{n_{01}+q_{1}-p_{1}}}{[\textbf{R}_{0}(\textbf{y}(\ell))]^{n_{01}+q_{0}-p_{0}}}}}$

####  7.3.1 The Full Jeffrey’s (no independence) $q_{k}=p_{k},$, for $k=i,j$
for ANOVA

In accordance with our earlier calculations pertaining to the General
Empirical SS Bayes factor for ANOVA, we apply the methodology proposed by
Spiegelhalter and Smith (1982) [15]. This involves the selection of the full
Jeffrey’s prior where $q_{i}=p_{i}$. Through the utilization of an imaginary
training sample, which renders maximum support for $M_{0}$, we ascertain that
$F{{}_{y}(\ell)}=0$. Consequently, we obtain the General Theoretical SS Bayes
Factor for ANOVA under the Full Jeffrey’s prior

$B_{01}^{GESSA_{FJP}}=\sqrt{\frac{m+1}{2}\big{[}\displaystyle\prod_{i=1}^{m}n_{i}\big{/}n\big{]}\Big{(}\frac{\textbf{R}_{1}(\textbf{y})}{\textbf{R}_{0}(\textbf{y})}\Big{)}^{n}\min_{\ell=1,...,L}\Big{(}\frac{\textbf{R}_{0}(\textbf{y}(\ell))}{\textbf{R}_{1}(\textbf{y}(\ell))}\Big{)}^{m+1}}.$

Now, expressing the ratio of the residual sum of squares as an F-statistics
results in:

$B_{01}^{GESSA_{FJP}}=\sqrt{\frac{m+1}{2}\big{[}\displaystyle\prod_{i=1}^{m}n_{i}\big{/}n\big{]}\Big{(}1+\frac{(m-1)}{(n-m)}F_{y}\Big{)}^{-n}\max_{\ell=1,...,L}\Big{(}1+\frac{(p_{1}-p_{0})}{(n_{01}-p_{1})}F_{y(\ell)}\Big{)}^{m+1}}$

The General Theoretical SS Bayes factor for ANOVA results in:

$B_{10}^{GSSA_{FJP}}=\sqrt{\frac{2}{m+1}\big{[}\displaystyle
n\big{/}\prod_{i=1}^{m}n_{i}\big{]}\Big{(}1+\frac{(m-1)}{(n-m)}F_{\textbf{y}}\Big{)}^{n}\sup_{y(\ell)\in
D}\Big{[}\frac{\textbf{R}_{1}(\textbf{y}(\ell))}{\textbf{R}_{0}(\textbf{y}(\ell))}\Big{]}^{m+1}}$

Observing that $\textbf{R}_{0}(y(\ell))\geq\textbf{R}_{1}(y(\ell))$, it
follows that:

$\sup_{y(\ell)\in
D}\frac{[\textbf{R}_{1}(\textbf{y}(\ell))]}{[\textbf{R}_{0}(\textbf{y}(\ell))]}=1\Longrightarrow\sup_{y(\ell)\in
D}\Big{[}\frac{\textbf{R}_{1}(\textbf{y}(\ell))}{\textbf{R}_{0}(\textbf{y}(\ell))}\Big{]}^{m+1}=1,\hskip
5.0pt,\forall m>0$

Consequently, the General Theoretical SS Bayes Factor, employing the full
Jeffrey’s prior for ANOVA models, simplifies to:

$\boxed{B_{10}^{GSSA_{FJP}}=\sqrt{\frac{2}{m+1}\big{[}\displaystyle
n\big{/}\prod_{i=1}^{m}n_{i}\big{]}\Big{(}1+\frac{(m-1)}{(n-m)}F_{\textbf{y}}\Big{)}^{n}}=B_{10}^{SS}.}$

Alternatively, this Bayes factor can be expressed as:

$\boxed{B_{10}^{GSSA_{FJP}}=\sqrt{\frac{2}{m+1}\big{[}\displaystyle
n\big{/}\prod_{i=1}^{m}n_{i}\big{]}\Big{[}\frac{\textbf{R}_{0}(\textbf{y})}{\textbf{R}_{1}(\textbf{y})}\Big{]}^{n}}.}$

####  7.3.2 The Reference prior ($q_{k}=0,k=i,j$) for ANOVA

In this section, we delve into the theoretical SS Bayes factor for the ANOVA
scenario, employing the reference prior obtained by setting $q_{k}=0$ for
$k=i,j$. The theoretical SS Bayes Factor for ANOVA, specifically under the
reference prior context, is formulated as:

$\boxed{B_{10}^{GSSA_{RP}}=\sqrt{\frac{|\textbf{A}_{1}^{T}(\ell)\textbf{A}_{1}(\ell)|}{|\textbf{A}_{0}^{T}(\ell)\textbf{A}_{0}(\ell)|}\frac{|\textbf{A}_{0}^{T}\textbf{A}_{0}|}{|\textbf{A}_{1}^{T}\textbf{A}_{1}|}\frac{[\textbf{R}_{0}(\textbf{y})]^{n+1}}{[\textbf{R}_{1}(\textbf{y})]^{n-(m+1)}}\sup_{y(\ell)\in
D}\Big{(}\frac{1}{\textbf{R}_{0}(\textbf{y}(\ell))}\Big{)}^{m}}}$

Upon inspection of the aforementioned equation, it becomes clear that under
this particular prior, the theoretical SS Bayes factor tends to either zero or
infinity. Consequently, the computation of the theoretical SS BF with this
prior is rendered unfeasible. This implies that for reference priors, the SS
bounds do not exist. This limitation is not previously highlighted in the
works of Smith and Spiegelhalter (1980)[14] and (1982) [15], thereby
underscoring an important consideration for future research in this domain.

####  7.3.3 The Modified Jeffrey’s Prior ($q_{j}=p_{j}-p_{i},q_{i}=0$) for
ANOVA

The General Empirical SS Bayes Factor for ANOVA, under the Modified Jeffrey’s
prior results in:

$\boxed{B_{10}^{GESSA_{MJP}}=\sqrt{\frac{2}{m+1}\big{[}n\big{/}\displaystyle\prod_{i=1}^{m}n_{i}\big{]}\Big{(}\frac{\textbf{R}_{0}(\textbf{y})}{\textbf{R}_{1}(\textbf{y})}\Big{)}^{n-1}\max_{\ell=1,...,L}\Big{(}\frac{\textbf{R}_{1}(\textbf{y}(\ell))}{\textbf{R}_{0}(\textbf{y}(\ell))}\Big{)}^{m}}.}$

Similarly, the general Theoretical SS Bayes Factor for ANOVA, under the
Modified Jeffrey’s prior results in:

$\boxed{B_{10}^{GSSA_{MJP}}=\sqrt{\frac{2}{m+1}\big{[}n\big{/}\displaystyle\prod_{i=1}^{m}n_{i}\big{]}\Big{(}\frac{\textbf{R}_{0}(\textbf{y})}{\textbf{R}_{1}(\textbf{y})}\Big{)}^{n-1}\sup_{y(\ell)\in
D}\Big{(}\frac{\textbf{R}_{1}(\textbf{y}(\ell))}{\textbf{R}_{0}(\textbf{y}(\ell))}\Big{)}^{m}}.}$

Alternatively, this Bayes Factor can be expressed in terms of the $F$
statistics as:

$\boxed{B_{10}^{GSSA_{MJP}}=\sqrt{\frac{2}{m+1}\big{[}n\big{/}\displaystyle\prod_{i=1}^{m}n_{i}\big{]}\Big{(}1+\frac{(m-1)}{(n-m)}F_{y}\Big{)}^{n-1}}.}$

###  7.4 Conclusions

In conclusion, the ANOVA case analysis has provided valuable insights
regarding the use of different priors in the context of Full Jeffreys,
Modified Jeffreys, and Reference Prior. The SS bounds obtained for Full
Jeffreys and Modified Jeffreys are informative, useful, and exhibit close
proximity to each other. Importantly, it has been demonstrated that the Bayes
Factor can be expressed as a function of the F-Statistics under both priors.
However, in the case of the Reference Prior, the SS bound fails to provide
informative results, highlighting its lack of usefulness. Notably, it is
intriguing to observe the sensitivity of the SS bound to the initial objective
prior, which further supports the conclusions drawn by Berger and Pericchi
(1996) [3] and others that the Modified Jeffreys prior is a superior choice,
particularly within the Linear Gaussian Model and potentially in a broader
context. Additionally, it may be of interest to explore the extent to which
the Modified Jeffreys prior possesses additional properties, such as matching,
besides its inherent simplicity and the one-on-one relationship it establishes
with the F-Statistics.

##  Chapter 8 Separated Hypothesis Testing

Traditional hypothesis testing typically involves comparing a single
hypothesis to a pre-determined null hypothesis, assuming a single underlying
distribution. However, in certain scenarios, there arises a need to examine
and compare two distinct distributions simultaneously. This is the essence of
separated hypothesis testing, a concept that expands upon conventional
approaches by accommodating situations where the hypothesis being tested
corresponds to two separate populations or data sources. By addressing this
unique requirement, separated hypothesis testing enables researchers to gain
valuable insights into the potential disparities or similarities between two
distinct distributions, thereby enhancing the accuracy and comprehensiveness
of statistical inference. In this chapter, we delve into the theoretical
underpinnings, methodology, and applications of separated hypothesis.

We explored the concept of separated hypothesis testing through practical
examples involving the comparison of different probability distributions.
Specifically, we consider the scenarios of Poisson vs. Geometric and Poisson
vs. Negative Binomial distributions. We computed IBF bounds for both
comparisons, ensuring the reliability and consistency of our results. To
validate our findings, we conducted simulations where we progressively
augmented the sample size, observing a numerical convergence toward accurate
decisions in our experiments. Each example is accompanied by a comprehensive
discussion of the methodologies and code used in our analyses. The figures
presented serve as compelling visual evidence substantiating the validity of
our conclusions.

###  8.1 The Poisson vs. Geometric Separated Hypothesis Test

Let’s consider comparing two distinct models: the Poisson and the Geometric
distributions, where $\textbf{y}=(y_{1},...,y_{n})$ represents independently
and identically distributed (i.i.d.) observations. This frames our model
selection problem as follows:

$M_{0}:Y\sim Poisson(\lambda)\mbox{ vs }M_{1}:Y\sim Geometric(\theta)$

The densities for each model are:

$f_{0}(y|\lambda)=\frac{e^{-\lambda}\lambda^{y}}{y!},\lambda>0,\hskip
3.0pt\mbox{ and }\hskip
3.0ptf_{1}(y|\theta)=\theta(1-\theta)^{y},\theta\in(0,1),y=0,1,2,...$

The Jeffreys’s priors and likelihoods associated with $M_{0}$ and $M_{1}$
yield the following outcomes:

$\pi_{0}^{J}(\lambda)=\frac{c_{0}}{\lambda^{1/2}},\hskip
10.0pt\pi_{1}^{J}(\theta)=\frac{c_{1}}{\theta(1-\theta)^{1/2}}$
$f_{0}(\textbf{y}|\lambda)=\frac{e^{-n\lambda}\lambda^{\sum_{i=1}^{n}y_{i}}}{\displaystyle\prod_{i=1}^{n}y_{i}!},\hskip
3.0pt\mbox{ and }\hskip
3.0ptf_{1}(\textbf{y}|\theta)=\theta^{n}(1-\theta)^{\sum_{i=1}^{n}y_{i}}$

Given that both models entail a single parameter, the minimal training sample
size required is one. Utilizing the priors and likelihoods, we proceed to
compute the marginals and subsequently derive the Bayes factors as follows:

$m_{0}(\textbf{y})=\frac{c_{0}}{\displaystyle\prod_{i=1}^{n}y_{i}!}\displaystyle\int_{0}^{\infty}e^{-n\lambda}\lambda^{\sum_{i=1}^{n}y_{i}}\frac{1}{\lambda^{1/2}}d\lambda$
$m_{1}(\textbf{y})=c_{1}\displaystyle\int_{0}^{1}\theta^{n}(1-\theta)^{\sum_{i=1}^{n}y_{i}}\frac{1}{\theta(1-\theta)^{1/2}}d\theta$

After straightforward algebraic manipulation, we find that the integral in
$m_{0}(\textbf{y})$ follows a Gamma distribution, while the integral in
$m_{1}(\textbf{x})$ corresponds to a Beta distribution, wherein:

$m_{0}(\textbf{y})=\frac{c_{0}\Gamma(\sum
y_{i}+\frac{1}{2})}{\displaystyle\prod_{i=1}^{n}y_{i}!}(1/n)^{(\sum
y_{i}+\frac{1}{2})}$
$m_{1}(\textbf{y})=c_{1}\displaystyle\int_{0}^{1}\theta^{n-1}(1-\theta)^{\sum_{i=1}^{n}y_{i}-1/2}d\theta=c_{1}\frac{\Gamma(n)\Gamma(\sum
y_{i}+1/2)}{\Gamma(n+\sum y_{i}+1/2)}$

Thus, the Bayes Factor results in:

$B_{01}^{N}(\textbf{y})=\frac{c_{0}}{c_{1}}\frac{\Gamma(n+\sum
y_{i}+1/2)}{\displaystyle\prod_{i=1}^{n}y_{i}!\Gamma(n)n^{(\sum
y_{i}+\frac{1}{2})}}.$

The Empirical IBF upper bound can be obtained by computing:

$\overline{B}^{I^{*}}_{10}(\textbf{y})=B_{10}^{N}(\textbf{y})\displaystyle\sup_{y(\ell)\in
D_{n,k}}\frac{m_{0}(y(\ell))}{m_{1}(y(\ell))}$

Evidently, the following marginals are required:

$m_{0}(y(\ell))=c_{0}\frac{\Gamma(y(\ell)+1/2)}{y(\ell)!},\hskip
15.0ptm_{1}(y(\ell))=c_{1}\frac{\Gamma(1)\Gamma(y(\ell)+1/2)}{\Gamma(1+y(\ell)+1/2)}$

By using those marginals, we can calculate the Bayes factor for $y(\ell)$ as:

$B_{01}^{N}(y(\ell))=\frac{c_{0}}{c_{1}}\frac{1}{\Gamma(y(\ell)+1+1/2)}$

The above Bayes factor, as a function of $y(\ell)$ is decreasing. Therefore,
calculating the theoretical supremum, which occurs at $y(\ell)=0$, results in:

$\sup_{y(\ell)\in
D}B_{01}^{N}(y(\ell))=\frac{c_{0}}{c_{1}}\frac{1}{\Gamma(1+1/2)}$

Utilizing the identity $\Gamma(n+1/2)=(2n!\sqrt{\pi})/(4^{n}n!)$, we derived:

$\sup_{y(\ell)\in
D}B_{01}^{N}(y(\ell))=\frac{c_{0}}{c_{1}}\frac{2}{\sqrt{\pi}}$

This results in the Theoretical IBF upper bound:

$\boxed{\overline{B}_{10}^{I}(\textbf{y})=\frac{\displaystyle\prod_{i=1}^{n}y_{i}!(n-1)!n^{(n\bar{y}+\frac{1}{2})}}{\Gamma(n(1+\bar{y})+1/2)}\frac{2}{\sqrt{\pi}}}$

And the Theoretical IBF lower bound:

$\boxed{\underline{B}_{01}^{I}(\textbf{y})=\frac{\Gamma(n(1+\bar{y})+1/2)}{\displaystyle\prod_{i=1}^{n}y_{i}!(n-1)!n^{(n\bar{y}+\frac{1}{2})}}\frac{\sqrt{\pi}}{2}}$

Similarly, the empirical IBF upper bound is:

$\boxed{\overline{B}_{10}^{I^{*}}(\textbf{y})=\frac{\displaystyle\prod_{i=1}^{n}y_{i}!(n-1)!n^{(n\bar{y}+\frac{1}{2})}}{\Gamma(n(1+\bar{y})+1/2)}\frac{1}{\Gamma(y_{(1)}+1+\frac{1}{2})},\hskip
15.0pty_{(1)}=\displaystyle\min_{i=1,...,n}\\{y_{i}\\}}$

And the empirical IBF lower bound is:

$\boxed{\underline{B}_{01}^{I^{*}}(\textbf{y})=\frac{\Gamma(n(1+\bar{y})+1/2)\Gamma(y_{(n)}+1+\frac{1}{2})}{\displaystyle\prod_{i=1}^{n}y_{i}!(n-1)!n^{(n\bar{y}+\frac{1}{2})}},\hskip
15.0pty_{(n)}=\displaystyle\max_{i=1,...,n}\\{y_{i}\\}}$

We aim to ascertain the reliability and consistency of our obtained results by
implementing an R script specifically designed to calculate the Bayes factor
in two distinct scenarios: when the null hypothesis ($H_{0}$) holds true and
when the alternative hypothesis ($H_{1}$) is true. By conducting this
analysis, we will demonstrate that when the samples are drawn from the Poisson
distribution, the computed Bayes factor approaches zero, indicating stronger
support for $H_{0}$. Conversely, when the samples originate from the Geometric
distribution, the Bayes factor exhibits values significantly greater than one,
favoring $H_{1}$. Furthermore, we will calculate the empirical IBF bounds,
providing numerical bounds for both $B_{01}^{N}(\textbf{y})$ and
$B_{10}^{N}(\textbf{y})$. The R script below generates one hundred simulations
of thirty independent samples from the Poisson distribution and from the
Geometric distribution. By employing this approach, we calculated the strength
of evidence in favor of the null hypothesis in both scenarios. As anticipated,
the obtained values are consistent with our prior knowledge about the data
that we simulated, indicating the expected outcome of growing indefinitely
when the null is false and going to zero under the null. We smoothed the
results by calculating the average value for each sample size over all of the
simulations.

R Program

⬇

1

2 #100 simulations of 30 samples each

3

4 #Null hypothesis true, lambda =1