# GRASP: A Goodness-of-Fit Test for Classification Learning

Adel Javanmard Data Sciences and Operations Department, University of Southern
CaliforniaA. Javanmard is partially supported by the Sloan Research Fellowship
in mathematics, an Adobe Data Science Faculty Research Award and the NSF
CAREER Awards DMS-1844481 and DMS-2311024. Mohammad Mehrabi11footnotemark: 1
The names of the authors are in alphabetical order.

###### Abstract

Performance of classifiers is often measured in terms of average accuracy on
test data. Despite being a standard measure, average accuracy fails in
characterizing the fit of the model to the underlying conditional law of
labels given the features vector ($Y|X$), e.g. due to model misspecification,
over fitting, and high-dimensionality. In this paper, we consider the
fundamental problem of assessing the goodness-of-fit for a general binary
classifier. Our framework does not make any parametric assumption on the
conditional law $Y|X$, and treats that as a black box oracle model which can
be accessed only through queries. We formulate the goodness-of-fit assessment
problem as a tolerance hypothesis testing of the form

$H_{0}:{\mathbb{E}}\Big{[}D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\Big{]}\leq\tau\,,$

where $D_{f}$ represents an $f$-divergence function, and $\eta(x)$,
$\widehat{\eta}(x)$ respectively denote the true and an estimate likelihood
for a feature vector $x$ admitting a positive label. We propose a novel test,
called GRASP for testing $H_{0}$, which works in finite sample settings, no
matter the features (distribution-free). We also propose model-X GRASP
designed for model-X settings where the joint distribution of the features
vector is known. Model-X GRASP uses this distributional information to achieve
better power. We evaluate the performance of our tests through extensive
numerical experiments.

##  1 Introduction

In classification learning, one is given a set of training data
$\\{(x_{i},y_{i})\\}_{i\leq n}$ (with $x_{i}$ representing multi-dimensional
feature vector and $y_{i}$ representing label variables), and aims to learn a
model which can be used to predict the labels on new feature vectors.
Classification algorithms are backbone of machine learning systems and
undoubtedly one of the prominent statistical learning tools in data
processing. There has been a plethora of classification methods proposed in
the literature ranging from logistic regression and generalized linear models
to more complex models such as boosting, random forests, and neural networks.
In practice, the performance of these methods is often assessed in terms of
accuracy on a test (hold-out) dataset, with the hope that it is a good
indicator of the predictive performance of the model on unseen data points.
Despite being used widely, the classification accuracy alone does not
necessarily characterize the deviations of the learnt model from the
underlying data generating process. Indeed relying solely on it as a measure
of performance can be misleading due to model misspecification, and over/under
fitting. This leads to the following fundamental question:

(*) How well does a classifier learn the ground truth data generating law
between the feature vector and the label?

Developing a statistical test for the above question has a myriad of
applications. It can flag the inherent and systematic flaws of a model, and
its poor generalization to unseen population. Further, it provides a more
holistic and honest assessment of the model performance, which is of paramount
importance with the rise of reproducibility issues in modern data analysis. In
particular, with the practice of data sharing, many datasets are used
routinely as benchmark to compare different models. However, perpetual use of
public datasets, without proper mechanism to access them to ensure validity of
inferences, causes spurious discoveries and overfitting; learning models which
performs well on benchmark datasets but generalize poorly to unseen datasets
[RAN14, JM18, DFH+15, DFH+17]. Another application of such test is for models
built by commercial “machine learning as a service” providers such as Google
and Amazon. They provide platforms where one can upload a dataset and a data
classification task and pay to construct a model. Therefore, it is important
to decide if the current model is sufficiently good (in a statistically sound
sense) or to continue training process, which incurs additional cost. Among
other applications, the proposed methodology can be used in $K$-fold framework
to choose the optimal model with respect to its goodness-of-fit, rather than
its empirical accuracy on the hold-out set.

In statistics, question (*) is often formulated as goodness-of-fit test.
However, most of the existing literature focus on specific parametric models,
such as logistic regression, and do not apply to more complicated models such
as neural network or random forest. An exception is the very recent seminal
work of [ZDY21]. We refer to Section 1.1 for a detailed discussion.

In this paper, we develop a goodness-of-fit test for a broad class of data
generating rules (unknown) and black-box models, with only query access. We
propose a novel methodology named GRASP (short for Goodness-of-fit with
Randomization and Scoring Procedure) for this end, which controls type I error
in finite sample settings, no matter the features, and does not make any
parametric assumption (distribution-free). We also propose a (less
conservative) variant test which comes with asymptotic validity. Both variants
show high power in identifying deviations of the classification procedure from
the true conditional law of the labels. We also consider model-X settings
where no knowledge of the conditional distribution labels ($Y|X$) is assumed,
but we do assume the joint distribution of the features $X$ is known, e.g., by
having access to abundant unlabeled data. This setting has been studied in
several recent work on variable selection and conditional testing; see
[CFJL18, BCS20, BCJW21, JM21] for a non-exhaustive list. We propose model-X
GRASP which leverages this information to obtain a better statistical power.

Our focus will be on the binary classification setup. We denote the feature
vectors by $x\in\mathcal{X}\subset{\mathbb{R}}^{d}$ and binary labels by
$y\in\\{0,1\\}$, and define the underlying conditional rule as
$\eta(x)={\mathbb{P}}(Y=+1|X=x)$. This simply reflects the likelihood of a
feature vector $x$ admitting the positive label. We have access to $\eta(x)$
only through a set of queries ${\cal D}$ consisting of $n$ samples drawn
i.i.d. from a common law $\mathcal{P}$ over $\mathcal{X}\times\\{0,1\\}$. We
are provided with an estimate model
$\widehat{\eta}:\mathcal{X}\rightarrow[0,1]$, and we would like to assess the
performance of this model. We assume that this model has been trained on a
dataset disjoint from ${\cal D}$ (${\cal D}$ has no share in the training
procedure.) The training mechanism is optional, e.g., can be a neural network
or decision tree, among many others. As a first step towards answering (*), we
impose the following hypothesis testing problem:

$\displaystyle H_{0}$
$\displaystyle:{\mathbb{E}}[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))]\leq\tau\,,$
(1)

where expectation is with respect to the distribution of $X$. The
$f$-divergence $D$ is a given metric to measure the distance between
$\widehat{\eta}(\cdot)$ and $\eta(\cdot)$ (See equation (2) for formal
definition). Different divergence functions $f$ lead to different distance
measures. For instance, setting $f(t)=\frac{1}{2}|t-1|$ gives us the total
variation distance ${\mathbb{E}}\left[|\widehat{\eta}(X)-\eta(X)|\right]$. By
letting $\tau$ be zero, one can test for the perfect alignment of the test
model $\widehat{\eta}$ and the ground truth rule $\eta$.

###  1.1 Related works

Goodness-of-fit. Once a model is fitted to data, it is important to assess the
quality of the fit. Several methods have been developed for testing goodness-
of-fit of generalized linear models under the low-dimensional setting ($d\ll
n$) with a focus on logistic regression and multinomial models [HL80, LCVH91,
Tsi80, SW91, LWY02, HH02, OR92, Far96]. For the high-dimensional setting,
[SB18] proposed a framework for testing goodness-of-fit of high-dimensional
linear models by using parametric bootstrap to calibrate the estimate model
for scaled residuals. [JSBS20] proposed the generalized residual prediction
(GRP) test for goodness-of-fit testing of high-dimensional generalized linear
models. The aforementioned works focus on a class of parametric models,
particularly for linear and generalized linear models. For a non-parametric
setting, [ZDY21] proposed the binary adaptive goodness-of-fit test (BAGofT),
where it tests for the perfect match of the test model and the ground-truth
conditional probability, in an asymptotic regime.

While the existing methods for testing goodness-of-fit mostly focus on
contexts where the estimate model has asymptotic convergence to the ground
truth law, we propose a method that is flexible for arbitrarily complex
classifiers independent of their predictive performance. In addition, it can
be used under the high-dimensional setting ($d\gg n$) as well and still
generates a high resolution p-value for moderately large number of samples. In
addition, we consider the tolerance testing scheme which is more general than
the perfect match testing (i.e. $\eta=\widehat{\eta}$, a.s.), making it more
useful in practice. Our proposed method also allows to consider a variety of
metrics, including the average absolute distance, and the excess cross
entropy, among many others.

We would like to highlight a distinction between our problem and the
conventional setup for testing goodness-of-fit. In our setup, the model
estimate $\widehat{\eta}$ is learned on the training data and is subsequently
evaluated on an independent test data. In hypothesis (3), the model estimate
$\widehat{\eta}$ is fixed and the randomness stems solely from the test data.
This differs from other goodness-of-fit setups, such as those described in
[ZDY21, JSBS20], where the statistical inference takes into account the
variability of the training set. For example, the phrase ‘finite sample size’
in those work refers to the size of the training data, while in our setting it
refers to the size of the test data.

Prediction error. Estimating the prediction performance of a model is one of
the core tasks in data-driven applications [HTFF09]. In particular, it can
serve as a benchmark for model selection. Covariance penalty (CP) and cross
validation (CV) are two of widely used methods to estimate the prediction
performance of a regression model [Aka98, Mal00, Efr86, Efr04]. However, It
has been shown recently that such methods are not statistically accurate in
measuring the out-of-the sample performance of models. In fact, [RT19, Wag20]
showed that covariance penalty and AIC statistic [Aka98] methods reflect the
in-sample prediction error. In addition, [BHT21] revealed potential
inaccuracies of cross validation, and argued that it indeed measures _average_
prediction accuracy over many _hypothetical_ datasets. In this work, we take
another perspective on model’s prediction error. For a broad class of models
(e.g., neural networks, decision trees, boosting algorithms, etc), and for the
widely used negative likelihood loss (a.k.a. cross entropy), our framework
allows to compare the model loss with the optimal oracle loss. Formally, for
the estimate probabilistic model $\widehat{p}_{Y|X}(y|x)$, for specific choice
of $f$-divergence function, we can test the null hypothesis testing:

$\displaystyle{\mathbb{E}}_{(x,y)\sim\mathcal{P}}[-\log\widehat{p}_{Y|X}(y|x)]-\min_{\widetilde{p}_{Y|X}}\left\\{{\mathbb{E}}_{(x,y)\sim\mathcal{P}}[-\log\widetilde{p}_{Y|X}(y|x)]\right\\}\leq\tau\,.$

Note that the minimum for the second term is achieved at
$\tilde{p}_{Y|X}=p_{Y|X}$.

Calibration. With the rise of deploying machine learning systems in real-life,
the confidence of these systems in their predictions is of a great importance.
Classification procedures often output a confidence value
$\widehat{\pi}\in[0,1]$ along with their predicted value $\widehat{y}$, which
is supposed to indicate the model certainty about $\widehat{y}$ being equal to
$y$ (correct prediction). Calibration of a model refers to aligning such
certainty with its long-run accuracy. A rather surprising observation made
recently is that many modern machine learning methods are not well calibrated
[GPSW17, nix19, VWA+19, WLZ19, KLM19]. Expected calibration error (ECE)
[NCH15] is a well-known metric to measure the calibration of models, where it
is formally defined as
${\mathbb{E}}_{\widehat{\pi}}[|{\mathbb{P}}(\widehat{y}=y|\widehat{\pi})-\widehat{\pi}|]$.

In [LHHD22], authors proposed a method for testing for the perfect calibration
of generic predictive models in multi-class classification problems. The
initial intuition behind model calibration is to test for the closeness of
classification models to the ground truth law, but mathematically has a
different formulation for the goodness-of-fit testing problem considered in
the current work. In an extreme instance, a model is well calibrated if
${\mathbb{E}}[\eta(X)|\widehat{\eta}(X)=\eta]=\eta$ (regression setting), but
in (1) with e.g, the average absolute difference metric, a perfect fit
($\tau$=0) corresponds to $\widehat{\eta}(X)=\eta(X)$, almost surely.

Hypothesis testing for nonparametric regression functions. In [MS18], the
authors construct confidence sets for the regression function in nonparametric
binary regression with an unknown design density. It is assumed that the the
regression function $\eta(.)$ and the marginal probability density function of
the features belong to a continuous class of Sobolev type spaces. Other than
results on adaptive parameter estimation, this work provides a framework for
testing the null hypothesis that the regression function is equal to 1/2
versus its alternative, while allowing the marginal density function to be a
general function in a Sobolev type space. Furthermore, it is shown that the
complexity of the null hypothesis does not affect the minimal rate of
separation between the null and the alternative. More on this line of
research, [LS99] considers the detection problem for a response function $f$
in a stochastic model ${\rm d}X(t)=f(t){\rm d}t+{\varepsilon}{\rm d}W(t)$ with
$W(t)$ being the standard Wiener process. The detection problem corresponds to
testing the null hypothesis $f=0$ against its alternative. In addition, [IS09]
considers simple hypothesis testing of the form $f=f_{0}$ in a multivariate
setup $y=f(x)+{\varepsilon}$ with isotropic Gaussian noise in $[0,1]^{d}$. We
refer to [IIS03] for more details on nonparametric hypothesis testing on
stochastic Gaussian models.

###  1.2 Summary of contributions and organization

In this paper, we introduce a novel method for tolerance testing of average
distance of classifiers to the underlying conditional law of labels. Our
proposed method, called GRASP can be used for arbitrarily complex black-box
models (with only query access), with no parametric assumption whatsoever, and
no matter the features. GRASP can be used in high-dimensional setting where
the features dimension can exceed the sample size. We propose two variants of
GRASP , one controls the type I error in finite sample settings, and the other
(less conservative) version controls the type I error in asymptotic setting.
Through extensive numerical experiments, we show that GRASP achieves high
statistical power.

For model-X settings, where the distribution of features can be well
approximated, we propose model-X GRASP which is built upon similar ideas as in
the distribution-free version but harnesses the knowledge of features
distribution to improve statistical power.

Here is an outline of the next sections:

  * •

Section 2: We first provide a brief review on $f$-divergence functions, and
then formulate the goodness-of-fit test as a tolerance hypothesis testing
problem. We end this section by a short review on some convex analysis
definitions.

  * •

Section 3: We start by focusing on the distribution-free regime, and propose
high-level intuitions behind the GRASP algorithm. We next move to formally
introduce GRASP procedure along with its test statistics and decision rules.
Next we prove that the type I error of GRASP (size of the test) can be
controlled in finite sample settings, for arbitrary data generative rules, and
general classifiers. We also introduce a less conservative version of GRASP
that has asymptotic control on the type I error. The GRASP test uses a score
function $T$ in forming the test statistics. It is worth noting that the size
of the test is controlled for arbitrary score functions. The choice of score
function though impacts the power of the test and we will discuss some choices
in Section 4.3. We conclude Section 3 by characterizing one-sided confidence
intervals and $p$-values for hypothesis (1).

  * •

Section 4: We move to the model-X setting, where abundant unlabeled data
points are available. We propose model-X GRASP that uses this data to learn
the features distribution and utilize it for a statistically more powerful
procedure. Similar to the distribution-free setting, we show that the size of
the test is controlled under the pre-determined level, for finite number of
data points. Further, a less conservative decision rule is introduced which
has asymptotic control on the type I error.

In Section 4.3 we discuss the role of score function $T$ on the power. We
derive the form of the optimal score function, which depends on the data
generative law $\eta(x)$. We discuss two approaches: $(i)$ _model-agnostic_
which replaces $\eta(x)$ by $1/2$ (random guessing) in the optimal score
function; $(ii)$ _GAN-based_ approach which uses generative adversarial
networks (GANs) to estimate the required densities to use in the optimal score
function.

  * •

Section 5: We provide extensive numerical experiments to evaluate the
performance of the distribution-free and model-X GRASP (type I control), the
power of these tests and the advantage of model-X framework in achieving a
higher power, as well as the impact of the score function $T$ on the power.

###  1.3 Notation

For an integer $k$ let $[k]$ stand for the set $\\{1,2,\dots,k\\}$. We denote
the distribution of a Bernoulli random variable with success parameter $p$
with $\mathsf{Bern}(p)$, and let $\mathsf{multi}(p_{1},\dots,p_{L})$ denote a
multinoulli distribution with $L$ categories, where the probability of
observing category $\ell\in[L]$ is $p_{\ell}$. For the case of
$p_{1}=\dots=p_{L}=1/L$ we use the shorthand $\mathsf{multi}(L)$. In addition,
for positive real values $a<b$, let $\mathsf{Unif}([a,b])$ indicate the
uniform distribution on the interval $[a,b]$. For a random variable $X$, we
write $\mathcal{L}(X)$ to refer to the probability density function of $X$. We
denote the density function of a chi-squared distribution with $L$ degrees of
freedom by $\chi^{2}_{L}$, and let $\chi^{2}_{L}(\beta)$ represents the
$\beta$-th quantile of an $\chi^{2}_{L}$ distribution. We use upper case
letters for random variables, and lower case letters to indicate deterministic
values, e.g, realizations of a random variable. We drop the subscript under
the expectation with showing the corresponding distribution, whenever it is
clear from the context.

##  2 Problem formulation

Under a binary classification setting, for a given test model
$\widehat{\eta}:\mathcal{X}\to[0,1]$, we are interested in measuring its
average distance to the oracle model $\eta(x)={\mathbb{P}}(y=+1|x)$. The test
model $\widehat{\eta}(\cdot)$ can be any arbitrarily complex predictive model,
such as a fitted logistic regression, or the last layer of a trained neural
network. We focus on a class of distances between $\eta$ and $\widehat{\eta}$
which are inspired by $f$-divergence of distributions. We start by the
definition of $f$-divergence of two density functions.

######  Definition 2.1.

($f$-divergence) Consider a convex and continuous function
$f:{\mathbb{R}}\to{\mathbb{R}}$. For two probability density functions $p,q$
that are defined with respect to the Lebesgue measure $\mu$ over
$\mathcal{X}\subset{\mathbb{R}}^{d}$, define the $f$-divergence between them
as the following

