[datatype=bibtex, overwrite] [fieldset=url, null] [fieldset=urldate, null]
[fieldset=doi, null] [fieldset=issn, null] [fieldset=isbn, null]
[fieldset=editor, null] [fieldset=publisher, null] [fieldset=location, null]

# Composite Goodness-of-fit Tests with Kernels

Oscar Key
University College London
<EMAIL_ADDRESS>
&Tamara Fernandez
Adolfo Ibañez University
<EMAIL_ADDRESS>
Arthur Gretton
University College London
<EMAIL_ADDRESS>
&François-Xavier Briol
University College London
<EMAIL_ADDRESS>
###### Abstract

Model misspecification can create significant challenges for the
implementation of probabilistic models, and this has led to development of a
range of inference methods which directly account for the issue. However,
these method tend to lose efficiency and should only be used when the model is
really misspecified. Unfortunately, there is a lack of generally applicable
methods to test whether this is the case or not. One set of tools which can
help are goodness-of-fit tests, which can test whether a dataset has been
generated from a _fixed_ distribution. Kernel-based tests have been developed
to for this problem, and these are popular due to their flexibility, strong
theoretical guarantees and ease of implementation in a wide range of
scenarios. In this paper, we extend this line of work to the more challenging
_composite goodness-of-fit_ problem, where we are instead interested in
whether the data comes from any distribution in some parametric family.

##  1 Introduction

One approach to answering the question of whether a model is misspecified is
goodness-of-fit testing [13, Chapter 14]. Given a _fixed_ distribution
$\operatorname{\mathbb{P}}$ and observations $x_{1},\ldots,x_{n}$ generated
from an unknown distribution $\operatorname{\mathbb{Q}}$, goodness-of-fit
tests compare the null hypothesis
$H_{0}:\operatorname{\mathbb{P}}=\operatorname{\mathbb{Q}}$ against the
alternative hypothesis
$H_{1}:\operatorname{\mathbb{P}}\neq\operatorname{\mathbb{Q}}$. Tests using
the kernel Stein discrepancy (KSD) as test statistic are popular for this task
because they can be applied to a wide-range of data types and can accommodate
models with unnormalised densities [15, 4].

Although these tests are very useful in practice, we will often be interested
in answering the more complex question of whether our data was generated by
_any element of some parametric family_ of distributions
$\\{\operatorname{\mathbb{P}}_{\theta}\\}_{\theta\in\Theta}$ with parameter
space $\Theta$. Specifically, the null hypothesis (corresponding to a well-
specified model) is $H^{C}_{0}:\exists\theta^{*}\in\Theta$ such that
$\operatorname{\mathbb{P}}_{\theta^{*}}=\operatorname{\mathbb{Q}}$, and the
alternative (corresponding to a misspecified model) is
$H^{C}_{1}:\operatorname{\mathbb{Q}}\not\in\\{\operatorname{\mathbb{P}}_{\theta}\\}_{\theta\in\Theta}$.
This type of test is known as a _composite goodness-of-fit test_ , and can be
much more challenging to construct since $\theta^{*}$ is usually unknown.

Our paper fills an important gap in the literature by proposing the first set
of _kernel-based composite hypothesis tests_ applicable to a wide range of
parametric models. These are in contrast to previously introduced composite
tests which are limited to very specific parametric families [12, 5]. To
devise these new tests, we make use of recently developed minimum distance
estimators based on the KSD [1]. A key challenge is that the dataset is used
twice, both to estimate the parameter and the test statistic, and this is done
without splitting it into estimation and test sets. To achieve the correct
level, the test must take account of this dependence, both via a more in-depth
theoretical analysis and by using a suitable method for approximating the
threshold. In this initial version of the work, we only consider the second
aspect by using the parametric bootstrap [23] to achieve the correct level
without data-splitting.

##  2 Background: Kernel Stein Discrepancies for Testing and Estimation