$D_{f}(p\|q)=\int qf\left(\frac{p}{q}\right){\rm d}\mu\,.$

Specializing this definition to Bernoulli distributions, with parameters
$a,b$, we obtain the following definition:

$D_{f}(\mathsf{Bern}(a)\|\mathsf{Bern}(b))=bf\left(\frac{a}{b}\right)+(1-b)f\left(\frac{1-a}{1-b}\right)\,.$
(2)

For the rest of this paper, we will focus on the class of measures
$D_{f}(\mathsf{Bern}(\eta(x))\|\mathsf{Bern}(\widehat{\eta}(x)))$ parametrized
by function $f$ as distance between the test model $\widehat{\eta}$ and the
ground truth model $\eta$. For some nonnegative value $\tau$, we consider the
following hypothesis testing problem:

$\displaystyle H_{0}$
$\displaystyle:{\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\leq\tau$
$\displaystyle H_{A}$
$\displaystyle:{\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]>\tau\,.$
(3)

with $H_{0}$ representing the null hypothesis and $H_{A}$ the alternative.

######  Lemma 2.2.

The followings hold:

  1. 1.

By choosing $f(t)=\frac{1}{2}|t-1|$ (total variation distance) we get

${\mathbb{E}}\left[D_{\mathsf{TV}}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]={\mathbb{E}}[|\widehat{\eta}(X)-\eta(X)|]\,.$

  2. 2.

For the choice of $f(t)=t\log t$ (KL divergence), we get

${\mathbb{E}}\left[D_{\mathsf{KL}}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]=\mathsf{CE}(\widehat{\eta})-\mathsf{CE}(\eta)\,,$

where for a model $\widehat{\eta}(x):\mathcal{X}\to[0,1]$, its cross entropy
loss is given by

$\mathsf{CE}(\widehat{\eta})=-{\mathbb{E}}[\eta(X)\log\widehat{\eta}(X)+(1-\eta(X))\log(1-\widehat{\eta}(X))]\,.$

The minimum cross-entropy loss is achieved when $\widehat{\eta}=\eta$.

  3. 3.

For the choice of $f(t)=(\sqrt{t}-1)^{2}$ (Hellinger distance), we get

$\displaystyle{\mathbb{E}}\left[D_{\mathsf{H}}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]={\mathbb{E}}\left[\Big{(}\sqrt{\eta(X)}-\sqrt{\widehat{\eta}(X)}\Big{)}^{2}+\Big{(}\sqrt{1-\eta(X)}-\sqrt{1-\widehat{\eta}(X)}\Big{)}^{2}\right]\,.$

The Proof of Lemma 2.2 is given in Section A.1. It is worth noting that by
considering different $f$, our framework allows to for a variety distance
measures between oracle model $\eta$ and the estimate model $\widehat{\eta}$.
We conclude this section by two definitions that will be used later in Section
4.

######  Definition 2.3.

(subdifferential) For a convex function $f:{\mathbb{R}}\to{\mathbb{R}}$, the
subdifferential $\partial f(t)$ at a point $t$ is given by the following set
of real values

$\partial f(t)=\\{u\in{\mathbb{R}}:f(s)-f(t)\geq u(s-t)\,,\forall
s\in{\mathbb{R}}\\}\,.$

In addition, for differentiable $f$, we have $\partial f(t)=f^{\prime}(t)$.

######  Definition 2.4.

(conjugate dual) The conjugate dual function of
$f:{\mathbb{R}}\to{\mathbb{R}}$ is defined as

$f^{*}(t)=\sup_{s\in{\mathbb{R}}}\;(st-f(s))\,.$

In addition, for convex lower semi-continuous $f$, we can write $f$ in terms
of its conjugate dual as

$f(s)=\sup_{t\in{\mathbb{R}}}\;(st-f^{*}(t))\,.$

##  3 Distribution-free setting

In this section, we develop a new methodology for testing the null hypothesis
in (3) without imposing any structure on the conditional law $\eta(x)$,
covariates $x$ distribution, or predictive model $\widehat{\eta}(x)$. In order
to provide a high level intuition behind the main idea of the procedure, we
first introduce a sampling scheme which characterizes the $f$-divergence
between the models $\eta$ and $\widehat{\eta}$ as the conditional
$f$-divergence distance of a sampled value and $\mathsf{Unif}[0,1]$
distribution.

######  Proposition 3.1.

For $(x,y)\sim\mathcal{P}$ let

$w=\begin{cases}u_{1}\sim\mathsf{Unif}[0,\widehat{\eta}(x)],&\quad y=1\,,\\\
u_{2}\sim\mathsf{Unif}[\widehat{\eta}(x),1],&\quad y=0\,.\\\ \end{cases}$ (4)

Then we have

${\mathbb{E}}_{X}\left[D_{f}(\mathcal{L}_{W|X}\|\mathsf{Unif}([0,1]))\right]={\mathbb{E}}_{X}[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))]\,.$

As a special case, if $\widehat{\eta}=\eta$ then
$W|X\sim\mathsf{Unif}([0,1])$.

The proof of Proposition 3.1 is given in Section A.2. Proposition 3.1 implies
that deviation of sampled $w$ from the uniform $[0,1]$ can be counted as
evidence for large distance of $\eta$ and $\widehat{\eta}$. In the next
section we elaborate the formal process to construct a set of statistics that
will be used later to define the decision rules for testing null hypothesis
(3).

###  3.1 Test statistic

Algorithm 1 describes the procedure for constructing the test statistic
$\mathbf{V}_{n,L}$. The construction consists in two main steps:
Counterfeit sampling. For each data point $(x,y)$, we follow the procedure
described in (4) and construct the sampled value $w$. We then construct
randomizations $\widetilde{w}_{1},\dots\widetilde{w}_{M}$ from the uniform
distribution $[0,1]$.
Score and label. In this step, by using a score function
$T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$, we try to score the original sample
$T(x,w)$ and corresponding values $T(x,\widetilde{w}_{j})$. Then the original
data point $(x,y)$ will be labeled based on the relative location of $T(x,w)$
among the counterfeit values $\\{T(x,\widetilde{w}_{j})\\}_{j=1:M}$. The final
output statistic $\mathbf{V}_{n,L}\in{\mathbb{R}}^{L}$ denote the count of
each label among the whole $n$ samples. Note that the number of labels $L$ is
an input parameter in Algorithm 1. The construction of test statistic
$\mathbf{V}_{n,L}$ is inspired by the PCR test proposed by [JM21], where a
similar structure (counterfeit sampling-scoring-labeling) is used but for a
different task, namely to test conditional independence between variables.
Also in [JM21], the counterfeits are drawn from a conditional distribution
(depending on variables under test), while here the counterfeits are drawn
from $\mathsf{Unif}([0,1])$, no matter the features or labels.

0: $n$ data points $(x_{j},y_{j})\in\mathcal{X}\times\\{0,+1\\}$, the model
$\widehat{\eta}:\mathcal{X}\to[0,1]$, a score function
$T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$, and integers $M,L\geq 1$ such that
$M+1=KL$ for some integer $K$.

0: Test statistics $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]\in{\mathbb{R}}^{L}$.

for  _$j\in\\{1,2,...,n\\}$_ do

  * •

Let

$w_{j}=\begin{cases}u_{1}\sim\mathsf{unif}[0,\widehat{\eta}(x_{j})],&\quad
y_{j}=1\,,\\\ u_{2}\sim\mathsf{unif}[\widehat{\eta}(x_{j}),1],&\quad
y_{j}=0\,.\\\ \end{cases}$

  * •

Draw $M$ i.i.d. samples $\widetilde{w}_{j}^{(1)},...,\widetilde{w}_{j}^{(M)}$
from $\mathsf{Unif}[0,1]$.

  * •

Use score function $T$ to score the initial sample $(x_{j},w_{j})$ and its $M$
constructed counterfeits
$\\{(x_{j},\widetilde{w}_{j}^{(1)}),...,(x_{j},\widetilde{w}_{j}^{(M)})\\}$:

$\displaystyle T_{j}$ $\displaystyle=T(x_{j},w_{j})$
$\displaystyle\widetilde{T}_{j}^{(i)}$
$\displaystyle=T(x_{j},\widetilde{w}_{j}^{(i)})\,,\quad\text{for }i\in[M]\,.$

  * •

Let $R_{j}$ denote the rank of $T_{j}$ among
$\\{T_{j},\widetilde{T}_{j}^{(1)},...,\widetilde{T}_{j}^{(M)}\\}$:

$R_{j}=1+\sum\limits_{i=1}^{M}{\mathbb{I}}_{\left\\{T_{j}\geq\widetilde{T}_{j}^{(i)}\right\\}}$

  * •

Assign label $L_{j}\in\\{1,2,...,L\\}$ to sample $j$ if $(L_{j}-1)K+1\leq
R_{j}\leq KL_{j}$.

end for

for  _$\ell\in\\{1,2,...,L\\}$_ do

  * •

Let $V_{\ell}$ be the number of samples with label $\ell$, i.e.
$V_{\ell}=\Big{|}\big{\\{}j\in\\{1,2,...,n\\}:L_{j}=\ell\big{\\}}\Big{|}\,.$

end for

Algorithm 1 Construction of GRASP test statistic (Distribution-free setting)

###  3.2 Decision rule

We introduce two decision rules based on the statistics
$\mathbf{V}_{n,L}=[V_{1},\dots V_{L}]$ given by Algorithm 1. For this end, we
first construct the following two test statistics that will be used later for
the decision rules:

$\displaystyle\begin{split}U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})=\min_{p\in{\mathbb{R}}^{L}}\quad&\frac{1}{n}\sum\limits_{\ell=1}^{L}\frac{\left(V_{\ell}-np_{\ell}\right)^{2}}{p_{\ell}+\frac{1}{L}}\\\
\text{s.t}\quad&p_{\ell}\geq
0\,,\quad\forall\ell\in[L]\,,\quad\sum\limits_{\ell=1}^{L}p_{\ell}=1\,,\quad\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})\leq\tau\,.\\\
\\\
{U}_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})=\min_{p\in{\mathbb{R}}^{L}}\quad&\frac{1}{n}\sum\limits_{\ell=1}^{L}\frac{\left(V_{\ell}-np_{\ell}\right)^{2}}{p_{\ell}}\\\
\text{s.t}\quad&p_{\ell}\geq
0\,,\quad\forall\ell\in[L]\,,\quad\sum\limits_{\ell=1}^{L}p_{\ell}=1\,,\quad\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})\leq\tau\,.\end{split}$
(5)

Note that the choice of function $f$ and the value of $\tau$ in the test
statistics ${U}_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$ and
${U}_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ are determined by the null
hypothesis (3). We next consider the following two decision rules based on
(5). At the significance level $\alpha$, the decision rules are given by

$\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}=\begin{cases}1~{}~{}(\text{reject
}H_{0})\,,&\text{if
}{U}_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})\geq\chi^{2}_{L-1}(1-\alpha)\,,\\\
0~{}~{}(\text{fail to reject }H_{0})\,,&\text{otherwise }\,.\\\ \end{cases}$
(6)

The next decision rule is given by

$\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}=\begin{cases}1~{}~{}(\text{reject
}H_{0})\,,&\text{ if }{U}_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\geq
L+\sqrt{\frac{2L}{\alpha}}\,,\\\ 0~{}~{}(\text{fail to reject
}H_{0})\,,&\text{otherwise }\,.\\\ \end{cases}$ (7)

The rationale behind the $\mathsf{asym}$ and $\mathsf{finite}$ names comes
from the fact that we show that later they have asymptotic and finite-sample
control guarantees for type I error, respectively.

As we show in the next lemma, the asymptotic rule is less conservative than
the finite rule.

######  Proposition 3.2.

For $L$ sufficiently large, the asymptotic decision rule (6) rejects more than
the finite decision rule (7). More precisely, for $L\geq 55$ and every
$\alpha\in(0,1)$, if $\Phi_{n,L,\alpha,\tau}^{\mathsf{finite}}=1$ then
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}=1$.

We refer to Section A.3 for the proof of Proposition 3.2.

We conclude this section by providing some insight behind the test statistics.
Consider the case of perfect alignment ($\tau=0$), for which we showed $w_{j}$
are uniform in $[0,1]$. Therefore, under the null each sample is identically
distributed as its counterfeits and its label, assigned in the last step of
Algorithm 1, follows $\mathsf{multi}(L)$ distribution. The test statistics in
(5) reduces to the Pearson’s chi-square test statistic. When $\tau>0$, we will
show in Theorem 3.5 that the labels follow a multinoulli distribution
$\mathsf{multi}(p_{1},\dotsc,p_{L})$ which should be close to the uniform
multinoulli distribution $\mathsf{multi}(L)$ in the sense that
$\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})\leq\tau$. However, the nominal
probabilities $p_{\ell}$ are unknown and so in our construction of test
statistics (5), we form an optimization problem over $p_{\ell}$ to impose this
constraint and then consider the minimum Pearson’s chi-square test value over
the feasible probability vectors $(p_{1},\dotsc,p_{L})$.

###  3.3 Size of the test

We will show that the proposed decision rules control the size of the test
under the pre-assigned level $\alpha$. We first state the following technical
assumption which posits a density function for specific conditional law for
$T(X,W)$ and $T(X,\widetilde{W})$.

######  Assumption 3.3.

Consider the following conditional cumulative distribution functions:

$\displaystyle F_{o}(t;x)$ $\displaystyle={\mathbb{P}}(T(X,W)\leq t|X=x)\,,$
$\displaystyle F_{c}(t;x)$ $\displaystyle={\mathbb{P}}(T(X,\widetilde{W})\leq
t|X=x)\,.$

where $(X,W)$ is is given by (4), and $\widetilde{W}$ is drawn independently
from $\mathsf{Unif}[0,1]$. Assume that the density functions of these cdfs
exist and indicate them by $f_{o}(t;x)$ and $f_{c}(t;x)$. In addition, a new
measure on the interval $[0,1]$ can be considered with cdf
$\omega(u;x):=F_{o}(F_{c}^{-1}(u;x);x)$. Suppose that this measure is
absolutely continuous with respect to the Lebesgue measure, and denote its
density function (Radon–Nikodym derivative) by $\psi(u;x)$.

######  Proposition 3.4 (distribution-free).

Let $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]$ be the output of Algorithm 1, then
$\mathbf{V}_{n,L}$ has a multinomial distribution with $L$ outcomes and
nominal probabilities $p_{1},\dots,p_{L}$. In addition, under Assumption 3.3
the following holds:

$\displaystyle\frac{1}{L}\sum\limits_{\ell=1}^{L}{f(Lp_{\ell})}$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\big{)}\right]$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}_{W|X}\|\mathcal{L}_{\widetilde{W}|X}\big{)}\right]$
$\displaystyle={\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}_{W|X}\|\mathsf{Unif}([0,1])\big{)}\right]\,.$

The proof of this Proposition is given in Section A.4. It is easy to observe
that the second inequality in Proposition 3.4 can be achieved by using the
score function $T(x,w)=w$. The next theorem is an immediate consequence of
Proposition 3.1 and Proposition 3.4.

######  Theorem 3.5 (distribution-free).

Let $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]$ be the output of Algorithm 1 that has
a multinomial distribution with nominal probabilities $p_{1},\dots p_{L}$, as
per Proposition 3.4. Then for a score function
$T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$ satisfying Assumption 3.3, the
following upper bound holds:

$\displaystyle\frac{1}{L}\sum_{\ell=1}^{L}{f(Lp_{\ell})}\leq{\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,.$

Note that by definition of $f$-divergence for discrete distributions we have

$D_{f}(\mathsf{multi}(p_{1},\dots,p_{L})||\mathsf{multi}(L))=\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})\,.$

Therefore, Theorem 3.5 implies that under the null hypothesis (3) the
$f$-divergence between the multinoulli distribution
$\mathsf{multi}(p_{1},\dots,p_{L})$ and the uniform multinoulli distribution
$\mathsf{multi}(L)$ should be bounded by $\tau$. As can be observed from (5),
our decision rules are based on optimization problems over probability vectors
$\\{p_{\ell}\\}$, which minimizes a chi-squared type statistic subject to the
constraint $\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})\leq\tau$.

In other words, by virtue of Theorem 3.5, we can reduce the null hypothesis
regarding the quantity of interest
${\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]$
to a hypothesis on the quantity $\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})$
which is more amenable to test.

Algorithm 2 outlines a simple version of Algorithm 1 with the score function
$T(x,w)=w$. In this case, scores of counterfeits will be uniform random
variables on $[0,1]$, and we label the original scores simply by partitioning
the $[0,1]$ interval into $L$ subintervals of equal size. It is worth noting
that the second inequality of Proposition 3.4 becomes an equality for the
choice of score function $T(x,w)=w$. At first glance, this is expected to
always results in a higher power compared to other choices of $T$. This
argument is not valid though, since the gap in the first inequality of 3.5 is
undetermined, and so we keep both versions (score function $T(x,w)=w$, and
general score function) in the paper.

The next result provides control over the size of our test with the two
decision rules (7) and (6), in the finite sample and asymptotic settings.

0: $n$ data points $(x_{j},y_{j})\in\mathcal{X}\times\\{0,+1\\}$, the model
$\widehat{\eta}:\mathcal{X}\to[0,+1]$, and an integer $L\geq 1$.

0: Test statistic $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]\in{\mathbb{R}}^{L}$.

for  _$j\in\\{1,2,...,n\\}$_ do

  * •

Let

$w_{j}=\begin{cases}u_{1}\sim\mathsf{unif}[0,\widehat{\eta}(x_{j})],&\quad
y_{j}=1\,,\\\ u_{2}\sim\mathsf{unif}[\widehat{\eta}(x_{j}),1],&\quad
y_{j}=0\,.\\\ \end{cases}$

  * •

Assign label $L_{j}\in\\{1,2,...,L\\}$ to sample $j$ if $\frac{L_{j}-1}{L}\leq
w_{j}\leq\frac{L_{j}}{L}$.

end for

for  _$\ell\in\\{1,2,...,L\\}$_ do

  * •

Let $V_{\ell}$ be the number of samples with label $\ell$, i.e.
$V_{\ell}=\Big{|}\big{\\{}j\in\\{1,2,...,n\\}:L_{j}=\ell\big{\\}}\Big{|}\,.$

end for

Algorithm 2 GRASP test statistic (distribution-free setting and with the score
function $T(x,w)=w$)

######  Theorem 3.6.

Let $\mathbf{V}_{n,L}=[V_{1},\dots,V_{L}]$ be the output of either Algorithms
1 or 2. Consider the decision rules $\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}$
and $\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}$, which are given respectively
in (6) and (7), with the test statistics
$U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$ and
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ as per (5). Under the null
hypothesis (3), we have:

${\mathbb{P}}\left(\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}=1\right)\leq\alpha\,,~{}~{}\lim\sup_{n\to\infty}{\mathbb{P}}\left(\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}=1\right)\leq\alpha\,.$

We refer to Section A.5 for the proof of Theorem 3.6. This implies that
deploying the decision rule $\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}$
controls the type I error at level $\alpha$, for every finite $n$. In
addition, the decision rule $\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}$ has an
asymptotic control over the type I error at the significance level $\alpha$.
Let us stress that the result of Theorem 3.6 is valid for every choice of
parameters $K,L$ and score function $T$.

###  3.4 P-values and one-sided confidence intervals

Considering the definition of rules $\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}$
and $\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}$ and the results of Theorem 3.6,
we construct the following p-values for the hypothesis testing problem in (3).

$\displaystyle
p_{n,L,\tau}^{\mathsf{finite}}=\begin{cases}1\,,&U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\leq
L\,,\\\
1\wedge\frac{2L}{(U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})-L)^{2}}\,,&\text{
otherwise}\,.\end{cases}$

$p_{n,L,\tau}^{\mathsf{asym}}=1-F_{L-1}(U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L}))\,,$

where $F_{L-1}(t)$ denote the cdf of a Chi-squared random variable with $L-1$
degrees of freedom. Super-uniformity of these $p$-values, under the null
hypothesis, follows simply from Theorem 3.6. Formally, for every $t\in[0,1]$
we have

$\displaystyle{\mathbb{P}}\left(p_{n,L,\tau}^{\mathsf{finite}}\leq
t\right)\leq t\,,\quad\forall n,L\geq 1\,.$
$\displaystyle\lim_{n\to\infty}\sup{\mathbb{P}}\left(p_{n,L,\tau}^{\mathsf{asym}}\leq
t\right)\leq t\,,\quad\forall L\geq 1\,.$

We next use the duality between confidence intervals and hypothesis testing to
construct a one-sided confidence interval for the quantity of interest
${\mathbb{E}}_{x}[D_{f}(\mathsf{Bern}(\eta(x))\|\mathsf{Bern}(\widehat{\eta}(x)))]$.
The intuition behind this construction comes from the fact that for a fixed
value $\mathbf{V}_{n,L}$, the test statistics
$U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$ and
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ are nonincreasing in $\tau$.

######  Proposition 3.7.

For $\alpha\in(0,1)$, let

$\displaystyle\tau^{\mathsf{finite}}_{n,L,\alpha}$
$\displaystyle=\sup\left\\{\tau\geq
0:U^{\mathsf{finite}}_{\tau}(\mathbf{V}_{n,L})\geq
L+\sqrt{2L/\alpha}\right\\}\,,$
$\displaystyle\tau^{\mathsf{asym}}_{n,L,\alpha}$
$\displaystyle=\sup\left\\{\tau\geq
0:U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})\geq\chi^{2}_{L-1}(1-\alpha)\right\\}\,.$

Then the followings hold

$\displaystyle{\mathbb{P}}\Big{(}{\mathbb{E}}\Big{[}D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\Big{]}\geq\tau^{\mathsf{finite}}_{n,L,\alpha}\Big{)}$
$\displaystyle\geq 1-\alpha\,,$
$\displaystyle\lim_{n\to\infty}\inf{\mathbb{P}}\Big{(}{\mathbb{E}}\Big{[}D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\Big{]}\geq\tau^{\mathsf{asym}}_{n,L,\alpha}\Big{)}$
$\displaystyle\geq 1-\alpha\,.$

Note that the probabilities in the above equations are with respect to the
randomness in $\tau^{\mathsf{finite}}_{n,L,\alpha}$ and
$\tau^{\mathsf{asym}}_{n,L,\alpha}$, which stem from the randomness in
$\mathbf{V}_{n,L}$.

We refer to Section A.7 for the Proof of Proposition 3.7.

###  3.5 Choice of the score function

In this section, we provide some insights on the choice of score function $T$
in Algorithm 1. First, from Theorem 3.5 we have

$\frac{1}{L}\sum_{\ell=1}^{L}{f(Lp_{\ell})}\leq{\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,.$

Note that in Algorithm 1 we test for
$\frac{1}{L}\sum_{\ell=1}^{L}{f(Lp_{\ell})}\leq\tau$ (which holds under null),
and therefore the statistical power of our method in the first place depends
on the gap between the quantity of primary interest
${\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]$
and $\frac{1}{L}\sum_{\ell=1}^{L}{f(Lp_{\ell})}$.

We next explore scenarios under which the gap in the chain of inequalities
presented in Proposition 3.4 can be tightened. Recall the chain of
inequalities summarized below:

$\displaystyle\frac{1}{L}\sum\limits_{\ell=1}^{L}{f(Lp_{\ell})}$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\big{)}\right]$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}_{W|X}\|\mathcal{L}_{\widetilde{W}|X}\big{)}\right]$
$\displaystyle={\mathbb{E}}\left[D_{f}\big{(}\mathcal{L}_{W|X}\|\mathsf{Unif}([0,1])\big{)}\right]$
$\displaystyle={\mathbb{E}}_{X}[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))]\,.\,,$
(8)

where the last step is proved in Proposition 3.1. Our goal is to develop a
score function that narrows the gap between the left and right-hand sides in
the above chain of inequalities. Note that the second inequality becomes tight
when $T(x,w)=w$, a simple score function outlined in Algorithm 2. While it is
a straightforward choice, it does not account for possible dependence between
$x$ and $w$, and it is also not clear how this choice would impact the first
inequality.

To examine the first inequality, we consider the limit case of $K,L\to\infty$
and see if the gap in the first inequality becomes tight in this asymptotic
case. Our next result answers this question in the negative. In contrast, we
later answer the similar question in the positive for the model-X setup.
Further details can be found in Proposition 4.3.

######  Proposition 3.8.

Recall the density function $\psi(u,X)$ from Assumption 3.3. Consider the
similar setup of Proposition 3.4 with the following two additional
assumptions: i) For a positive value $C$ we have $\psi(u;X)\leq C$ a.s. over
[0,1], and ii) $\psi(u,X)$ is differentiable on $(0,1)$ and there exists a
positive constant $B$ such that a.s. we have $|\frac{\partial}{\partial
u}\psi(u,X)|\leq B$. Then the following holds as $K,L$ grow to infinity:

$\displaystyle\lim_{L\to\infty}\lim_{K\to\infty}\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})$
$\displaystyle=\int_{0}^{1}f({\mathbb{E}}_{X}[\psi(u;X)]){\rm d}u$
$\displaystyle\leq{\mathbb{E}}_{X}\left[\int_{0}^{1}f(\psi(u;X)){\rm
d}u\right]$
$\displaystyle={\mathbb{E}}_{X}\left[D_{f}\big{(}\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\big{)}\right]\,.$
(9)

We refer to Section A.8 for the proof of Proposition 3.8.

The above Proposition highlights an important issue: even with $K,L\to\infty$,
we can still have a gap in the first inequality in (8) for any strictly convex
divergence function $f$, no matter the choice of score function $T$. One can
try to minimize this gap for some $T$, but apart from being a challenging
question it is not clear how this choice would impact the second inequality in
(8).

In Section 5, Experiment 5.2, we evaluate the performance of our test using
$T(x,w)=w$. We repeat the same numerical study in Experiment 5.3 using another
score function which also depends on $x$. Concretely, we regress $w$ on to
$x$, and take $T(x,w)=|w-x^{\sf T}{\widehat{\theta}}|$ the residuals. As we
see the former choice of $T(x,w)$ slightly outperforms the regression-type one
in power, for different choices of $f$-divergence. These experiments indicate
that the simple choice of $T(x,w)=w$ can be competitive and in general other
more complicated score score which takes into account the dependence between
$w$ and $x$ may achieve a lower power.

###  3.6 Testing perfect fit of the model

By setting $\tau=0$ in hypothesis (1), our framework allows to test whether we
have perfect fit for the model at hand $\widehat{\eta}$. We next provide a
modification of our test tailored for this special case, which has
demonstrated higher power in our numerical experiments. The main adjustment to
the framework is that the score functions are now defined at _dataset level_ ,
rather than _sample level_. Let $\mathbf{x},\mathbf{w}$ respectively stand for
$\\{x_{i}\\}_{i=1:n}$, $\\{w_{i}\\}_{i=1:n}$ with sampling process for $w_{i}$
described in (4). In addition, for $j\in[M]$ we define $\mathbf{w}_{j}$ to be
a vector of size $n$ with iid samples drawn from $\mathsf{Unif}([0,1])$. By
Proposition 3.1, under the null, we have $W|X\sim\mathsf{Unif}([0,1])$ and
therefore,
$T(\mathbf{x},\mathbf{w}),T(\mathbf{x},\mathbf{w}_{1}),\dotsc,T(\mathbf{x},\mathbf{w}_{M})$
are exchangeable. Using this observation, we construct the following $p$-value
for the null:

$p=\frac{1+\sum\limits_{j=1}^{M}{\mathbb{I}}\left(T(\mathbf{x},\mathbf{w})\geq
T(\mathbf{x},\mathbf{w}_{j})\right)}{M+1}\,,$ (10)

Using the exchangeability property, in the next proposition we show that the
above $p$-value is super uniform under null and hence we have control on the
size of the test.

The above construction of $p$-value is inspired by the Conditional
Randomization Test (CRT) [CFJL18]. Although the choice of score function $T$
is optional, it is recommended to use a score function that captures the
variation between $w$ and $x$, and hence $T(\mathbf{x},\mathbf{w})$ is smaller
than most of its counterfeits $T(\mathbf{x},\mathbf{w}_{j})$, resulting in a
small $p$-value under the alternative hypothesis. One simple choice is to use
the residual (e.g., mean squared error) of a linear regression model when
regressing variable $w$ on $x$. Alternatively, more complex predictive methods
such as Lasso, random forest, or neural networks can be used to fit a model,
and the residual of the fitted model can be reported as the score value.

Our next proposition shows the super-uniformity of p-value (10) under the
null, and its proof is deferred to Section A.9.

######  Proposition 3.9.

Under the null hypothesis (3) with $\tau=0$, the p-value (10) is super
uniform, i.e. for every $t\in[0,1]$ we have ${\mathbb{P}}(p\leq t)\leq t$.

In Example (5.5), we evaluate the performance of our test in a setting where
the size of the test data is significantly smaller than the training size, and
is comparable to the feature dimension. The modification made in defining the
test score at the data-set level allows us to achieve non-trivial power even
for small size test data. For the score function, we regress $w$ on
$\mathbf{x}$ using a three-layer neural network and define the score function
as the mean-squared-error of this model across the dataset. As we discuss,
although the test error is not an indicative measure of goodness-of-fit, our
test returns a significantly smaller $p$-value for the model that is closer to
the underlying conditional law. We refer to Example (5.5) for further details.

##  4 Model-X setting

In many applications, we may have access to a large amount of unsupervised
data (covariate data without corresponding labels) in addition to a limited
number of labeled observations. Motivated by this trend, [CFJL18] proposed
model-X setup where it assumes no knowledge of the conditional distribution of
labels $Y|X$, but assumes that the joint distribution of the covariates is
known, or can be well approximated. Model-X setup is also relevant in
experimental design, where we control the covariate distribution and so it is
known to us, e.g., in randomized controlled trials, as well as sensitivity
analysis of quantitive models [SRA+08]. This modeling assumption has paved the
way to address several statistical problems which are open or notoriously hard
otherwise. For example, [CFJL18] used this setup to extend the knockoff
framework of [BC15] to high-dimensional regime and provided a methodology for
variable selection with controlled false discovery rate in this regime. Also,
it develops a conditional randomization test using the knowledge of covariates
distribution (see [TVZ+22, BWBS20, JM21] for other related work on this
topic.)

Using similar ideas as in GRASP , in this section we propose model-X GRASP ,
which uses the covariates distribution to improve the power in flagging the
differences between the estimate model $\widehat{\eta}$ and the oracle model
$\eta$.

Our next proposition is analogous to Proposition 3.1 and extends it to the
model-X setup.

######  Proposition 4.1.

For $(x,y)\sim\mathcal{P}$ let

$w=\begin{cases}u_{1}\sim\mathsf{Unif}[0,\widehat{\eta}(x)],&\quad y=1\,,\\\
u_{2}\sim\mathsf{Unif}[\widehat{\eta}(x),1],&\quad y=0\,.\\\ \end{cases}$ (11)

Then the following holds:

$D_{f}\left(\mathcal{L}_{X,W}\|\mathcal{P}_{X}\times\mathsf{Unif}([0,1])\right)={\mathbb{E}}\left[D_{f}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,.$

The proof of Proposition 4.1 is given in Section A.6. Proposition 4.1 implies
that the distance of interest
${\mathbb{E}}[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))]$
can be seen as the $f$-divergence between multivariate random variables
$(X,W)$ and $({\widetilde{X}},\widetilde{W})$ with
$X,{\widetilde{X}}\sim\mathcal{P}_{X}$, independently and
$\widetilde{W}\sim\mathsf{Unif}([0,1])$.

###  4.1 Test statistic and decision rule

Algorithm 3 describes the procedure for constructing statistics
$\mathbf{V}_{n,L}$ under model-X setup for testing the null hypothesis in (3).
The overall procedure is similar to the distribution-free counterpart proposed
in Algorithm 1. The major difference is that here the counterfeits are of the
form $(\tilde{x},\widetilde{w})$, with the covariate component $\tilde{x}$
drawn from $\mathcal{P}_{X}$ and $\widetilde{w}$ drawn from
$\mathsf{Unif}[0,1]$. In contrast, in the distribution-free version (cf.
Algorithm 1) the covariate $x$ was fixed for an original sample and its
counterfeits.

0: $n$ data points $(x_{j},y_{j})\in\mathcal{X}\times\\{0,+1\\}$, the model
$\widehat{\eta}:\mathcal{X}\to[0,+1]$, a score function
$T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$, and integers $M,L\geq 1$ such that
$M+1=KL$ for some integer $K$.

0: Test statistics $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]\in{\mathbb{R}}^{L}$.

for  _$j\in\\{1,2,...,n\\}$_ do

  * •

Let

$w_{j}=\begin{cases}u_{1}\sim\mathsf{unif}[0,\widehat{\eta}(x_{j})],&\quad
y_{j}=1\,,\\\ u_{2}\sim\mathsf{unif}[\widehat{\eta}(x_{j}),1],&\quad
y_{j}=0\,.\\\ \end{cases}$

  * •

Draw $M$ i.i.d. samples
$\widetilde{w}_{j}^{(1)},\dots,\widetilde{w}_{j}^{(M)}$ from
$\mathsf{Unif}[0,1]$ and $\tilde{x}_{j}^{(1)},\dots,\tilde{x}_{j}^{(M)}$ from
$\mathcal{P}_{x}$.

  * •

Use score function $T$ to score the initial sample $(x_{j},w_{j})$ and its $M$
constructed counterfeits
$\\{(\tilde{x}_{j}^{(1)},\widetilde{w}_{j}^{(1)}),...,(\tilde{x}_{j}^{(M)},\widetilde{w}_{j}^{(M)})\\}$:

$\displaystyle T_{j}$ $\displaystyle=T(x_{j},w_{j})$
$\displaystyle\widetilde{T}_{j}^{(i)}$
$\displaystyle=T(\tilde{x}_{j}^{(i)},\widetilde{w}_{j}^{(i)})\,,\quad\text{for
}i\in[M]\,.$

  * •

Let $R_{j}$ denote the rank of $T_{j}$ among
$\\{T_{j},\widetilde{T}_{j}^{(1)},...,\widetilde{T}_{j}^{(M)}\\}$:

$R_{j}=1+\sum\limits_{i=1}^{M}{\mathbb{I}}_{\left\\{T_{j}\geq\widetilde{T}_{j}^{(i)}\right\\}}$

  * •

Assign label $L_{j}\in\\{1,2,...,L\\}$ to sample $j$ if $(L_{j}-1)K+1\leq
R_{j}\leq KL_{j}$.

end for

for  _$\ell\in\\{1,2,...,L\\}$_ do

  * •

Let $V_{\ell}$ be the number of samples with label $\ell$, i.e.
$V_{\ell}=\Big{|}\big{\\{}j\in\\{1,2,...,n\\}:L_{j}=\ell\big{\\}}\Big{|}\,.$

end for

Algorithm 3 Construction of GRASP test statistic (Model-X setting)

Decision rule. Let $\mathbf{V}_{n,L}$ be the statistic returned by Algorithm
3. We construct the test statistics
$U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})$ and
$U^{\mathsf{finite}}_{\tau}(\mathbf{V}_{n,L})$ following the same formulation
as in (5). The obtained test statistics are then used in (7) and (6) to give
the decision rules $\Phi_{n,L,\alpha,\tau}^{\mathsf{finite}}$ and
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$. In the next section, we show that
these rules control the size of our test
($\Phi_{n,L,\alpha,\tau}^{\mathsf{finite}}$ in finite sample settings and
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$ in asymptotic regimes).