We will now briefly review the use of KSD for testing and estimation. Denote
by $\mathcal{X}$ the data space and $\mathcal{P}(\mathcal{X})$ the set of all
Borel distributions on $\mathcal{X}$. For simplicity, we will focus on
$\mathcal{X}=\operatorname{\mathbb{R}}^{d}$, but note that KSD s have been
developed for other data types [27, 28, 11, 6, 5, 25, 26]. The KSD is a
function
$\mathcal{P}(\mathcal{X})\times\mathcal{P}(\mathcal{X})\rightarrow\mathbb{R}^{+}$
which measures the similarity between two distributions
$\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}}\in\mathcal{P}(\mathcal{X})$.
Although it is not a probability metric, it is closely connected to the class
of integral pseudo-probability metrics (IPMs) [17], which measure similarity
as follows:

$\displaystyle
D_{\mathcal{F}}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=\sup_{f\in\mathcal{F}}\left|\mathbb{E}_{X\sim\operatorname{\mathbb{P}}}[f(X)]-\mathbb{E}_{Y\sim\operatorname{\mathbb{Q}}}[f(Y)]\right|$
(1)

Let $\mathcal{H}_{k}$ be the reproducing kernel Hilbert space (RKHS)
associated to the symmetric positive definite kernel
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ [2] and let
$\mathcal{H}^{d}_{k}=\mathcal{H}_{k}\times\ldots\times\mathcal{H}_{k}$ denote
the $d$-dimensional tensor product of $\mathcal{H}_{k}$. The _Langevin kernel
Stein discrepancy_ [18, 4, 15, 8] is obtained by considering an IPM with

$\displaystyle\mathcal{F}_{\operatorname{KSD}}=\\{\mathcal{S}_{\operatorname{\mathbb{P}}}[f]\>|\>\left\lVert
f\right\rVert_{\mathcal{H}_{k}^{d}}\leq 1\\}\text{ where
}\mathcal{S}_{\operatorname{\mathbb{P}}}[f](x)=f(x)\cdot\nabla\log
p(x)+\nabla\cdot f(x).$

The operator $\mathcal{S}_{\operatorname{\mathbb{P}}}$ is called the Langevin
Stein operator and $p$ is the Lebesgue density of $\operatorname{\mathbb{P}}$.
Using $\mathcal{F}_{\operatorname{KSD}}$, Equation 1 simplifies to:

$\displaystyle\operatorname{KSD}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})$
$\displaystyle:=\sup_{\|f\|_{\mathcal{H}_{k}^{d}}\leq
1}\left|\mathbb{E}_{X\sim\operatorname{\mathbb{Q}}}[\mathcal{S}_{\operatorname{\mathbb{P}}}[f](X)]\right|=\sqrt{\mathbb{E}_{X,Y\sim\operatorname{\mathbb{Q}}}[h(X,Y)]}$
$\displaystyle\text{ where }h(x,y)$ $\displaystyle={}k(x,y)\nabla\log
p(x)\cdot\nabla\log p(y)+\nabla_{x}\cdot\nabla_{y}k(x,y)$
$\displaystyle\qquad+\nabla\log p(x)\cdot\nabla_{y}k(x,y)+\nabla\log
p(y)\cdot\nabla_{x}k(x,y).$

Let $\operatorname{\mathbb{Q}}_{n}=\frac{1}{n}\sum_{i=1}^{n}\delta(x_{i})$
where $\delta(x)$ denotes a Dirac measure at $x\in\mathcal{X}$. The KSD 2 can
be straightforwardly computed in this case using a single V-statistic
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}}_{n})=\frac{1}{n^{2}}\sum_{i,j=1}^{n}h(x_{i},x_{j})$
at a cost of $O(n^{2}d)$. Under mild regularity conditions, the KSD is a
statistical divergence meaning that
$\operatorname{KSD}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})=0$ if
and only if $\operatorname{\mathbb{P}}=\operatorname{\mathbb{Q}}$; see [4,
Theorem 2.1] & [1, Proposition 1]. The KSD is convenient in the context of
unnormalised models since it depends on $\operatorname{\mathbb{P}}$ only via
$\nabla\log p$, which can be evaluated without knowledge of the normalisation
constant of $p$. We now briefly recall how it can be used for goodness-of-fit
testing and estimation.