###  4.2 Size of the model-X GRASP test

We first start with the following assumption, which assumes a density function
for the random variable $T(X,W)$. This is the unconditional version of
Assumption 3.3, and is used in our analysis of the test in a model-X setting.

######  Assumption 4.2.

Consider the following cumulative distribution functions:

$\displaystyle F_{o}(t)$ $\displaystyle={\mathbb{P}}(T(X,W)\leq t)\,,$
$\displaystyle F_{c}(t)$
$\displaystyle={\mathbb{P}}(T({\widetilde{X}},\widetilde{W})\leq t)\,,$

where $x,\tilde{x}\sim\mathcal{P}_{X}$ independently, $w$ is given by (11) and
$\widetilde{w}\sim\mathsf{Unif}[0,1]$. Assume that the density functions of
cdfs $F_{o}$ and $F_{c}$ exist and show them by $f_{o}(t)$ and $f_{c}(t)$,
respectively. In addition, consider a measure on the interval $[0,1]$ defined
with the cdf $\omega(u):=F_{o}(F_{c}^{-1}(u))$. Suppose that this measure is
absolutely continuous with respect to the Lebesgue measure, and denote its
density function (Radon–Nikodym derivative) by $\psi(u)$.

In our next result we show that the $f$-divergence between the multinomial
distribution of $\mathbf{V}_{n,L}$ (output of algorithm 3) and the uniform
multinomial distribution is bounded by the distance of the random variables
$T(x,w)$ and $T(\tilde{x},\widetilde{w})$.

######  Proposition 4.3.

Let $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]$ be outputs of Algorithm 3, then
$\mathbf{V}_{n,L}$ has a multinomial distribution with $L$ outcomes and
nominal probabilities $p_{1},\dots,p_{L}$. Under the setting of Assumption
4.2, for a score function $T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$, the
following holds

$\displaystyle\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})$ $\displaystyle\leq
D_{f}(\mathcal{L}(T(X,W))\|\mathcal{L}(T({\widetilde{X}},\widetilde{W})))\,.$

In addition, if the function $\psi(u)$ from Assumption 4.2 is continuous, the
above inequality becomes equality as $K,L$ grow to infinity:

$\lim_{L\to\infty}\lim_{K\to\infty}\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})=D_{f}(\mathcal{L}(T(X,W)),\mathcal{L}(T({\widetilde{X}},\widetilde{W})))\,.$

The proof of this proposition is given in Section A.10. The next Theorem
follows from the combination of Propositions 4.1 and 4.3, along with the data
processing inequality for $f$-divergence.

######  Theorem 4.4 (Model-X).

Let $\mathbf{V}_{n,L}=[V_{1},...,V_{L}]$ be the output of Algorithm 3 that has
a multinomial distribution with nominal probabilities $p_{1},\dots p_{L}$ as
per Proposition 4.3. Then for a score function
$T:\mathcal{X}\times[0,1]\to{\mathbb{R}}$ satisfying Assumption 4.2 the
following holds

$\displaystyle\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]\,.$

The proof of this theorem is given in Section A.11. The result of this theorem
allows us to perform statistical inference on the complex expression
${\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(x))||\mathsf{Bern}(\widehat{\eta}(x)))\right]$
by focusing on the deviation of $\mathbf{V}_{n,L}$ from the uniform
multinomial distribution, similar to the distribution-free version.

We are now ready to show that the two decision rules (7) and (6) control the
size of our testing in the model-X setup, respectively in the finite sample
and asymptotic regimes.

######  Theorem 4.5.

Let $\mathbf{V}_{n,L}=[V_{1},\dots,V_{L}]$ be the output of Algorithm 3.
Consider the decision rules $\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}$ and
$\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}$ which are given respectively by (6)
and (7), with test statistics $U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$ and
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ in (5). Then under the null
hypothesis (3), this holds:

${\mathbb{P}}\left(\Phi^{\mathsf{finite}}_{n,L,\alpha,\tau}=1\right)\leq\alpha\,,~{}~{}\lim\sup_{n\to\infty}{\mathbb{P}}\left(\Phi^{\mathsf{asym}}_{n,L,\alpha,\tau}=1\right)\leq\alpha\,.$

The proof of this theorem follows along the same lines as in the proof of
Theorem 3.6 .

###  4.3 Choice of the score function

In summary, our analysis in the previous section showed the following chain of
inequalities:

$\displaystyle\frac{1}{L}\sum_{\ell=1}^{L}f(Lp_{\ell})$ $\displaystyle\leq
D_{f}(\mathcal{L}(T(X,W))\|\mathcal{L}(T({\widetilde{X}},\widetilde{W})))$
$\displaystyle\leq
D_{f}(\mathcal{L}(X,W)\|\mathcal{L}({\widetilde{X}},\widetilde{W}))$
$\displaystyle={\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]\,,$
(12)

where the first inequality follows from Proposition 4.3, the second one is an
application of the data processing inequality, and the third line is the claim
proved in Proposition 4.1.

In order to increase the power of the test (make it less conservative), our
guide is to choose score functions for which the gap between the leftmost side
and the rightmost side in the chain of inequality is minimum.

For the first inequality, as shown in Proposition 4.3 (second part), the gap
vanishes for large values of $L$ and $K$. Focusing on the second inequality,
we seek score functions $T$ such that

$D_{f}(\mathcal{L}(T(X,W))\|\mathcal{L}(T({\widetilde{X}},\widetilde{W})))=D_{f}(\mathcal{L}(X,W)\|\mathcal{L}({\widetilde{X}},\widetilde{W}))\,.$
(13)

A trivial choice is $T(x,w)=(x,w)$, however, we would like to have scores with
one-dimensional range, since we need to compare the scores values of the
original samples and the counterfeits.

We proceed our discussion with providing a brief background on the variational
representation of $f$-divergences, with the following statement borrowed from
[NWJ10, Lemma 1].

######  Lemma 4.6.

For two probability density functions $q,p$ on
$\mathcal{X}\subset{\mathbb{R}}^{d}$ and a set of measurable functions
$\mathcal{G}$, we have

$D_{f}(p\|q)\geq\sup_{\varphi\in\mathcal{G}}\int(\varphi
p-f^{*}(\varphi)q){\rm d}\mu\,,$ (14)

where $f^{*}$ is the conjugate dual of $f$. In addition, the equality is
achieved if the subdifferential $\partial f(p/q)$ contains an element of
$\mathcal{G}$.

For the reader’s convenience, we provide the proof of Lemma 4.6 in Section
A.12.

######  Proposition 4.7.

For two probability density functions $p,q$ over
$\mathcal{X}\subset{\mathbb{R}}^{d}$ and function $f$ as per Definition 2.1,
let $\sigma\in\partial f(p/q)$. Then, for $X_{p}\sim p$ and $X_{q}\sim q$ the
$f$-divergence between distributions of $\sigma(X_{p})$ and $\sigma(X_{q})$ is
equal to the $f$-divergence of distributions $p$ and $q$. Formally,

$D_{f}(\mathcal{L}(\sigma(X_{p}))\|\mathcal{L}(\sigma(X_{q})))=D_{f}(p\|q)\,.$

The proof of this proposition is given in A.13.

Going back to our discussion on the choice of score function $T(x,w)$, observe
that by an application of Proposition 4.7, condition (13) is satisfied for
$T(x,w)\in\partial
f\left(p_{X,W}(x,w)/p_{{\widetilde{X}},\widetilde{W}}(x,w)\right)$. Since
$\widetilde{W}\sim\mathsf{Unif}[0,1]$, independent from
${\widetilde{X}}\sim\mathcal{P}_{X}$, we have

$\frac{p_{X,W}(x,w)}{p_{{\widetilde{X}},\widetilde{W}}(x,w)}=p_{W|X}(w|x)\,.$

In addition, from the construction of $w$, cf. Algorithm 3, we have

$\displaystyle{\mathbb{P}}(w|x)$
$\displaystyle={\mathbb{P}}(y=1|x){\mathbb{P}}(w|x,y=1)+{\mathbb{P}}(y=0|x){\mathbb{P}}(w|x,y=0)$
$\displaystyle=\frac{\eta(x)}{\widehat{\eta}(x)}{\mathbb{I}}(w\leq\widehat{\eta}(x))+\frac{1-\eta(x)}{1-\widehat{\eta}(x)}{\mathbb{I}}(\widehat{\eta}(x)\leq
w)\,.$ (15)

Putting things together, the optimal score function is given by

$T(x,w)\in\partial
f\left(\frac{\eta(x)}{\widehat{\eta}(x)}{\mathbb{I}}(w\leq\widehat{\eta}(x))+\frac{1-\eta(x)}{1-\widehat{\eta}(x)}{\mathbb{I}}(\widehat{\eta}(x)\leq
w)\right)\,.$

We next note that our test statistics in Algorithm 3 is based on the _relative
ranking_ of the score values, not the absolute values. For an increasing
function $g$, the rank values obtained by the score functions $T$ and $g\circ
T$ are the same. In particular, since $f$ is convex its derivative is a non-
decreasing function, and so we can simply consider the ratio of densities.
When $f$ is not strictly-convex, we may have ties in the ordering with respect
to the above derivation. We use the following rule which also breaks ties.

$T^{\mathsf{opt}}(x,w)=\frac{\eta(x)}{\widehat{\eta}(x)}{\mathbb{I}}(w\leq\widehat{\eta}(x))+\frac{1-\eta(x)}{1-\widehat{\eta}(x)}{\mathbb{I}}(\widehat{\eta}(x)\leq
w)\,.$ (16)

As can be seen $T^{\mathsf{opt}}(x,w)$ involves the true model $\eta(x)$ which
is unknown. To cope with this issue, we consider two routes:

  * •

Model-agnostic approach: We use $T^{\mathsf{opt}}(x,w)$ with replacing
$\eta(x)=\frac{1}{2}$.

  * •

GAN-based approach: We propose an optimization inspired by generative
adversarial networks (GANs) to approximate the ratio of densities
$p_{X,W}(x,w)/p_{{\widetilde{X}},\widetilde{W}}(x,w)$.

Our next subsection gives a brief background overview on GANs and the details
of our GAN-based approach.

####  4.3.1 Generative Adversarial Nets (GANs)

The GANs framework was introduced by [GPAM+14] to generate samples from the
population of given data samples $x\in\mathcal{X}\subset{\mathbb{R}}^{d}$.
This framework can be perceived as a game between a generator that tries to
learn data distribution $p_{X}$ and produce fake samples that are
statistically close to the original data, and a discriminator that tries to
discriminate samples of the generator from the original data samples.

It is a repeated game where each player aims to improve her model with respect
to her objective. Starting from noise sample $z\in{\cal
Z}\subset{\mathbb{R}}^{k}$ with density function $p_{Z}$, the generator
produce fake samples in the data space $\mathcal{X}$ via a mapping
$G(z;\theta_{g})$, where $G:{\cal Z}\to\mathcal{X}$ is a differentiable neural
network with parameters $\theta_{g}$. The discriminator, on the other side,
has access to another neural network $D(x;\theta_{d})$ that represents the
likelihood of data $x$ coming from the original data distribution $p_{X}$.
This game can be cast as the following min-max optimization problem

$\min_{\theta_{g}}\max_{\theta_{d}}\left\\{{\mathbb{E}}_{X}\Big{[}\log
D(X;\theta_{d})\Big{]}+{\mathbb{E}}_{Z}\Big{[}\log\big{(}1-D(G(Z;\theta_{g});\theta_{d})\big{)}\Big{]}\right\\}\,.$
(17)

As evident from the above formulation, the discriminator aims to increase the
likelihood of the original samples coming from $p_{X}$, while simultaneously
lower the likelihood of fake sample $G(z;\theta_{g})$ coming from the same
distribution. The generator, on the other hand, aims to find a model
$G(z;\theta_{g})$ which generates samples with high likelihood of coming from
$p_{X}$. As shown in [GPAM+14], the above optimization problem has a global
optimum at $p_{g}=p_{X}$, where $p_{g}$ stands for the distribution of samples
$G(z;\theta_{g})$ with $z\sim p_{Z}$.

We follow the GANs framework by considering $\mathcal{L}(X,W)$ as the original
data distribution. In addition, we assume that the generator at each round
produces the randomizations $(\tilde{x},\widetilde{w})$, and update the
discriminator to discriminate $(x,w)$ from $(\tilde{x},\widetilde{w})$. The
major distinction with the original GAN setup is dropping a separate network
for the generator as the generated samples are always coming from
$\mathcal{L}({\widetilde{X}},\widetilde{W})$. Optimization problem (17) then
reduces to:

$\max_{\theta_{d}}\left\\{{\mathbb{E}}_{(X,W)}\Big{[}\log
D((X,W);\theta_{d})\Big{]}+{\mathbb{E}}_{({\widetilde{X}},\widetilde{W})}\Big{[}\log\big{(}1-D(({\widetilde{X}},\widetilde{W});\theta_{d})\big{)}\Big{]}\right\\}\,.$
(18)

The next result connects the optimal discriminator $D(x)$ (solution of (18))
and the likelihood ratio function. This proposition is borrowed from [GPAM+14,
Proposition 1].

######  Proposition 4.8.

Let $D^{*}(x,w)$ be the maximizer of optimization problem (18). Then we have

$D^{*}(x,w)=\frac{p_{X,W}(x,w)}{p_{X,W}(x,w)+p_{{\widetilde{X}},\widetilde{W}}(x,w)}\,.$

Since we only have one network, we only update the parameters of the
discriminative network $D$. For this end, we adopt the same optimization
procedure used in [GPAM+14, NCT16], where at each step given a mini-batch of
samples $\left\\{(x^{1},w^{1}),\dots,(x^{m},w^{m})\right\\}$ and
$\left\\{(\tilde{x}^{1},\widetilde{w}^{1}),\dots,(\tilde{x}^{m},\widetilde{w}^{m})\right\\}$
we use the following stochastic gradient ascent:

$\theta_{d}^{(t+1)}=\theta_{d}^{(t)}+\frac{\beta_{t}}{m}\sum_{i=1}^{m}\nabla_{\theta_{d}}\left[\log
D\Big{(}(x^{i},w^{i});\theta_{d}^{(t)}\Big{)}+\log\Big{(}1-D\big{(}(\tilde{x}^{i},\widetilde{w}^{i});\theta_{d}^{(t)}\big{)}\Big{)}\right]\,,$

with step size $\beta_{t}$. The gradients of the above optimization problem
can be efficiently computed using back-propagation algorithms. Let
${\widehat{\theta}}_{d}$ be the resulting estimate after convergence.
Considering Proposition 4.8 and our discussion below Proposition 4.7, we
consider the following score function

$T(x,w)=\frac{D((x,w),{\widehat{\theta}}_{d})}{1-D((x,w),{\widehat{\theta}}_{d})}\,.$

##  5 Numerical Experiments

In this section, we evaluate the performance of our proposed methodology on
various synthetic datasets. Consider a binary classification problem under the
logistic regression setting. Throughout the experiments, we let the feature
vectors have an isotropic Gaussian distribution and the conditional
probability be given by

$\eta(x)=\frac{1}{1+\exp(-x^{\sf
T}\theta_{0})}\,,~{}x\sim\mathsf{N}(0,I_{d})\,.$ (19)

We construct $\theta_{0}$ by drawing one time realization from the
$\mathsf{N}(0,\sigma^{2}I_{d})$ distribution and this $\theta_{0}$ value is
fixed for the rest of the experiments. We set the feature dimension $d=200$
and standard deviation $\sigma=0.25$. Let $\widehat{\eta}(x)$ denote the model
that we want to perform the goodness-of-fit testing. We consider
$\widehat{\eta}(x)$ to be another logistic model with parameter $\theta_{1}$.
Concretely, $\widehat{\eta}(x)={1/[1+\exp(-x^{\sf T}\theta_{1})]}$. Obviously
$\widehat{\eta}(x)=\eta(x)$ if $\theta_{1}=\theta_{0}$.

######  Example 5.1.

( Size of the test) _We begin with studying the size of our proposed test
statistics. We focus on the $\theta_{1}=\theta_{0}$ setting, which implies
that
${\mathbb{E}}[D_{f}(\mathsf{Bern}(\eta(X))||\mathsf{Bern}(\widehat{\eta}(X))]=0$
for every divergence function $f$. We set $\tau=0$ in the null hypothesis (3)
and consider three significance levels $\alpha=0.05,0.1$ and $0.15$. The
rejection rate of $H_{0}$ is reported as the average over $200$ independent
experiments. Since $\tau=0$ the optimization problem in (5) is independent of
the choice of $f$. Its only feasible solution is given by
$p_{1}=\dots=p_{L}=1/L$, and so the test statistics amount to_

$U_{0}^{\mathsf{asym}}(\mathbf{V}_{n,L})=\frac{L}{n}\sum_{\ell=1}^{L}\left(\mathbf{V}_{\ell}-\frac{n}{L}\right)^{2}\,,\quad
U_{0}^{\mathsf{finite}}(\mathbf{V}_{n,L})=\frac{2L}{n}\sum_{\ell=1}^{L}\left(\mathbf{V}_{\ell}-\frac{n}{L}\right)^{2}\,.$

_Table 1 summarizes the rejection rates for the distribution-free Algorithm 2
and for different values of $(n,L)$. Each cell of the table has an ordered
pair, where the first and second entry respectively corresponds to the
asymptotic decision rule (6) and the finite decision rule (7)._

_We next consider model-X setting for a similar experiment. We consider the
model-agnostic approach in choosing the score function $T$ discussed at the
end of Section 4.3, namely,_

$T(x,w)={\mathbb{I}}(w\leq\widehat{\eta}(x))~{}\frac{1}{2\widehat{\eta}(x)}+{\mathbb{I}}(\widehat{\eta}(x)\leq
w)~{}\frac{1}{2(1-\widehat{\eta}(x))}\,.$ (20)

_Table 2 presents the rejection rate of the model-X Algorithm 3 for three
significance levels $\alpha=0.05,0.1,0.15$ for different number of samples
($n$), number of labels ($L$) and randomizations per label $(K)$. Each cell of
Tables 2 consists of an ordered pair, where the first and second numbers
respectively correspond to the asymptotic decision rule (6) and the finite
decision rule (7). Reported numbers are averaged over $200$ independent
experiments. _

| $\alpha=0.05$ | $\alpha=0.1$ | $\alpha=0.15$
---|---|---|---
$n$ $L$ | $50$ | $100$ | $50$ | $100$ | $50$ | $100$
5000 | (0.055,0) | (0.03,0) | (0.105,0) | (0.075,0) | (0.16,0) | (0.1,0)
20000 | (0.07,0) | (0.05,0) | (0.1,0) | (0.1,0) | (0.185,0) | (0.16,0)
50000 | (0.06,0) | (0.05,0) | (0.105,0) | (0.105,0) | (0.15,0) | (0.145,0)

Table 1: Size of the distribution-free GRASP, cf. Algorithm 2. In this
experiment, $\widehat{\eta}(x)=\eta(x)$ is considered with $\eta(x)$ being a
logistic model (19). We consider both decision rules (asymptotic (6)) and
(finite (7)). The first and the second entry of each pair in the table
represent the rejection rate of these decision rules, respectively. The
numbers are averaged over $200$ independent experiments for three significance
levels $\alpha=0.05,0.1,0.15$ and for different sample sizes ($n$) and number
of labels ($L$).

| $\alpha=0.05$ | $\alpha=0.1$ | $\alpha=0.15$
---|---|---|---
$(n,K)$ $L$ | $50$ | $100$ | $50$ | $100$ | $50$ | $100$
(5000,1) | (0.055,0) | (0.04,0) | (0.115,0) | (0.08,0) | (0.165,0) | (0.11,0)
(5000,5) | (0.025,0) | (0.07,0) | (0.12,0) | (0.12,0) | (0.185,0) | (0.125,0)
(10000,1) | (0.085,0) | (0.035,0) | (0.125,0) | (0.095,0) | (0.18,0) | (0.105,0)
(10000,5) | (0.065,0) | (0.06,0) | (0.07,0) | (0.095,0) | (0.17,0) | (0.125,0)

Table 2: Size of the model-X GRASP, cf. Algorithm 3. In this experiment,
$\widehat{\eta}(x)=\eta(x)$ is considered with $\eta(x)$ being a logistic
model (19). Algorithm 3 is run with the score function $T(x,w)$ given in (20).
We consider both decision rules (asymptotic (6)) and (finite (7)). The first
and the second entry of each pair in the table represent the rejection rate of
these decision rules, respectively. The numbers are averaged over $200$
independent experiments for three significance levels $\alpha=0.05,0.1,0.15$
and for different sample sizes ($n$) and number of labels ($L$), and number of
randomizations per label ($K$).

######  Example 5.2.

(Power of the test) _In this experiment, we consider the null hypothesis ( 3)
for a logistic model $\widehat{\eta}(x)$ with parameter
$\theta_{1}=-\theta_{0}$. Formally, we have_

$\widehat{\eta}(x)=\frac{1}{1+\exp(-\theta_{1}^{\sf T}x)}\,,\text{ with
}\theta_{1}=-\theta_{0}\,.$ (21)

_In this experiment, we consider three divergence functions (Kl, TV, and
Hellinger given in Lemma 2.2) between the true model $\eta(\cdot)$ and the
test model $\widehat{\eta}(\cdot)$. The computed nominal values of these
divergences are $\tau_{0}^{\mathsf{KL}}=2.7819,\tau_{0}^{\mathsf{TV}}=0.7330,$
and $\tau_{0}^{\mathsf{H}}=0.9576$. More precisely, we have_

$\displaystyle\tau_{0}^{\mathsf{KL}}$
$\displaystyle={\mathbb{E}}\left[D_{\mathsf{KL}}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,,$
(22) $\displaystyle\tau_{0}^{\mathsf{TV}}$
$\displaystyle={\mathbb{E}}\left[D_{\mathsf{TV}}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,,$
(23) $\displaystyle\tau_{0}^{\mathsf{H}}$
$\displaystyle={\mathbb{E}}\left[D_{\mathsf{H}}\Big{(}\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\Big{)}\right]\,.$
(24)

_In order to analyze the power of proposed methods, we choose $\tau$ that are
smaller than $\tau_{0}$ values, consequently we expect that the null
hypothesis (3) must be rejected. For each divergence function, we consider
four $\tau$ values and report the average number of times (out of $50$
independent experiments) that the proposed method rejects the null hypothesis
$H_{0}$. Table 3 reports the results for distribution-free GRASP, outlined in
Algorithm 2, for a variety of choices of $(n,L)$. Each cell of Table 3
consists of an ordered pair where the first and second coordinate respectively
stands for the asymptotic decision rule (6) and the finite decision rule (7)
at significance level $\alpha=0.1$. _

_To analyze the statistical power of model-X GRASP, as outlined in Algorithm
3, we adopt a similar setup with the test model $\widehat{\eta}(x)$ given by
(21), and run the model-X GRASP with the score function (20). Table 4 presents
the results for different values of the sample size ($n$), number of labels
($L$), and randomizations per label $(K)$, at the significance level
$\alpha=0.1$. Each cell of Table 4 has an ordered pair with the first numbers
indicates the rejection rate with the asymptotic decision rule (6) and the
second number indicates the rejection rate with the finite decision rule (7).
Reported numbers are averaged over $50$ independent experiments. _

_Since the problem setting for the last two power analysis experiments
(distribution-free GRASP and model-X GRASP ) are the same, by a simple
comparison between Table 3 and Table 4, it can be observed that for each of
the $f$-divergence functions, the model-X GRASP algorithm achieves a higher
statistical power than the distribution-free GRASP procedure. This is expected
as the model-X GRASP utilizes the covariate distribution to obtain higher
statistical power. _

| $\mathsf{KL}$ | $\mathsf{TV}$ | $\mathsf{Hellinger}$
---|---|---|---
| $\tau_{0}^{\mathsf{KL}}=2.7819$ | $\tau_{0}^{\mathsf{TV}}=0.7330$ | $\mathsf{\tau_{0}^{\mathsf{H}}}=0.9576$
$(n,L)$ $\tau$ | $0.72$ | $0.82$ | $0.96$ | $1.02$ | $0.4$ | $0.44$ | $0.48$ | $0.52$ | $0.28$ | $0.32$ | $0.36$ | $0.4$
$(5000,50)$ | (1,1) | (1,0.02) | (0,0) | (0,0) | (1,1) | (1,0.52) | (0.18,0) | (0,0) | (1,1) | (1,0.08) | (0.08,0) | (0,0)
$(20000,50)$ | (1,1) | (1,1) | (0,0) | (0,0) | (1,1) | (1,1) | (1,0.86) | (0,0) | (1,1) | (1,1) | (1,0.04) | (0,0)
$(50000,50)$ | (1,1) | (1,1) | (0.54,0) | (0,0) | (1,1) | (1,1) | (1,1) | (0,0) | (1,1) | (1,1) | (1,0.04) | (0,0)
$(5000,100)$ | (1,1) | (1,0) | (0,0) | (0,0) | (1,1) | (1,0) | (0.02,0) | (0,0) | (1,0.86) | (1,0) | (0.06,0) | (0,0)
$(20000,100)$ | (1,1) | (1,1) | (1,0) | (0,0) | (1,1) | (1,1) | (1,0) | (0,0) | (1,1) | (1,1) | (1,0) | (0,0)
$(50000,100)$ | (1,1) | (1,1) | (1,1) | (0.44,0) | (1,1) | (1,1) | (1,1) | (0,0) | (1,1) | (1,1) | (1,1) | (0.06,0)

Table 3: Statistical power of distribution-free GRASP, as outlined in
Algorithm 2, with the true model $\eta(x)$ given by (19) and the test model
$\widehat{\eta}(x)$ in (21). In this experiment, we consider three
$f-$divergences (KL, TV, and Hellinger as per Lemma 2.2) for different choices
of sample size ($n$) and number of labels ($L$). Each cell has two numbers,
where the first and second number respectively refer to the rejection rates
under the asymptotic decision rule (6) and the finite decision rule (7), at
significance level $\alpha=0.1$. The reported numbers are averaged over $50$
independent experiments

| $\mathsf{KL}$ | $\mathsf{TV}$ | $\mathsf{Hellinger}$
---|---|---|---
| $\tau_{0}=2.7819$ | $\tau_{0}=0.7330$ | $\mathsf{\tau_{0}=0.9576}$
$(n,K,L)$ $\tau$ | $1.5$ | $1.7$ | $1.9$ | $2$ | $0.6$ | $0.63$ | $0.66$ | $0.7$ | $0.58$ | $0.62$ | $0.70$ | $0.8$
$(5000,5,50)$ | (1,1) | (1,0.12) | (0.82,0) | (0,0) | (1,1) | (1,1) | (1,0) | (0.04,0) | (1,1) | (1,0.98) | (1,0) | (0.5,0)
$(5000,5,100)$ | (1,1) | (1,0.38) | (1,0) | (0.28,0) | (1,1) | (1,0) | (1,0) | (0,0) | (1,1) | (1,0) | (1,0) | (0.02,0)
$(10000,5,50)$ | (1,1) | (1,1) | (1,0) | (0.1,0) | (1,1) | (1,1) | (1,1) | (0.84,0) | (1,1) | (1,1) | (1,0) | (1,0)
$(10000,5,100)$ | (1,1) | (1,1) | (1,0) | (1,0) | (1,1) | (1,1) | (1,0) | (0,0) | (1,1) | (1,1) | (1,0) | (1,0)
$(5000,1,50)$ | (1,1) | (1,0) | (0,0) | (0,0) | (1,1) | (1,0.78) | (1,0) | (0,0) | (1,1) | (1,0) | (1,0) | (0,0)
$(5000,1,100)$ | (1,1) | (1,0) | (0.38,0) | (0,0) | (1,0.96) | (1,0) | (0.96,0) | (0,0) | (1,0.44) | (1,0) | (1,0) | (0,0)
$(10000,1,50)$ | (1,1) | (1,0.1) | (0.02,0) | (0,0) | (1,1) | (1,1) | (1,0.38) | (0.04,0) | (1,1) | (1,1) | (1,0) | (0.74,0)
$(10000,1,100)$ | (1,1) | (1,1) | (1,0) | (0.06,0) | (1,1) | (1,1) | (1,0) | (0,0) | (1,1) | (1,1) | (1,0) | (0.76,0)

Table 4: Statistical power of model-X GRASP, as outlined in Algorithm 3, with
the true model $\eta(x)$ given by (19) and the test model $\widehat{\eta}(x)$
in (21). We consider three $f-$divergences (KL, TV, and Hellinger as per Lemma
2.2) for different choices of sample size ($n$) and number of labels ($L$),
and number of randomizations per label ($K$). The model-X algorithm is run
with the score function $T(x,w)$ given in (20). Each cell contains an ordered
pair, where the first and second number respectively refer to the rejection
rates under the asymptotic decision rule (6) and the finite decision rule (7),
at significance level $\alpha=0.1$. The reported numbers are averaged over
$50$ independent experiments.

######  Example 5.3.

(Score function for distribution-free GRASP) _In this experiment, we evaluate
the performance of distribution-free GRASP with other choices of score
function. So far, we have used the score function $T(x,w)=w$. However, we want
to examine whether one can achieve higher statistical power by selecting a
score function that accounts for variability between $x$ and $w$. To this end,
we perform a similar experiment as the one reported in Table 3 for
$(n,L)=(5000,50)$, but with a new score function. We fit a linear model with
response value $w$ and feature vectors $x$ on an auxiliary dataset of size
$4000$. We denote the fitted linear model by $\hat{\theta}$ and consider the
score function $T(x,w)=|w-x^{\sf T}{\widehat{\theta}}|$. We report the results
for three f-divergence functions (TV, KL-divergence, and Hellinger distance)
for both asymptotic and finite-sample decision rules in Figure 1. We averaged
the reported numbers over $50$ experiments at a significance level of
$\alpha=0.1$. By comparing the rejection rates with the reported numbers in
the first row of Table 3, we find that the previous score function $T(x,w)=w$
performs slightly better than the fitted model. For instance, for KL
divergence, using $T(x,w)=w$ with both asymptotic and finite decision rules
achieves full power (rejection rate 1.00) at $\tau=0.72$, while using the
fitted score function already achieves trivial power (rejection rate $0$) at
$\tau=0.7$. Similar observations can be made for other divergence functions as
well. Likewise, for TV and Hellinger divergence, we also observe that score
function $T(x,w)=w$ outperforms the regression-based score function slightly._

Figure 1: Statistical power of distribution free GRASP for trained score
function $T(x,w)=|w-x^{\sf T}{\widehat{\theta}}|$. The linear model
${\widehat{\theta}}$ is trained on an auxiliary dataset of size $4000$. The
true model $\eta(x)$ is given in (19) and the test model $\widehat{\eta}(x)$
is given in (21). In this case, the nominal values are
$\tau_{0}^{\mathsf{KL}}=2.7819$, $\tau_{0}^{\mathsf{TV}}=0.7330,$ and
$\tau_{0}^{\mathsf{H}}=0.9576$ (see (22)). For three $f-$divergence (KL, TV,
Hellinger) we run distribution free GRASP with two decision rules
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$ and
$\Phi_{n,L,\alpha,\tau}^{\mathsf{finite}}$ given in (6) and (7). For each
setting, we consider the null-hypothesis (3) with different $\tau$ values.
Here, sample size is $n=5000$ with number of labels $L=50$, and the
significance level is $\alpha=0.1$. Reported rates are averaged over $50$
experiments.

######  Example 5.4.

(The choice of number of labels $L$) _In this experiment, we analyze the
statistical power of the proposed methods for goodness-of-fit testing, while
allowing a wide range of values for the inside parameter $L$ ( number of
labels). We consider the previous setting with the true model $\eta(x)$ as in
(19) and the test model $\widehat{\eta}(x)$ as in (21). We consider the null
hypothesis (3) for three $f$ divergence functions $\mathsf{KL},\mathsf{TV}$
and $\mathsf{H}$. In addition, the two decision rules (6) and (7) are being
deployed. In this setting, the nominal values of divergence functions are
$\tau_{0}^{\mathsf{KL}}=2.7819$, $\tau_{0}^{\mathsf{TV}}=0.7330,$ and
$\tau_{0}^{\mathsf{H}}=0.9576$, where $\tau_{0}$ values are given in (22). _

_We consider the same range of $L$ values for different $f$-divergence (KL,
TV, Hellinger) and decision rules ($\Phi_{n,L,\tau}^{\mathsf{asym}}$ and
$\Phi_{n,L,\tau}^{\mathsf{finite}}$). However, we choose different values of
$\tau$ in the hypothesis (3), depending on the setting to better capture the
effect of $L$. For distribution-free GRASP we consider the null hypothesis (3)
with
$\tau_{\mathsf{asym}}^{\mathsf{KL}}=0.84,\tau_{\mathsf{finite}}^{\mathsf{KL}}=0.73,\tau_{\mathsf{asym}}^{\mathsf{TV}}=0.46,\tau_{\mathsf{finite}}^{\mathsf{TV}}=0.42,\tau_{\mathsf{asym}}^{\mathsf{H}}=0.34,\tau_{\mathsf{finite}}^{\mathsf{H}}=0.30$.
Figure 2 depicts the statistical power curves under these settings versus the
number of labels $L$. The sample size is set to $n=5000$ and the significance
level to $\alpha=0.1$. The obtained numbers are averaged over $50$ independent
experiments._

_We repeat the same experiment for model-X GRASP with the score function (
20). Here, we consider
$\tau_{\mathsf{asym}}^{\mathsf{KL}}=1.87,\tau_{\mathsf{finite}}^{\mathsf{KL}}=1.6,\tau_{\mathsf{asym}}^{\mathsf{TV}}=0.67,\tau_{\mathsf{finite}}^{\mathsf{TV}}=0.62,\tau_{\mathsf{asym}}^{\mathsf{H}}=0.75$,
and $\tau_{\mathsf{finite}}^{\mathsf{H}}=0.6$. The sample size is set to
$n=5000$ with randomizations per label $K=1$. The significance level is set to
$\alpha=0.1$ and the reported numbers are averaged over $50$ independent
experiments. The behavior of the statistical power with respect to the change
in the number of labels ($L$) can be seen in Figure 3. As observed from the
curves in Figures 2 and 3, the power favors a middle range of $L$ values
before and after which the power starts to decline. _

Figure 2: Statistical power of distribution free GRASP versus the number of
labels $L$. The true model $\eta(x)$ is given in (19) and the test model
$\widehat{\eta}(x)$ is given in (21). In this case, the nominal values are
$\tau_{0}^{\mathsf{KL}}=2.7819$, $\tau_{0}^{\mathsf{TV}}=0.7330,$ and
$\tau_{0}^{\mathsf{H}}=0.9576$ (see (22)). For three $f-$divergence (KL, TV,
Hellinger) we run distribution free GRASP with two decision rules
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$ and
$\Phi_{n,L,\alpha,\tau}^{\mathsf{finite}}$ given in (6) and (7). For each
setting, we consider the null-hypothesis (3) with the following $\tau$ values:
$\tau_{\mathsf{asym}}^{\mathsf{KL}}=0.84,\tau_{\mathsf{finite}}^{\mathsf{KL}}=0.73,\tau_{\mathsf{asym}}^{\mathsf{TV}}=0.46,\tau_{\mathsf{finite}}^{\mathsf{TV}}=0.42,\tau_{\mathsf{asym}}^{\mathsf{H}}=0.34$,
and $\tau_{\mathsf{finite}}^{\mathsf{H}}=0.30$. Here, sample size is $n=5000$,
the significance level is $\alpha=0.1$ and the reported rates are averaged
over $50$ experiments.  Figure 3: Statistical power of model-X GRASP versus
the number of labels $L$. The true model $\eta(x)$ is given in (19) and the
test model $\widehat{\eta}(x)$ is given in (21). In this case, the nominal
values are $\tau_{0}^{\mathsf{KL}}=2.7819$, $\tau_{0}^{\mathsf{TV}}=0.7330,$
and $\tau_{0}^{\mathsf{H}}=0.9576$ (see (22)). For three $f-$divergence (KL,
TV, Hellinger) we run model-X GRASP with the score function $T(x,w)$ as per
(20) and the two decision rules $\Phi_{n,L,\tau}^{\mathsf{asym}}$ and
$\Phi_{n,L,\tau}^{\mathsf{finite}}$ given in (6) and (7). For each setting, we
consider the null-hypothesis (3) with the following $\tau$ values:
$\tau_{\mathsf{asym}}^{\mathsf{KL}}=1.87,\tau_{\mathsf{finite}}^{\mathsf{KL}}=1.6,\tau_{\mathsf{asym}}^{\mathsf{TV}}=0.67,\tau_{\mathsf{finite}}^{\mathsf{TV}}=0.62,\tau_{\mathsf{asym}}^{\mathsf{H}}=0.75$,
and $\tau_{\mathsf{finite}}^{\mathsf{H}}=0.6$. Here, sample size is $n=5000$,
with $K=1$ randomizations per sample, and the significance level $\alpha=0.1$.
The reported rates are averaged over $50$ experiments.

######  Example 5.5.

(Statistical evidence for misfitted models) _In this experiment, we test for
perfect fit of models to the underlying conditional law, by considering
hypothesis testing problem ( 3) with $\tau=0$. We consider feature vectors
drawn independently from $\mathsf{N}(0,I_{d})$, where $d=300$. We use a
planted model setup, where the data generating law is formulated by a two-
layer neural network. Specifically, we consider a fully-connected neural
network with $N=400$ ReLU neurons in the first layer and a sigmoid activation
function as the output unit. The conditional distribution $\mathcal{L}(Y|X)$
is given by:_

$Y|X\sim\mathsf{Bern}\left({\mathsf{Sigmoid}}\big{(}\theta_{0}^{\sf
T}{\mathsf{ReLU}}(W_{0}X)\big{)}\right)\,,$ (25)

_where $W_{0}\in{\mathbb{R}}^{N\times d}$ is the weight matrix,
$\theta_{0}\in{\mathbb{R}}^{N}$ is the weight vector of the second layer, and
${\mathsf{ReLU}}$ is the ReLU activation function applied entry-wise. The
entries of the planted parameters $(W_{0},\theta_{0})$ are drawn independently
from $\mathsf{N}(0,1)$._

_We compare two neural networks, $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$,
trained on a similar dataset ${\cal
D}^{\mathsf{train}}=\\{(x_{i},y_{i})\\}_{i=1:T}$, but with different weight
initialization schemes. Specifically, we initialize the weights of
$\mathcal{N}_{1}$ using the Kaiming initialization method, while the weights
of $\mathcal{N}_{2}$ are initialized by adding independent noise from a
standard normal distribution with mean 0 and standard deviation $10^{-2}$ to
each entry of $W_{0}$ and $\theta_{0}$, We train both models on a training
dataset of size $T=2^{16}$ (approximately 65K samples), and evaluate them on a
test dataset of size 128._

_We first compare the predictive performance of $\mathcal{N}_{1}$ and
$\mathcal{N}_{2}$ on test datasets. For this end, we consider $500$
independent draws of evaluation sets of size $128$ and report the average test
accuracy of models on the drawn datasets. Figure 4(a) shows the boxplots of
computed test accuracies for both models $\mathcal{N}_{1}$ and
$\mathcal{N}_{2}$ over all $500$ experiments. The boxplots show that
$\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ have very similar accuracies, despite
the fact that the $\ell_{2}$-distance between the model parameters of
$\mathcal{N}_{1}$, and $\mathcal{N}_{2}$ from the ground-truth model (25) are
significantly different (362.12 for $\mathcal{N}_{1}$ and 80.87 for
$\mathcal{N}_{2}$). This indicates that the perturbed initialization of
$\mathcal{N}_{2}$ has resulted in a better local minimum in the optimization
landscape compared to $\mathcal{N}_{1}$, despite their similar empirical
accuracies, and therefore $\mathcal{N}_{2}$ is a better fit to
$\mathcal{L}(Y|X)$ than $\mathcal{N}_{1}$. This further highlights the fact
that solely focusing on test accuracy is not sufficient to differentiate
between the models. _

_We evaluate the performance of GRASP to determine the significance of
evidence against the assumption that $\mathcal{N}_{1},\mathcal{N}_{2}$ are
perfectly fitted to $\mathcal{L}(Y|X)$. We construct the $p$-values for the
two models according to (10), with $M=500$. For the score function, we fit a
three-layer neural network to regress $w$ on $\mathbf{x}$, namely
$\mathbf{v}^{\sf
T}\mathsf{ReLU}(\mathbf{A}_{2}\mathsf{ReLU}(\mathbf{A}_{1}\mathbf{x}))$, with
$\mathbf{v}\in{\mathbb{R}}^{100\times 1}$,
$\mathbf{A}_{2}\in{\mathbb{R}}^{100\times 300}$,
$\mathbf{A}_{1}\in{\mathbb{R}}^{300\times 300}$. We define score function
$T(\mathbf{x},\mathbf{w})$ as the MSE of this fitted model over the evaluation
data set of size 128. For the $500$ independent draws of the evaluation sets
of size $128$, we use the same datasets used earlier for test accuracies, and
compute the GRASP p-values for $\mathcal{N}_{1}$, $\mathcal{N}_{2}$, as well
as the nominal model ($\mathcal{L}(Y|X)$). Figures 4(b), 4(c), 4(d) depict the
Q-Q plot of GRASP p-values (across the 500 experiments) for models
$\mathcal{N}_{1}$, $\mathcal{N}_{2}$, and the nominal model
$\mathcal{L}(Y|X)$, respectively. It can be observed that the GRASP p-values
for the nominal model follow a uniform distribution, while for the other two
models the Q-Q plots deviate from the uniform distribution. In addition, the
$\mathcal{N}_{2}$ p-values have higher deviation from uniform distribution
than the $\mathcal{N}_{1}$ p-values, and this aligns with our expectation
given that $\mathcal{N}_{2}$ is closer to the ground-truth model in $\ell_{2}$
distance than $\mathcal{N}_{1}$. By hypothesis testing (3) (for $\tau=0$) at
the significance level $\alpha=0.1$, the nominal model has rejection rate
$9.8\%$, and models $\mathcal{N}_{1},\mathcal{N}_{2}$ have average rejection
rates of $49.2\%,21.2\%$, respectively. _

(a) Boxplot for test accuracies of $\mathcal{N}_{1}$ and $\mathcal{N}_{2}$

(b) QQ plot of GRASP p-values for $\mathcal{N}_{1}$

(c) QQ plot of GRASP p-values for $\mathcal{N}_{2}$

(d) QQ plot of GRASP p-values for nominal model (25)

Figure 4:  Test accuracies and GRASP $p$-values for the two models
$\mathcal{N}_{1}$ and $\mathcal{N}_{2}$ over $500$ draws of evaluation
datasets of size $128$ (under the set-up of Example 5.5). While both models
have very similar test accuracies, we observe that the GRASP $p$-values of
$\mathcal{N}_{1}$ show a significanlty higher deviation from the uniform
distribution, compared to the GRASP $p$-values of $\mathcal{N}_{2}$. For
testing the perfect fit ($\tau=0$), at the significance level $\alpha=0.1$,
the nominal model has rejection rate $9.8\%$, and models
$\mathcal{N}_{1},\mathcal{N}_{2}$ have average rejection rates of
$49.2\%,21.2\%$, respectively.

######  Example 5.6.

(GAN-based versus model-agnostic score functions) _In this experiment, we
compare the performance of model-X GRASP algorithm for two choices of score
functions. Formally, the considered score functions are GAN-based approach and
the model-agnostic given in ( 20). We consider the hypothesis testing problem
(3) for three different $f$-divergence functions TV, KL and the Hellinger
distance, and two decision rules, asymptotic (6) and finite (7). For each
setting (6 overall), we run the experiment with one of the score functions (12
experiments in total). The true model $\eta(x)$ is given by (19) and the test
model $\widehat{\eta}(x)$ is given by (21), with $\theta_{0}$ a draw of
$\mathsf{N}(0,I_{d})$, the feature dimension $d=200$ and
$\theta_{1}=-2\theta_{0}$. We let the number of samples be $n=5000$, number of
labels $L=50$, and $K=1$ randomizations per label. Figure 5 depict the
performance of model-X GRASP for GAN based score function and agnostic (AGN)
score function for three $f$-divergence functions, and the two decision rules
$\mathsf{asym}$ and $\mathsf{finite}$. The reported rates are averaged over
$50$ instances. It is observed that in each setting, the GAN-based score
function achieves a higher statistical power for larger $\tau$ values than its
counterpart used with an agnostic score function, for both the asymptotic and
the finite decision rules. _

(a) KL

(b) TV

(c) Hellinger

Figure 5: Comparison between the statistical power of model-X GRASP with a
GAN-based score function and the agnostic score function (20). In this
experiment, the test size is $n=5000$, number of labels $L=50$, and
randomizations per label $K=1$. Hypothesis testing problem (3) is considered
for a variety of $\tau$ values, three $f$-divergence functions, and two
asymptotic and finite decision rules (6) and (7). In each setting, it can be
observed that the GAN-based score function achieves a higher statistical power
than the agnostic score function. The reported numbers are average over $50$
experiments.

_For the discriminator of the GAN-based score, we use a neural network with
input dimension $201$ (as inputs are $(x,w)$ with $x\in{\mathbb{R}}^{200}$ and
$w\in{\mathbb{R}}$). The first hidden layer consists of $256$ neurons with
ReLU activation function. The second and third hidden layers have respectively
128 and 64 neurons with ReLU activations. The output is a single neuron with
sigmoid activation function. The cross-entropy loss is considered for the
training process, and the network is trained on $64$k inputs. In order to
prevent overfitting, the dropout probability 0.3 is considered for neural
connections. _

##  6 Solving the optimization problem for GRASP decision rules

The optimization problems (5) have convex objectives and constraints
($f$-divergence ball), and therefore we have a convex optimization problem. We
use an iterative optimization procedure to find the optimal solution. It is
worth highlighting that vanilla projected gradient descent requires projection
onto the $f$-divergence ball, which can be computationally complex for large
$L$ or when selecting certain $f$ divergence functions. To circumvent this
issue, we explore an alternative approach by leveraging conditional gradient
methods (Frank-Wolfe [FW56]) that avoid the need for such projections. This
adoption of conditional gradient method is also motivated by the observation
that minimizing a linear objective over an $f$-divergence ball has a simpler
to analyze dual formulation [BTDHDW+13, ND16]. We first provide a brief
overview of conditional gradient methods.

We consider a convex differentiable objective function
$g:{\mathbb{R}}^{L}\to{\mathbb{R}}$ and a compact convex set
$\mathcal{P}\subset{\mathbb{R}}^{L}$. Given the optimization problem of
minimizing $g$ over $\mathcal{P}$, the following iterative rule is considered
by conditional gradient methods:

$p_{t+1}=\gamma_{t}q_{t}+(1-\gamma_{t})p_{t}\,,\quad
q_{t}=\arg\min_{q\in\mathcal{P}}q^{\sf T}\nabla g(p_{t})\,,$ (26)

where $\gamma_{t}$ is the step size and can beset as
$\gamma_{t}=\frac{1}{t+2}$.

Getting back to initial optimization problems (5), we have

$U^{\mathsf{finite}}_{\tau}(\mathbf{V}_{n,L})=\min_{p\in\mathcal{U}_{\tau}}\frac{1}{n}\sum_{\ell=1}^{L}\frac{(V_{\ell}-np_{\ell})^{2}}{p_{\ell}+1/L}\,,\quad
U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})=\min_{p\in\mathcal{U}_{\tau}}\frac{1}{n}\sum_{\ell=1}^{L}\frac{(V_{\ell}-np_{\ell})^{2}}{p_{\ell}}\,,$
(27)

where the constraint set $\mathcal{U}_{\tau}$ is given by

${\mathcal{U}}_{\tau}=\left\\{p\in{\mathbb{R}}^{L}\text{ s.t. }\quad p\geq
0,\quad\sum_{\ell=1}^{L}p_{\ell}=1,\quad\sum_{\ell=1}^{L}f(Lp_{\ell})\leq
L\tau\right\\}\,.$

We next focus on solving the linearization of objective functions (26) over
the $f$-divergence ball, which is needed in the update rule (26). This problem
has been studied before by [BTDHDW+13, Sha17] for general continuous
distributions. For the reader’s convenience, we state the result for the
specific case of the problem over discrete distributions.

######  Proposition 6.1.

For $x\in{\mathbb{R}}^{L}$ and $\tau\geq 0$, let

$\displaystyle(\lambda^{*},\eta^{*})\in\arg\min_{\lambda\geq
0,\eta\in{\mathbb{R}}}\left[\lambda\tau+\eta+\frac{\lambda}{L}\sum_{\ell=1}^{L}f^{*}\Big{(}\frac{-\eta-
x_{\ell}}{\lambda}\Big{)}\right]\,.$

If $q^{*}\in{\mathbb{R}}^{L}$ is such that

$\frac{-x_{\ell}-\eta^{*}}{L\lambda^{*}}\in\partial
f(q^{*}_{\ell})\,,\quad\forall\ell\in[L]\,,$

then $q^{*}$ is the minimizer of $q^{\sf T}x$ over ${\mathcal{U}}_{\tau}$.

It is easy to observe that the objective function in Proposition 6.1 is over a
half-space in a two-dimensional space with decoupled constraints, and can be
solved efficiently. Combining the result in Proposition 6.1 with (26)
completes the iterative procedure to solve optimization problems (27). Since
objective functions (27) are different for asymptotic and finite decision
rules, therefore we have distinct values for gradients ($x$ in $q^{\sf T}x$ as
per Proposition 6.1). Namely, we have

$\nabla
g^{\mathsf{asym}}(p)=n-\frac{\mathbf{V}_{n,L}(\ell)}{np_{\ell}^{2}}\,,\quad\nabla
g^{\mathsf{finite}}(p)=\frac{(np_{\ell}-\mathbf{V}_{n,L}(\ell))\left(np_{\ell}+\mathbf{V}_{n,L}(\ell)+\frac{2n}{L}\right)}{n(p_{\ell}+\frac{1}{L})}\,,$

where $g^{\mathsf{asym}}$ and $g^{\mathsf{finite}}$ are associate objective
functions in (27). Algorithms 4 and 5 summarizes the iterative procedures for
solving finite and asymptotic optimization problems (27), respectively.

1: (initialization): $p_{\ell}(0)\leftarrow 1/L\,,$ for all $\ell\in[L]$.

2: for iteration $h=1,\dotsc,H$ do

3:
$x_{\ell}(h)=\frac{(np_{\ell}(h)-\mathbf{V}_{\ell})\left(np_{\ell}(h)+\mathbf{V}_{\ell}+\frac{2n}{L}\right)}{n(p_{\ell}(h)+\frac{1}{L})}$

4: $(\lambda(h),\eta(h))=\arg\min\limits_{\lambda\geq
0,\eta\in{\mathbb{R}}}\left\\{\lambda\tau+\eta+\frac{\lambda}{L}\sum\limits_{\ell=1}^{L}f^{*}\Big{(}\frac{-\eta-
x_{\ell}(h)}{\lambda}\Big{)}\right\\}$

5: let $q_{\ell}(h)$ be such that
$\frac{-\eta(h)-x_{\ell}(h)}{L\lambda(h)}\in\partial f(q_{\ell}(h))\,,$ for
all $\ell\in[L]$

6: $p(h+1)=\gamma(h)q(h)+(1-\gamma(h))p(h)$

7: return
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\leftarrow\frac{1}{n}\sum\limits_{\ell=1}^{L}\frac{(\mathbf{V}_{\ell}-np_{\ell}(H))^{2}}{(p_{\ell}(H)+\frac{1}{L})^{2}}$

Algorithm 4 Iterative procedure to solve optimization problem of
$U^{\mathsf{finite}}_{\tau}(\mathbf{V}_{n,L})$

1: (initialization): $p_{\ell}(0)\leftarrow 1/L\,,$ for all $\ell\in[L]$.

2: for iteration $h=1,\dotsc,H$ do

3: $x_{\ell}(h)=n-\frac{\mathbf{V}_{\ell}^{2}}{np_{\ell}(h)^{2}}\,,$ for all
$\ell\in[L]$

4: $(\lambda(h),\eta(h))=\arg\min\limits_{\lambda\geq
0,\eta\in{\mathbb{R}}}\left\\{\lambda\tau+\eta+\frac{\lambda}{L}\sum\limits_{\ell=1}^{L}f^{*}\Big{(}\frac{-\eta-
x_{\ell}(h)}{\lambda}\Big{)}\right\\}$

5: let $q_{\ell}(h)$ be such that
$\frac{-\eta(h)-x_{\ell}(h)}{L\lambda(h)}\in\partial f(q_{\ell}(h))\,,$ for
all $\ell\in[L]$

6: $p(h+1)=\gamma(h)q(h)+(1-\gamma(h))p(h)$

7: return
$U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})\leftarrow\frac{1}{n}\sum\limits_{\ell=1}^{L}\frac{(\mathbf{V}_{\ell}-np_{\ell}(H))^{2}}{p_{\ell}(H)^{2}}$