#### Goodness-of-fit testing with KSD

In goodness-of-fit testing, we would like to test
$H_{0}:\operatorname{\mathbb{P}}=\operatorname{\mathbb{Q}}$ against
$H_{1}:\operatorname{\mathbb{P}}\neq\operatorname{\mathbb{Q}}$. A natural
approach is to compute
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}})$
and check whether this quantity is zero (i.e. $H_{0}$ holds) or not (i.e.
$H_{1}$ holds) [15, 4]. Of course, since we only have access to
$x_{1},\ldots,x_{n}$ instead of $\operatorname{\mathbb{Q}}$, this idealised
procedure is replaced by the evaluation of
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}}_{n})$.
The question then becomes whether or not this quantity is further away from
zero than expected under $H_{0}$ given the sampling error associated with a
dataset of size $n$.

To determine whether $H_{0}$ should be rejected, we need to select an
appropriate threshold $c_{\alpha}\in\operatorname{\mathbb{R}}$, which will
depend on the level of the test $\alpha\in[0,1]$. More precisely, $c_{\alpha}$
should be set to the $(1-\alpha)$-quantile of the distribution of
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}},\operatorname{\mathbb{Q}}_{n})$
under $H_{0}$. This distribution will usually be unknown a-priori, but can be
approximated using a bootstrap method. A common example is the wild bootstrap
[21, 14], which was specialised for kernel tests by [3].

#### Minimum distance estimation with KSD

As shown by [1], the KSD can also be used for parameter estimation through
_minimum distance estimation_ [24]. Given a parametric family
$\\{\operatorname{\mathbb{P}}_{\theta}\\}_{\theta\in\Theta}$ indexed by
$\Theta\subseteq\operatorname{\mathbb{R}}^{p}$. Given
$x_{1},\ldots,x_{n}\sim\operatorname{\mathbb{Q}}$, a natural estimator is

$\displaystyle\hat{\theta}_{n}:=\operatorname*{arg\,min}_{\theta\in\Theta}\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{Q}}_{n}).$

Under regularity conditions, the estimator approaches
$\theta^{*}:=\operatorname*{arg\,min}_{\theta\in\Theta}\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{Q}})$
as $n\rightarrow\infty$. The use of KSD was later extended by [9, 16, 7].
These estimators are also closely related to score-matching estimators [10];
see [1, Theorem 2] for details.

Input: $\operatorname{\mathbb{P}}_{\hat{\theta}}$,
$\operatorname{\mathbb{Q}}_{n}$, $\alpha$, $b$

for  _$k\in\\{1,\ldots,b\\}$_ do

$w^{(k)}\sim\operatorname{Rademacher}(n)$;

$\Delta^{(k)}=\frac{1}{n}\sum_{i,j=1}^{n}w^{(k)}_{i}w^{(k)}_{j}h(x_{i},x_{j})$

$c_{\alpha}=\operatorname{quantile}(\\{\Delta^{(1)},\ldots,\Delta^{(b)}\\},1-\alpha)$;

$\Delta=n\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}},\operatorname{\mathbb{Q}}_{n})$;

if $\Delta\geq c_{\alpha}$, reject the null, else do not reject.

Algorithm 1 Wild bootstrap test

Input: $\operatorname{\mathbb{P}}_{\theta}$, $\operatorname{\mathbb{Q}}_{n}$,
$\hat{\theta}_{n}$, $\alpha$, $b$

for  _$k\in\\{1,\ldots,b\\}$_ do

$\operatorname{\mathbb{Q}}_{n}^{(k)}=\frac{1}{n}\sum_{i=1}^{n}\delta_{y_{i}^{(k)}}$,
$\big{\\{}y_{i}^{(k)}\big{\\}}_{i=1}^{n}\sim\operatorname{\mathbb{P}}_{\hat{\theta}_{n}}$;