Algorithm 5 Iterative procedure to solve optimization problem of
$U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})$

## References

  * [Aka98] Hirotogu Akaike.  Information theory and an extension of the maximum likelihood principle.  In Selected papers of hirotugu akaike, pages 199–213. Springer, 1998.
  * [BC15] Rina Foygel Barber and Emmanuel J Candès.  Controlling the false discovery rate via knockoffs.  The Annals of Statistics, 43(5):2055–2085, 2015.
  * [BCJW21] Stephen Bates, Emmanuel Candès, Lucas Janson, and Wenshuo Wang.  Metropolized knockoff sampling.  Journal of the American Statistical Association, 116(535):1413–1427, 2021.
  * [BCS20] Rina Foygel Barber, Emmanuel J Candès, and Richard J Samworth.  Robust inference with knockoffs.  The Annals of Statistics, 48(3):1409–1431, 2020.
  * [BHT21] Stephen Bates, Trevor Hastie, and Robert Tibshirani.  Cross-validation: what does it estimate and how well does it do it?  arXiv preprint arXiv:2104.00673, 2021.
  * [Bir01] Lucien Birgé.  An alternative point of view on lepski’s method.  Lecture Notes-Monograph Series, pages 113–133, 2001.
  * [BTDHDW+13] Aharon Ben-Tal, Dick Den Hertog, Anja De Waegenaere, Bertrand Melenberg, and Gijs Rennen.  Robust solutions of optimization problems affected by uncertain probabilities.  Management Science, 59(2):341–357, 2013.
  * [BW+19] Sivaraman Balakrishnan, Larry Wasserman, et al.  Hypothesis testing for densities and high-dimensional multinomials: Sharp local minimax rates.  Annals of Statistics, 47(4):1893–1927, 2019.
  * [BWBS20] Thomas B Berrett, Yi Wang, Rina Foygel Barber, and Richard J Samworth.  The conditional permutation test for independence while controlling for confounders.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(1):175–197, 2020.
  * [CFJL18] Emmanuel Candès, Yingying Fan, Lucas Janson, and Jinchi Lv.  Panning for gold:‘model-x’knockoffs for high dimensional controlled variable selection.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(3):551–577, 2018.
  * [DFH+15] Cynthia Dwork, Vitaly Feldman, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Aaron Leon Roth.  Preserving statistical validity in adaptive data analysis.  In Proceedings of the forty-seventh annual ACM symposium on Theory of computing, pages 117–126, 2015.
  * [DFH+17] Cynthia Dwork, Vitaly Feldman, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Aaron Roth.  Guilt-free data reuse.  Communications of the ACM, 60(4):86–93, 2017.
  * [Efr86] Bradley Efron.  How biased is the apparent error rate of a prediction rule?  Journal of the American statistical Association, 81(394):461–470, 1986.
  * [Efr04] Bradley Efron.  The estimation of prediction error: covariance penalties and cross-validation.  Journal of the American Statistical Association, 99(467):619–632, 2004.
  * [Far96] CP Farrington.  On assessing goodness of fit of generalized linear models to sparse data.  Journal of the Royal Statistical Society: Series B (Methodological), 58(2):349–360, 1996.
  * [FW56] Marguerite Frank and Philip Wolfe.  An algorithm for quadratic programming.  Naval research logistics quarterly, 3(1-2):95–110, 1956.
  * [GPAM+14] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [GPSW17] Chuan Guo, Geoff Pleiss, Yu Sun, and Kilian Q Weinberger.  On calibration of modern neural networks.  In International conference on machine learning, pages 1321–1330. PMLR, 2017.
  * [HH02] David W Hosmer and Nils Lid Hjort.  Goodness-of-fit processes for logistic regression: simulation results.  Statistics in medicine, 21(18):2723–2738, 2002.
  * [HL80] David W Hosmer and Stanley Lemesbow.  Goodness of fit tests for the multiple logistic regression model.  Communications in statistics-Theory and Methods, 9(10):1043–1069, 1980.
  * [HTFF09] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman.  The elements of statistical learning: data mining, inference, and prediction, volume 2.  Springer, 2009.
  * [IIS03] Yuri Ingster, Jurij I Ingster, and IA Suslina.  Nonparametric goodness-of-fit testing under Gaussian models, volume 169.  Springer Science & Business Media, 2003.
  * [IS09] Yu I Ingster and Theofanis Sapatinas.  Minimax goodness-of-fit testing in multivariate nonparametric regression.  Mathematical Methods of Statistics, 18:241–269, 2009.
  * [JM18] Adel Javanmard and Andrea Montanari.  Online rules for control of false discovery rate and false discovery exceedance.  The Annals of statistics, 46(2):526–554, 2018.
  * [JM21] Adel Javanmard and Mohammad Mehrabi.  Pearson chi-squared conditional randomization test.  arXiv preprint arXiv:2111.00027, 2021.
  * [JSBS20] Jana Janková, Rajen D Shah, Peter Bühlmann, and Richard J Samworth.  Goodness-of-fit testing in high dimensional generalized linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 82(3):773–795, 2020.
  * [KLM19] Ananya Kumar, Percy S Liang, and Tengyu Ma.  Verified uncertainty calibration.  Advances in Neural Information Processing Systems, 32, 2019.
  * [Lan59] E Lann.  Testing statistical hypotheses.  Wiley, New York, 1959.
  * [LCVH91] S Le Cessie and JC Van Houwelingen.  A goodness-of-fit test for binary regression models, based on smoothing methods.  Biometrics, pages 1267–1282, 1991.
  * [LHHD22] Donghwan Lee, Xinmeng Huang, Hamed Hassani, and Edgar Dobriban.  T-cal: An optimal test for the calibration of predictive models.  arXiv preprint arXiv:2203.01850, 2022.
  * [LS99] Oleg V Lepski and Vladimir G Spokoiny.  Minimax nonparametric hypothesis testing: the case of an inhomogeneous alternative.  Bernoulli, pages 333–358, 1999.
  * [LWY02] DY Lin, LJ Wei, and Z Ying.  Model-checking techniques based on cumulative residuals.  Biometrics, 58(1):1–12, 2002.
  * [Mal00] Colin L Mallows.  Some comments on cp.  Technometrics, 42(1):87–94, 2000.
  * [MS18] Rajarshi Mukherjee and Subhabrata Sen.  Optimal adaptive inference in random design binary regression.  2018\.
  * [NCH15] Mahdi Pakdaman Naeini, Gregory Cooper, and Milos Hauskrecht.  Obtaining well calibrated probabilities using bayesian binning.  In Twenty-Ninth AAAI Conference on Artificial Intelligence, 2015\.
  * [NCT16] Sebastian Nowozin, Botond Cseke, and Ryota Tomioka.  f-gan: Training generative neural samplers using variational divergence minimization.  Advances in neural information processing systems, 29, 2016.
  * [ND16] Hongseok Namkoong and John C Duchi.  Stochastic gradient methods for distributionally robust optimization with f-divergences.  Advances in neural information processing systems, 29, 2016.
  * [nix19] Measuring Calibration in Deep Learning., volume 2, 2019.
  * [NWJ10] XuanLong Nguyen, Martin J Wainwright, and Michael I Jordan.  Estimating divergence functionals and the likelihood ratio by convex risk minimization.  IEEE Transactions on Information Theory, 56(11):5847–5861, 2010\.
  * [OR92] Gerhard Osius and Dieter Rojek.  Normal goodness-of-fit tests for multinomial models with large degrees of freedom.  Journal of the American Statistical Association, 87(420):1145–1152, 1992.
  * [RAN14] Saharon Rosset, Ehud Aharoni, and Hani Neuvirth.  Novel statistical tools for management of public databases facilitate community-wide replicability and control of false discovery.  Genetic epidemiology, 38(5):477–481, 2014.
  * [RT19] Saharon Rosset and Ryan J Tibshirani.  From fixed-x to random-x regression: Bias-variance decompositions, covariance penalties, and prediction error estimation.  Journal of the American Statistical Association, 2019.
  * [SB18] Rajen D Shah and Peter Bühlmann.  Goodness-of-fit tests for high dimensional linear models.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(1):113–135, 2018.
  * [Sha17] Alexander Shapiro.  Distributionally robust stochastic programming.  SIAM Journal on Optimization, 27(4):2258–2275, 2017.
  * [SRA+08] Andrea Saltelli, Marco Ratto, Terry Andres, Francesca Campolongo, Jessica Cariboni, Debora Gatelli, Michaela Saisana, and Stefano Tarantola.  Global sensitivity analysis: the primer.  John Wiley & Sons, 2008.
  * [SW91] John Q Su and LJ Wei.  A lack-of-fit test for the mean function in a generalized linear model.  Journal of the American Statistical Association, 86(414):420–426, 1991.
  * [Tsi80] Anastasios A Tsiatis.  A note on a goodness-of-fit test for the logistic regression model.  Biometrika, 67(1):250–251, 1980.
  * [TVZ+22] Wesley Tansey, Victor Veitch, Haoran Zhang, Raul Rabadan, and David M Blei.  The holdout randomization test for feature selection in black box models.  Journal of Computational and Graphical Statistics, 31(1):151–162, 2022.
  * [VWA+19] Juozas Vaicenavicius, David Widmann, Carl Andersson, Fredrik Lindsten, Jacob Roll, and Thomas Schön.  Evaluating model calibration in classification.  In The 22nd International Conference on Artificial Intelligence and Statistics, pages 3459–3467. PMLR, 2019.
  * [Wag20] Stefan Wager.  Cross-validation, risk estimation, and model selection: Comment on a paper by rosset and tibshirani.  Journal of the American Statistical Association, 115(529):157–160, 2020.
  * [WLZ19] David Widmann, Fredrik Lindsten, and Dave Zachariah.  Calibration tests in multi-class classification: A unifying framework.  Advances in Neural Information Processing Systems, 32, 2019.
  * [ZDY21] Jiawei Zhang, Jie Ding, and Yuhong Yang.  Is a classification procedure good enough?—a goodness-of-fit assessment tool for classification learning.  Journal of the American Statistical Association, pages 1–11, 2021\.

Supplementary material for “GRASP: A Goodness-of-Fit Test for Classification
Learning”

Adel Javanmard and Mohammad Mehrabi

##  Appendix A Proofs of theorems and lemmas

###  A.1 Proof of Lemma 2.2

Using (2) for $f(t)=\frac{1}{2}|t-1|$ we get

$\displaystyle{\mathbb{E}}\left[D_{\mathsf{TV}}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]$
$\displaystyle={\mathbb{E}}\left[\frac{\widehat{\eta}(X)}{2}\Big{|}\frac{\eta(X)}{\widehat{\eta}(X)}-1\Big{|}+\frac{1-\widehat{\eta}(X)}{2}\Big{|}\frac{1-\eta(X)}{1-\widehat{\eta}(X)}-1\Big{|}\right]$
$\displaystyle={\mathbb{E}}\left[|\widehat{\eta}(X)-\eta(X)|\right]\,.$

By a similar argument for $f(t)=t\log(t)$ we get

$\displaystyle{\mathbb{E}}\left[D_{\mathsf{KL}}(\mathsf{Bern}(\eta(X)\|\mathsf{Bern}(\widehat{\eta}(X)))\right]$
$\displaystyle={\mathbb{E}}\left[\widehat{\eta}(X)\frac{\eta(X)}{\widehat{\eta}(X)}\log\frac{\eta(X)}{\widehat{\eta}(X)}+(1-\widehat{\eta}(X))\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\log\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right]$
$\displaystyle={\mathbb{E}}\left[\eta(X)\log\frac{\eta(X)}{\widehat{\eta}(X)}+(1-\eta(X))\log\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right]\,.$

We then expand the expression inside logarithm to get

$\displaystyle{\mathbb{E}}\left[D_{\mathsf{KL}}(\mathsf{Bern}(\eta(X)\|\mathsf{Bern}(\widehat{\eta}(X)))\right]$
$\displaystyle=-{\mathbb{E}}\left[\eta(X)\log{\widehat{\eta}(X)}+(1-\eta(X))\log{(1-\widehat{\eta}(X))}\right]$
$\displaystyle\;+{\mathbb{E}}\left[\eta(X)\log{\eta(X)}+(1-\eta(X))\log{(1-\eta(X))}\right]$
$\displaystyle=\mathsf{CE}(\widehat{\eta})-\mathsf{CE}(\eta)\,.$

Next, we prove the result for the Hellinger distance. By using (2) for
$f(t)=(\sqrt{t}-1)^{2}$ we get

$\displaystyle{\mathbb{E}}\left[D_{\mathsf{H}}(\mathsf{Bern}(\eta(X)||\mathsf{Bern}(\widehat{\eta}(X)))\right]$
$\displaystyle={\mathbb{E}}\left[\widehat{\eta}(X)\left(\sqrt{\frac{\eta(X)}{\widehat{\eta}(X)}}-1\right)^{2}+(1-\widehat{\eta}(X))\left(\sqrt{\frac{1-\eta(X)}{1-\widehat{\eta}(X)}}-1\right)^{2}\right]$
$\displaystyle={\mathbb{E}}\left[\Big{(}\sqrt{\eta(X)}-\sqrt{\widehat{\eta}(X)}\Big{)}^{2}+\Big{(}\sqrt{1-\eta(X)}-\sqrt{1-\widehat{\eta}(X)}\Big{)}^{2}\right]\,.$

###  A.2 Proof of Proposition 3.1

In oder to compute
${\mathbb{E}}[D_{f}(\mathcal{L}(W|X)\|\mathsf{Unif}([0,1]))]$, we first need
to characterize the conditional distribution $\mathcal{L}(W|X)$. From the
mechanism to construct $w$, it is easy to obtain

$\displaystyle p_{W|X}(w|x)$
$\displaystyle={\mathbb{P}}(y=+1|x)p_{W|Y,X}(w|y=+1,x)+{\mathbb{P}}(y=0|x)p_{W|Y,X}(w|y=0,x)$
$\displaystyle=\frac{\eta(x)}{\widehat{\eta}(x)}{\mathbb{I}}(w\leq\widehat{\eta}(x))+\frac{1-\eta(x)}{1-\widehat{\eta}(x)}{\mathbb{I}}(w\geq\widehat{\eta}(x))\,.$
(28)

Using (28) in $f$-divergence definition 2.1 yields

$\displaystyle{\mathbb{E}}\left[D_{f}(\mathcal{L}(W|X)\|\mathsf{Unif}([0,1]))\right]$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f(p_{W|X}(w|X)){\rm d}w\right]$
$\displaystyle={\mathbb{E}}\left[\int\limits_{0}^{\widehat{\eta}(X)}f\left(\frac{\eta(X)}{\widehat{\eta}(X)}\right){\rm
d}t+\int\limits_{\widehat{\eta}(X)}^{1}f\left(\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right){\rm
d}t\right]\,.$
$\displaystyle={\mathbb{E}}\left[\widehat{\eta}(X)f\left(\frac{\eta(X)}{\widehat{\eta}(X)}\right)+(1-\widehat{\eta}(X))f\left(\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right)\right]$
$\displaystyle={\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(X)\|\mathsf{Bern}(\widehat{\eta}(X)))\right]\,.$

###  A.3 Proof of Proposition 3.2

From the optimization problems used in the definition of statistics
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ and
$U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$ in (5), it is easy to get
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\leq
U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})$. In fact, this is due to the extra
$1/L$ term in the denominator of the test statistic
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$. We also need to compare the
rejection thresholds. For this end, we first provide the following upper bound
on the quantiles of a chi-squared distribution with $L-1$ degrees of freedom.
By using ([Bir01], Lemma 8.1),

$\chi^{2}_{L-1}(1-\alpha)\leq
L-1+2\sqrt{(L-1)\log\frac{1}{\alpha}}+2\log\frac{1}{\alpha}\,.$

Also it is easy to get that for $L\geq 60$ the following holds (we prove this
later):

$L-1+2\sqrt{(L-1)\log\frac{1}{\alpha}}+2\log\frac{1}{\alpha}\leq
L+\sqrt{\frac{2L}{\alpha}}\,.$ (29)

Putting all together and by recalling the definition of
$\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$ we get

$\displaystyle\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}$
$\displaystyle={\mathbb{I}}(U_{\tau}^{\mathsf{asym}}(\mathbf{V}_{n,L})\geq\chi^{2}_{L-1}(1-\alpha))$
$\displaystyle\geq{\mathbb{I}}(U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\geq\chi^{2}_{L-1}(1-\alpha))$
$\displaystyle\geq{\mathbb{I}}\left(U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\geq
L+\sqrt{\frac{2L}{\alpha}}\right)=\Phi_{n,L,\alpha,\tau}^{\mathsf{asym}}\,.$

We only need to prove (29). By upper bounding $L-1$ with $L$ and multiplying
both sides of (29) by $0.5\sqrt{\alpha/L}$ the relation reduces to showing

$\sqrt{\alpha\log\frac{1}{\alpha}}+\sqrt{\frac{\alpha}{L}}\log\frac{1}{\alpha}\leq\frac{1}{\sqrt{2}}\,.$

By simple algebraic calculation, it is easy to see that
$\alpha\log{1/\alpha}\leq 1/e$ and $\sqrt{\alpha}\log{1/\alpha}\leq 2/e$ and
so

$\sqrt{\alpha\log\frac{1}{\alpha}}+\sqrt{\frac{\alpha}{L}}\log\frac{1}{\alpha}\leq\frac{1}{\sqrt{e}}+\frac{2}{e\sqrt{L}}\,.$

Next note that for $L\geq 55$ we have
$\frac{1}{\sqrt{e}}+\frac{2}{e\sqrt{L}}\leq\frac{1}{\sqrt{2}}$, which
completes the proof.

###  A.4 Proof of Proposition 3.4