$\hat{\theta}^{(k)}_{n}=\operatorname*{arg\,min}_{\theta\in\Theta}\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{Q}}_{n}^{(k)})$;

$\Delta^{(k)}=n\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}^{(k)}_{n}},\operatorname{\mathbb{Q}}_{n}^{(k)})$;

$c_{\alpha}=\operatorname{quantile}(\\{\Delta^{(1)},\ldots,\Delta^{(b)}\\},1-\alpha)$;

$\Delta=n\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}_{n}},\operatorname{\mathbb{Q}}_{n})$;

if $\Delta\geq c_{\alpha}$, reject the null, else do not reject.

Algorithm 2 Parametric bootstrap test

##  3 Methodology

We now consider a novel composite goodness-of-fit test where both estimation
and testing is based on the $\operatorname{KSD}^{2}$ with some fixed kernel
$k$. The new test contains the two following stages:

Stage 1 (Estimation):

$\hat{\theta}_{n}=\operatorname*{arg\,min}_{\theta\in\Theta}\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta},\operatorname{\mathbb{Q}}_{n})$.

Stage 2 (Testing):

reject if
$n\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}_{n}},\operatorname{\mathbb{Q}}_{n})\geq
c_{\alpha}$, using one of Algorithm 1 or Algorithm 2.

Note that we could possibly replace stage 1 with another estimator such as
maximum likelihood estimation, but this would require knowledge of any
normalisation constant of the likelihood which may not be possible in many
cases. To implement Stage 2, we require a bootstrap algorithm to estimate the
threshold $c_{\alpha}$. A natural first approach is to use the wild bootstrap,
since it is the current gold-standard for kernel-based goodness-of-fit tests.
Algorithm 1 gives the implemention of a composite test using the wild
bootstrap, where $\alpha$ is the desired test level and $b$ is the number of
bootstrap samples. Note that in our implementation the choice of sampling
$w^{(k)}$ from a Rademacher distribution assumes that $x_{1},\ldots,x_{n}$ are
independent.

If we make the unrealistic assumption that $\hat{\theta}_{n}=\theta^{*}$, then
we are back in the setting considered by [4, 15] and the wild bootstrap works
as follows. Under $H^{C}_{0}$, as $n\rightarrow\infty$, $\Delta^{(k)}$ and
$\Delta$ converge to the same distribution, thus $c_{\alpha}$ will converge to
the $(1-\alpha)$-quantile of $\Delta$ [3, Theorem 1,2], and the test will
reject $H^{C}_{0}$ with probability $\alpha$, as desired. Under $H^{C}_{1}$,
$\Delta$ will diverge while $\Delta^{(k)}$ will converge to some distribution,
hence the probability that $\Delta>c_{\alpha})$ goes to $1$ as
$n\rightarrow\infty$.

However, it is clear that we cannot assume that $\hat{\theta}_{n}=\theta^{*}$
in the finite data case. In fact, using the wild bootstrap in this fashion for
composite tests may result in an incorrect type I error rate under
$H^{C}_{0}$, and lost power under $H^{C}_{1}$. This is because the estimation
stage of the composite test introduces a second source of error which this
approach does not take account of when computing $c_{\alpha}$. The two sources
of error that the test encounters are as follows. Recall that in an idealised
setting we would use
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta^{*}},\operatorname{\mathbb{Q}})$
as the test statistic. The first source of error is introduced because we must
estimate this statistic with
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta^{*}},\operatorname{\mathbb{Q}}_{n})$,
as we only have access to a sample from $\operatorname{\mathbb{Q}}$. This
source of error also occurs in non-composite tests, and is accounted for
correctly by the wild bootstrap. The second source of error is specific to
composite tests, and occurs because we must further estimate
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\theta^{*}},\operatorname{\mathbb{Q}}_{n})$
with
$\operatorname{KSD}^{2}(\operatorname{\mathbb{P}}_{\hat{\theta}_{n}},\operatorname{\mathbb{Q}}_{n})$,
as we do not have access to $\theta^{*}$. Algorithm 1 fixes the parameter
estimate $\hat{\theta}_{n}$ and then applies the bootstrap, thus failing to
take account of the error in in $\hat{\theta}_{n}$ and potentially computing
an incorrect threshold. We can also view Algorithm 1 as a non-composite
goodness-of-fit test against the wrong null hypothesis. By estimating the
parameter and then applying the test, the test is not evaluating
$H_{0}:\operatorname{\mathbb{P}}_{\theta^{*}}=\operatorname{\mathbb{Q}}$, but
instead
$H_{0}^{\prime}:\operatorname{\mathbb{P}}_{\hat{\theta}_{n}}=\operatorname{\mathbb{Q}}$.