It is easy to observe that the joint distribution of $(V_{1},\dots,V_{L})$ is
a multinomial distribution with probabilities $(p_{1},\dots,p_{L})$ where
$p_{\ell}$ denotes the probability of sample $(X,W)$ admitting label $\ell$
for $\ell\in[L]$. As per Algorithm 1, let
$\widetilde{W}_{1},\dots,\widetilde{W}_{M}$ be iid random variables with
uniform distribution $\mathsf{Unif}[0,1]$. To lighten the notation, we use the
shorthands $T=T(X,W)$, and $T_{j}=T(X,\widetilde{W}_{j})$ for $j\in[M]$. In
this case, the rank value is given by

$R={1+\sum_{j=1}^{M}{\mathbb{I}}\big{(}T\geq T_{j}\big{)}}\,.$

We get

$\displaystyle{\mathbb{P}}((X,W)\text{ admits label }\ell)$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}{\mathbb{E}}[{\mathbb{P}}(R=j+1|X)]$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}{\mathbb{E}}\left[{\mathbb{P}}(T\text{
is exactly larger than }j\text{ of }T_{j\in[M]}|X)\right]\,.$

We know that conditioned on $X$, random variables $T$, $T_{j\in[M]}$ are
independent. This gives us

$\displaystyle{\mathbb{P}}((X,W)\text{ admits label }\ell)$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}{\mathbb{E}}\left[\int{\mathbb{P}}(\text{
exactly }j\text{ of }T_{j\in[M]}\text{ are smaller than }t|X){\rm
d}P_{T|X}(t|X)\right]$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}{\mathbb{E}}\left[\int\binom{M}{j}F_{c}(t;X)^{j}(1-F_{c}(t;X))^{M-j}\frac{\partial
F_{o}(t;X)}{\partial t}{\rm d}t\right]\,.$ (30)

The last relation follows the iid property of random variables $T_{j\in[M]}$
conditioned on $X$ along with the definitions of $F_{c}(t;x),F_{o}(t;x)$ as
per Assumption 3.3. By adopting the change of variable $u=F_{c}(t;X)$ in the
integral (30) we arrive at

${\mathbb{P}}((X,W)\text{ admits label
}\ell)=\sum_{j=K(\ell-1)}^{K\ell-1}{\mathbb{E}}\left[\int_{0}^{1}\binom{M}{j}u^{j}(1-u)^{M-j}\frac{\partial
F_{o}(F_{c}^{-1}(u;X);X)}{\partial u}{\rm d}u\right]\,.$

We next recall the function $\psi(u;X)=\frac{\partial
F_{o}(F_{c}^{-1}(u;X);X)}{\partial u}$ from Assumption 3.3, which allows us to
write the above probability as

$p_{\ell}=\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}{\mathbb{E}}\left[\int_{0}^{1}u^{j}(1-u)^{M-j}\psi(u;X){\rm
d}u\right]\,.\quad\forall\ell\in[L]\,.$ (31)

For $\ell\in[L]$, we define the polynomial $\beta_{\ell}:[0,1]\to{\mathbb{R}}$
as

$\beta_{\ell}(u)=\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}u^{j}(1-u)^{M-j}\,.$

By using the above in (31), we obtain

$p_{\ell}={\mathbb{E}}\left[\int_{0}^{1}\beta_{\ell}(u)\psi(u;X){\rm
d}u\right]\,.\quad\forall\ell\in[L]\,,$ (32)

and hence

$\sum_{\ell=1}^{L}\frac{f(p_{\ell}L)}{L}=\frac{1}{L}\sum_{\ell=1}^{L}f\left(L{\mathbb{E}}\Big{[}\int_{0}^{1}\beta_{\ell}(u)\psi(u;X){\rm
d}u\Big{]}\right)\,.$ (33)

We next upper bound the right-hand side in the above equation. Consider the
following probability density function over interval $[0,1]$:

$\varphi_{\ell}(u)=\frac{\beta_{\ell}(u)}{\int_{0}^{1}\beta_{\ell}(u){\rm
d}u}\,,\quad\forall u\in[0,1]\,.$ (34)

By an application of Jensen’s inequality and the fact that $f$ is a convex
function, we have

$f\left({\mathbb{E}}_{X}\Big{[}{\mathbb{E}}_{u\sim\varphi_{\ell}}[\psi(u;X)]\Big{]}\right)\leq{\mathbb{E}}_{X}\left[{\mathbb{E}}_{u\sim\varphi_{\ell}}\Big{[}f(\psi(u;X))\Big{]}\right]\,.$
(35)

Then, by rewriting the expectation of $u$ in terms of density function
$\varphi_{\ell}$ we get

$f\left({\mathbb{E}}_{X}\left[\int_{0}^{1}\varphi_{\ell}(u)\psi(u;X){\rm
d}u\right]\right)\leq{\mathbb{E}}_{X}\left[\int_{0}^{1}\varphi_{\ell}(u)f(\psi(u;X)){\rm
d}u\right]\,.$ (36)

Plugging (34) into (36) yields

$f\left({\mathbb{E}}\left[\frac{\int_{0}^{1}\beta_{\ell}(u)\psi(u;X){\rm
d}u}{\int_{0}^{1}\beta_{\ell}(u){\rm
d}u}\right]\right)\leq{\mathbb{E}}\left[\frac{\int_{0}^{1}\beta_{\ell}(u)f(\psi(u;X)){\rm
d}u}{\int_{0}^{1}\beta_{\ell}(u){\rm d}u}\right]\,.$ (37)

In addition, from the definition of $\beta_{\ell}(u)$ we have

$\displaystyle\int_{0}^{1}\beta_{\ell}(u){\rm d}u$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}\int_{0}^{1}u^{j}(1-u)^{M-j}{\rm
d}u$ $\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}B(j+1,M-j+1)$
$\displaystyle=\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}\frac{j!(M-j)!}{(M+1)!}=\frac{1}{L}\,,$
(38)

where $B(a,b)$ is the beta function. Using (38) in (37) brings us to

$f\left(L{\mathbb{E}}\Big{[}\int_{0}^{1}\beta_{\ell}(u)\psi(u;X){\rm
d}u\Big{]}\right)\leq
L{\mathbb{E}}\left[\int_{0}^{1}\beta_{\ell}(u)f(\psi(u;X)){\rm d}u\right]$
(39)

In the next step, combining (39) with (33) gives us

$\displaystyle\sum_{\ell=1}^{L}\frac{f(p_{\ell}L)}{L}$
$\displaystyle\leq\sum_{\ell=1}^{L}{\mathbb{E}}\left[\int_{0}^{1}\beta_{\ell}(u)f(\psi(u;X)){\rm
d}u\right]$ $\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f(\psi(u;X)){\rm
d}u\right]\,,$ (40)

where in the last relation we used the identity
$\sum_{\ell=1}^{L}\beta_{\ell}(u)=1$, because

$\displaystyle\sum_{\ell=1}^{L}\beta_{\ell}(u)$
$\displaystyle=\sum_{\ell=1}^{L}\sum_{j=K(\ell-1)}^{K\ell-1}\binom{M}{j}u^{j}(1-u)^{M-j}$
$\displaystyle=\sum_{j=0}^{M}\binom{M}{j}u^{j}(1-u)^{M-j}=(u+1-u)^{M}=1\,.$

We next use definition of $\psi(u;x)$ in (40) to obtain

$\displaystyle{\mathbb{E}}\left[\int_{0}^{1}f(\psi(u;X)){\rm d}u\right]$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f\Big{(}\frac{\partial
F_{o}(F_{c}^{-1}(u;X);X)}{\partial u}\Big{)}{\rm d}u\right]\,.$ (41)

In the next step, by using density functions $f_{o}(t;x),f_{c}(t;x)$ as per
Assumption 3.3 in (41) we arrive at

$\displaystyle{\mathbb{E}}\left[\int_{0}^{1}f(\psi(u;X)){\rm d}u\right]$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f\Big{(}\frac{f_{o}(F_{c}^{-1}(u;X);X)}{f_{c}(F_{c}^{-1}(u;X);X)}\Big{)}{\rm
d}u\right]$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f\Big{(}\frac{f_{o}(t;X)}{f_{c}(t;X)}\Big{)}f_{c}(t;X){\rm
d}t\right]$
$\displaystyle={\mathbb{E}}\left[D_{f}\Big{(}\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\Big{)}\right]\,,$
(42)

where the last relation comes from Assumption 3.3 and the definition of
$f$-divergence function. In the next step, by combining (42) and (40) we get

$\displaystyle\sum_{\ell=1}^{L}\frac{f(p_{\ell}L)}{L}\leq{\mathbb{E}}\Big{[}D_{f}\left(\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\right)\Big{]}\,.$
(43)

Further, by an application of the data processing inequality we have

$D_{f}\left(\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\right)\leq
D_{f}\left(\mathcal{L}(W|X)\|\mathcal{L}(\widetilde{W}|X)\right)\,.$ (44)

Finally, combining (44) and (43) gives us

$\displaystyle\sum_{\ell=1}^{L}\frac{f(p_{\ell}L)}{L}$
$\displaystyle\leq{\mathbb{E}}\Big{[}D_{f}\left(\mathcal{L}(T(X,W)|X)\|\mathcal{L}(T(X,\widetilde{W})|X)\right)\Big{]}$
$\displaystyle\leq{\mathbb{E}}\left[D_{f}(\mathcal{L}(W|X)\|\mathcal{L}(\widetilde{W}|X))\right]\,.$

###  A.5 Proof of Theorem 3.6

We know that $\mathbf{V}_{n,L}$ has a multinomial distribution. Let $q_{\ell}$
denote the probability of occurrence for category $\ell\in[L]$. Invoking the
result of [BW+19, Theorem 2] for multinomial hypothesis testing with truncated
chi-squared test statistics, we get

$\displaystyle\alpha$
$\displaystyle\geq{\mathbb{P}}\left(\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}-V_{\ell}}{\max\\{q_{\ell},\frac{1}{L}\\}}\geq
n\sqrt{\frac{2}{\alpha}\sum\limits_{\ell=1}^{L}\left(\frac{q_{\ell}}{\max\\{q_{\ell},1/L\\}}\right)^{2}}\right)\,.$
(45)

In the next step, we use $\max(q_{\ell},1/L)\geq 1/L$ along with
$\sum_{\ell=1}^{L}q_{\ell}=1$ to get

$\displaystyle\alpha$
$\displaystyle\geq{\mathbb{P}}\left(\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}-V_{\ell}}{\max\\{q_{\ell},\frac{1}{L}\\}}\geq
n\sqrt{\frac{2}{\alpha}L}\right)$
$\displaystyle={\mathbb{P}}\left(\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{\max\\{q_{\ell},\frac{1}{L}\\}}\geq\sum\limits_{\ell=1}^{L}\frac{V_{\ell}}{\max\\{q_{\ell},\frac{1}{L}\\}}+n\sqrt{\frac{2}{\alpha}L}\right)$
$\displaystyle\geq{\mathbb{P}}\left(\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{\max\\{q_{\ell},\frac{1}{L}\\}}\geq
L\sum\limits_{\ell=1}^{L}{V_{\ell}}+n\sqrt{\frac{2}{\alpha}L}\right)$
$\displaystyle={\mathbb{P}}\left(\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{\max\\{q_{\ell},\frac{1}{L}\\}}\geq
nL+n\sqrt{\frac{2}{\alpha}L}\right)$
$\displaystyle\geq{\mathbb{P}}\left(\frac{1}{n}\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{q_{\ell}+\frac{1}{L}}\geq
L+\sqrt{\frac{2}{\alpha}L}\right)\,.$ (46)

By deploying Theorem 3.5, under the null hypothesis (3) we have

$\frac{1}{L}\sum_{\ell=1}^{L}f(Lq_{\ell})\leq{\mathbb{E}}\left[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))\right]\leq\tau\,.$
(47)

This implies that $(q_{1},\dots,q_{L})$ is a feasible point for the
optimization problem in the definition of
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ in (5), and so

$\frac{1}{n}\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{q_{\ell}+\frac{1}{L}}\geq
U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\,.$ (48)

Plugging (48) into (46) yields

${\mathbb{P}}\left(U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})\geq
L+\sqrt{\frac{2}{\alpha}L}\right)\leq\alpha\,.$

This completes the proof by the definition of decision rule
$\Phi^{\mathsf{finite}}_{n,L,\tau}$.

We next proceed to the proof of the asymptotic result. Following a similar
argument shows that $(q_{1},\dotsc,q_{\ell})$ is also a feasible solution for
the optimization problem in the definition of
$U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})$, which implies that

$\sum_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{q_{\ell}}\geq
U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})\,.$ (49)

In the next step, we use the following asymptotic result on the Pearson’s chi-
squared test statistic (see e.g., [Lan59] Theorem 14.3.1):

$\lim_{n\to\infty}\sum\limits_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{nq_{\ell}}\overset{(d)}{\to}\chi^{2}_{L-1}\,,$
(50)

where $\chi^{2}_{L-1}$ denotes the chi-squared distribution with $L-1$ degrees
of freedom. As a direct result of (50) we have

$\lim_{n\to\infty}\sup{\mathbb{P}}\left(\sum\limits_{\ell=1}^{L}\frac{(V_{\ell}-nq_{\ell})^{2}}{q_{\ell}}\geq\chi^{2}_{L-1}(1-\alpha)\right)=\alpha\,.$

Combining (49) with the above relation yields

$\lim_{n\to\infty}\sup{\mathbb{P}}\left(U^{\mathsf{asym}}_{\tau}(\mathbf{V}_{n,L})\geq\chi^{2}_{L-1}(1-\alpha)\right)\leq\alpha\,.$

The claim then follows simply from the definition of decision rule
$\Phi^{\mathsf{asym}}_{n,L,\tau}$.

###  A.6 Proof of Proposition 4.1

Let ${\widetilde{X}}\sim\mathcal{P}_{X}$ and
$\widetilde{W}\sim\mathsf{Unif}([0,1])$, then we have

$\displaystyle
D_{f}\left(\mathcal{L}(X,W)\|\mathcal{P}_{X}\times\mathsf{Unif}([0,1])\right)$
$\displaystyle=D_{f}\left(\mathcal{L}(X,W)\|\mathcal{L}({\widetilde{X}},\widetilde{W})\right)$
$\displaystyle=\int
p_{{\widetilde{X}},\widetilde{W}}(x,w)f\left(\frac{p_{X,W}(x,w)}{p_{{\widetilde{X}},\widetilde{W}}(x,w)}\right){\rm
d}x{\rm d}w$ $\displaystyle=\int
p_{{\widetilde{X}}}(x)p_{\widetilde{W}|{\widetilde{X}}}(w|x)f\left(\frac{p_{X}(x)p_{W|X}(x)}{p_{{\widetilde{X}}}(x)p_{\widetilde{W}|{\widetilde{X}}}(w|x)}\right){\rm
d}x{\rm d}w\,.$

Since $p_{{\widetilde{X}}}=p_{X}$ and
$p_{\widetilde{W}|{\widetilde{X}}}(w|x)=1$(uniform distribution) we get

$\displaystyle
D_{f}\left(\mathcal{L}(X,W)\|\mathcal{L}({\widetilde{X}},\widetilde{W})\right)$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{1}f\left({p_{W|X}(w|X)}\right){\rm
d}w\right]\,.$ (51)

Also from the construction of $w$, it is easy to see that

$\displaystyle p_{W|X}(w|x)$
$\displaystyle={\mathbb{P}}(y=+1|x)p_{W|Y,X}(w|y=+1,x)+{\mathbb{P}}(y=0|x)p_{W|Y,X}(w|y=0,x)$
$\displaystyle=\frac{\eta(x)}{\widehat{\eta}(x)}{\mathbb{I}}(w\leq\widehat{\eta}(x))+\frac{1-\eta(x)}{1-\widehat{\eta}(x)}{\mathbb{I}}(w\geq\widehat{\eta}(x))\,.$
(52)

Plugging (A.6) into (51) brings us to

$\displaystyle
D_{f}\left(\mathcal{L}(X,W)\|\mathcal{L}({\widetilde{X}},\widetilde{W})\right)$
$\displaystyle={\mathbb{E}}\left[\int_{0}^{\widehat{\eta}(X)}f\left(\frac{\eta(X)}{\widehat{\eta}(X)}\right){\rm
d}w+\int_{\widehat{\eta}(X)}^{1}f\left(\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right){\rm
d}w\right]$
$\displaystyle={\mathbb{E}}\left[\widehat{\eta}(X)f\left(\frac{\eta(X)}{\widehat{\eta}(X)}\right)+(1-\widehat{\eta}(X))f\left(\frac{1-\eta(X)}{1-\widehat{\eta}(X)}\right)\right]$
$\displaystyle={\mathbb{E}}\left[D_{f}\left(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X))\right)\right]\,,$

which completes the proof.

###  A.7 Proof of Proposition 3.7

We adopt the shorthand $\tau_{1}$ for $\tau^{\mathsf{finite}}_{n,L,\alpha}$
and $\tau_{0}$ for
${\mathbb{E}}[D_{f}(\mathsf{Bern}(\eta(X))\|\mathsf{Bern}(\widehat{\eta}(X)))]$.
From the definitions, $\tau_{1}$ is a function of $\mathbf{V}_{n,L}$ and is
random. We only prove the result for the finite test statistics. The claim for
the asymptotic statistics follows from the same argument.

First, from the optimization problem used for
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ in (5) it is easy to observe
that for fixed values of $\mathbf{V}_{n,L}$ statistics
$U_{\tau}^{\mathsf{finite}}(\mathbf{V}_{n,L})$ is non-increasing in $\tau$.
Therefore,

$\displaystyle{\mathbb{P}}({\mathbb{E}}[D_{f}(\mathsf{Bern}(\eta(x))\|\mathsf{Bern}(\widehat{\eta}(x)))]\geq\tau_{1})$
$\displaystyle={\mathbb{P}}(\tau_{0}\geq\tau_{1})$
$\displaystyle=1-{\mathbb{P}}(\tau_{1}\geq\tau_{0})$ $\displaystyle\geq
1-{\mathbb{P}}(U^{\mathsf{finite}}_{\tau_{0}}(\mathbf{V}_{n,L})\geq