Figure 2 demonstrates the impact of ignoring this error. It shows the power of
our test with
$H^{C}_{0}:\operatorname{\mathbb{P}}_{\theta}=\operatorname{\mathcal{N}}(\mu,1)$
for $\mu\in\mathbb{R}$ and
$\operatorname{\mathbb{Q}}=\operatorname{\mathcal{N}}(1.3,\sigma^{2})$, as
$\sigma^{2}$ varies, see Appendix B for details. We can see that the wild
bootstrap test has lower power for most values of $\sigma^{2}$, including when
$\sigma^{2}=1$, thus $H^{C}_{0}$ holds, and it corresponds to the type I
error. It also demonstrates that the error depends on the size of $n$. The
estimator is consistent: as $n$ becomes large we expect the estimate of the
parameter to converge on the true value, thus minimizing the impact of the
error. However, for smaller $n$ it is necessary to use an alternative approach
which is able to take account of the additional error.

To take account of both types of error we apply the parametric bootstrap
described in Algorithm 2. This bootstrap approximates the distribution of the
test statistic under $H^{C}_{0}$ by repeatedly resampling the observations and
re-estimating the parameter, thus taking account of the error in the estimate
of the parameter. In comparison to the wild bootstrap, this is substantially
more computationally intensive because it requires repeatedly computing the
kernel matrix on fresh data, whereas the wild bootstrap only computes the
kernel matrix once. However, when $n$ is large and this extra computation is
likely to be an issue, the size of the estimation error should also be reduced
as the estimator converges on the true value of the parameter. Thus, under
this high data regime, the wild and parametric bootstraps may achieve
comparable powers, and it may be reasonable to use the cheaper wild bootstrap.
In the setting considered in Figure 2 we find that this is the case, with the
parametric bootstrap performing substantially better than the wild bootstrap
for $n=10$, but comparably for $n=200$.

##  4 Illustration: The Kernel Exponential Family Model

Figure 1:  Power of of the test when using the wild and parametric bootstraps
for $n=10$ (solid) and $n=200$ (dashed). The dashed horizontal line shows the
test level, $\alpha=0.05$. The error bars show one standard error over $4$
random seeds.

Figure 2:  Goodness-of-fit of a kernel exponential family model with
increasing $p$. The grey histogram shows the dataset, and the lines the
modelled density. Dashed lines indicate that the test rejected $H^{C}_{0}$,
and solid that it did not.

We apply our method to a density estimation task from [16], to test whether a
robust method was really necessary. The model is in the kernel exponential
family $p_{\theta}(x)\propto q(x)\exp(f(x))$, where $f$ is assumed to belong
to some RKHS but approximated by a finite sum of basis functions, and
$q(x)=\operatorname{\mathcal{N}}(0,3^{2})$ is a reference density.
Specifically, we follow [22] and consider
$f(x)=\sum_{i=1}^{p}\theta_{i}\phi_{i}(x)$, with
$\phi_{i}(x)=x^{i}/\sqrt{i!}\exp\left(-x^{2}/2\right),\theta\in\operatorname{\mathbb{R}}^{p}$.
The dataset is comprised of the velocities of $82$ galaxies [19, 20]. An open
question is how large $p$ needs to be for the model to have enough capacity to
fit the data. [16] set $p=25$, and we test this assumption. See Appendix B for
the full experiment configuration. Figure 2 shows the fit of the model for
increasing values of $p$, and whether the test rejected the null hypothesis.
We find that the test rejects $H^{C}_{0}$ for $p=1,2,3$, but does not reject
for $p=4,5,25$. This suggests that $p=25$ is a suitable choice, though it
would be reasonable to use a smaller value of $p$ which would decrease the
computational cost of inference.

##  5 Conclusion

In this initial work, we introduced a flexible composite goodness-of-fit test
which requires only an unnormalized density. In the full version of this paper
we will consider a general framework for composite goodness-of-fit tests,
which will also include a tests based on the maximum mean discrepancy. This
will allow us to tackle a much wider set of problems, including composite
testing for generative models with no computable density function. We will
also include theoretical results showing that the test statistic converges
under $H^{C}_{0}$, the test is consistent under $H^{C}_{1}$, and that the wild
and parametric bootstraps behave appropriately.

## References

  * [1] Alessandro Barp, Francois-Xavier Briol, Andrew Duncan, Mark Girolami and Lester Mackey  “Minimum Stein Discrepancy Estimators”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [2] A. Berlinet and C. Thomas-Agnan  “Reproducing Kernel Hilbert Spaces in Probability and Statistics”, 2004
  * [3] Kacper Chwialkowski, Dino Sejdinovic and Arthur Gretton  “A Wild Bootstrap for Degenerate Kernel Tests”  In _NIPS_ , 2014, pp. 3608–3616
  * [4] Kacper Chwialkowski, Heiko Strathmann and Arthur Gretton  “A Kernel Test of Goodness of Fit”  In _International Conference on Machine Learning_ , 2016, pp. 2606–2615  PMLR
  * [5] Tamara Fernandez, Nicolas Rivera, Wenkai Xu and Arthur Gretton  “Kernelized Stein Discrepancy Tests of Goodness-of-Fit for Time-to-Event Data”  In _International Conference on Machine Learning_ , 2020, pp. 3112–3122  PMLR
  * [6] T. Fernández and A. Gretton  “A Maximum-Mean-Discrepancy Goodness-of-Fit Test for Censored Data”  In _Artificial Intelligence and Statistics_ , 2019
  * [7] W. Gong, Y. Li and J.. Hernández-Lobato  “Sliced Kernelized Stein Discrepancy”  In _International Conference on Learning Representations_ , 2021
  * [8] Jackson Gorham and Lester Mackey  “Measuring Sample Quality with Steins Method”  In _Advances in Neural Information Processing Systems_ 28, 2015
  * [9] W. Grathwohl, K.-C. Wang, J.-H. Jacobsen, D. Duvenaud and R. Zemel  “Learning the Stein Discrepancy for Training and Evaluating Energy-Based Models Without Sampling”  In _International Conference on Machine Learning_ , 2020, pp. 9485–9499
  * [10] A. Hyvärinen  “Estimation of Non-Normalized Statistical Models by Score Matching”  In _Journal of Machine Learning Research_ 6, 2006, pp. 695–708
  * [11] Heishiro Kanagawa, Wittawat Jitkrittum, Lester Mackey, Kenji Fukumizu and Arthur Gretton  “A Kernel Stein Test for Comparing Latent Variable Models”, 2020  arXiv:1907.00586
  * [12] Jérémie Kellner and Alain Celisse  “A One-Sample Test for Normality with Kernel Methods”  In _Bernoulli_ 25.3, 2019, pp. 1816–1837
  * [13] Erich L. Lehmann and Joseph P. Romano  “Testing Statistical Hypotheses”, Springer Texts in Statistics, 2005
  * [14] Anne Leucht and Michael H. Neumann  “Dependent Wild Bootstrap for Degenerate U- and V-Statistics”  In _Journal of Multivariate Analysis_ 117, 2013, pp. 257–280
  * [15] Qiang Liu, Jason Lee and Michael Jordan  “A Kernelized Stein Discrepancy for Goodness-of-Fit Tests”  In _Proceedings of the 33rd International Conference on Machine Learning_ 48, Proceedings of Machine Learning Research, 2016, pp. 276–284
  * [16] Takuo Matsubara, Jeremias Knoblauch, François-Xavier Briol and Chris J. Oates  “Robust Generalised Bayesian Inference for Intractable Likelihoods”, 2021  arXiv:2104.07359
  * [17] A. Muller  “Integral Probability Metrics and Their Generating Classes of Functions”  In _Advances in Applied Probability_ 29.2, 1997, pp. 429–443
  * [18] Chris Oates and Mark Girolami  “Control Functionals for Quasi-Monte Carlo Integration”  In _Proceedings of the 19th International Conference on Artificial Intelligence and Statistics_ 51, Proceedings of Machine Learning Research, 2016, pp. 56–65
  * [19] M. Postman, J.. Huchra and M.. Geller  “Probes of Large-Scale Structure in the Corona Borealis Region.”  In _The Astronomical Journal_ 92, 1986, pp. 1238–1247
  * [20] Kathryn Roeder  “Density Estimation with Confidence Sets Exemplified by Superclusters and Voids in the Galaxies”  In _Journal of the American Statistical Association_ 85.411, 1990, pp. 617–624
  * [21] Xiaofeng Shao  “The Dependent Wild Bootstrap”  In _Journal of the American Statistical Association_ 105.489, 2010, pp. 218–235
  * [22] I. Steinwart, D. Hush and C. Scovel  “An Explicit Description of the Reproducing Kernel Hilbert Spaces of Gaussian RBF Kernels”  In _IEEE Transactions on Information Theory_ 52.10, 2006, pp. 4635–4643
  * [23] Winfried Stute, Wenceslao Gonzáles Manteiga and Manuel Presedo Quindimil  “Bootstrap Based Goodness-of-Fit-Tests”  In _Metrika_ 40.1, 1993, pp. 243–256
  * [24] J. Wolfowitz  “The Minimum Distance Method”  In _The Annals of Mathematical Statistics_ 28.1, 1957, pp. 75–88
  * [25] W. Xu and T. Matsuda  “A Stein Goodness-of-Fit Test for Directional Distributions”  In _Artificial Intelligence and Statistics_ , 2020
  * [26] W. Xu and G. Reinert  “A Stein Goodness of Fit Test for Exponential Random Graph Models”  In _Artificial Intelligence and Statistics_ , 2021
  * [27] J. Yang, Q. Liu, V. Rao and J. Neville  “Goodness-of-Fit Testing for Discrete Distributions via Stein Discrepancy”  In _International Conference on Machine Learning_ , 2018, pp. 5561–5570
  * [28] J. Yang, V. Rao and J. Neville  “A Stein-Papangelou Goodness-of-Fit Test for Point Processes”  In _International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 226–235

##  Appendix A Closed-form expression for KSD estimator

For exponential family models, such as the Gaussian and kernel exponential
families considered in this paper, there is a closed-form expression for the
KSD estimator. Let the density of a model in the exponential family be

$\displaystyle p_{\theta}(x)=\exp(\eta(\theta)\cdot t(x)-a(\theta)+b(x))\,,$

with $\eta:\Theta\to\mathbb{R}^{k}$ an invertible map,
$t:\mathbb{R}^{d}\to\mathbb{R}^{k}$ any sufficient statistic for some
$k\in\mathbb{N}_{1}$, $a:\Theta\to\mathbb{R}$ and
$b:\mathbb{R}^{d}\to\mathbb{R}$. Then the KSD estimator of theta is given by

$\displaystyle\hat{\theta}_{n}:$
$\displaystyle=\operatorname*{arg\,min}\operatorname{KSD}^{2}(\operatorname{\operatorname{\mathbb{P}}}_{\theta}\|\operatorname{\operatorname{\mathbb{Q}}}_{n})=\eta^{-1}\left(-\frac{1}{2}\Lambda_{n}^{-1}v_{n}\right)\,,$

where $\Lambda_{n}\in\mathbb{R}^{k\times k}$ and $\nu_{n}\in\mathbb{R}^{k}$
are defined as
$\Lambda_{n}:=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\Lambda(x_{i},x_{j})$
and $\nu_{n}:=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\nu(x_{i},x_{j})$.
These are based on functions
$\Lambda:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{k\times k}$ and
$\nu:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}^{k}$ which depend on the
specific model, defined as

$\displaystyle\Lambda(x,x^{\prime}):$ $\displaystyle=k(x,x^{\prime})\nabla
t(x)\nabla t(x^{\prime})^{\top}$ $\displaystyle\nu(x,x^{\prime}):$
$\displaystyle=k(x,x^{\prime})\nabla b(x)\nabla t(x^{\prime})^{\top}+\nabla
t(x)\nabla_{x^{\prime}}k(x,x^{\prime})+k(x,x^{\prime})\nabla
b(x^{\prime})\nabla t(x)^{\top}+\nabla
t(x^{\prime})\nabla_{x}k(x^{\prime},x)\,.$

The detailed derivation of this result can be found in [1, Appendix D3] or
[16].

##  Appendix B Experiment details

###  B.1 Figure 2

For each datapoint in the plot, the power is computed as follows:

  1. 1.

For four random seeds:

    1. (a)

For $2000$ repeats:

      1. i.

Sample $n$ points from
$\operatorname{\mathbb{Q}}=\operatorname{\mathcal{N}}(1.3,\sigma^{2})$

      2. ii.

Run the test

    2. (b)

Compute the power, the fraction of repeats during which the null hypothesis is
rejected

  2. 2.

Compute the mean and standard error over the seeds

We find the standard error to be sufficiently small that it is not visible in
the figure.

#### Estimator

We estimate only the mean, which we denote by $\theta$. The variance,
$\sigma^{2}$, is known and specified by the null hypothesis. We use the
closed-form KSD estimator as defined above, with

$\displaystyle\eta(\theta)=\frac{\mu}{\sigma^{2}},\qquad t(x)=x,\qquad
b(x)=\log\frac{1}{\sqrt{2\pi\sigma^{2}}}-\frac{x^{2}}{2\sigma^{2}}\,.$

#### Kernel

Gaussian kernel:

$\displaystyle k(x_{1},x_{2})=\exp\left(-\frac{\left\lVert
x_{1}-x_{2}\right\rVert_{2}^{2}}{2l^{2}}\right)\,.$

We choose $l=0.7$.

#### Bootstrap and test

The wild bootstrap (Algorithm 1) is configured with $b=500$. The parametric
bootstrap (Algorithm 2) is configured with $b=300$. We set $\alpha=0.05$ in
both cases.

###  B.2 Figure 2

#### Data normalisation

Following [16], we normalize the dataset by

$\displaystyle
y^{\prime}_{i}=\frac{y_{i}-\operatorname{mean}(y_{1},\ldots,y_{n})}{\frac{1}{2}\operatorname{std}(y_{1},\ldots,y_{n})}\qquad\forall
i\in\\{1,\ldots,n\\}\,.$

#### Estimator

We use the closed-form KSD estimator as defined above, with

$\displaystyle\eta(\theta)=\theta,\qquad
t(x)=(\phi_{1}(x),\ldots,\phi_{p}(x))^{\top},\qquad b(x)=\log q_{0}(x)\,.$

#### Kernel

IMQ kernel:

$\displaystyle k(x_{1},x_{2})=\left(1+\frac{\left\lVert
x_{1}-x_{2}\right\rVert_{2}^{2}}{l^{2}}\right)^{-\frac{1}{2}}\,.$

We select $l$ using the median heuristic,

$\displaystyle l_{\textrm{med}}$
$\displaystyle=\sqrt{\operatorname{median}(d_{ij}/2)}\textrm{ with
}d_{ij}=\left\lVert y_{i}-y_{j}\right\rVert_{2}^{2}\textrm{ for all
}i,j=1,\ldots,n\,,$

which in this case results in $l\approx 0.9$.

#### Bootstrap and test

Parametric bootstrap (Algorithm 2), with $b=400$ and $\alpha=0.05$.