# Large Deviation Analysis of Score-based Hypothesis Testing

Enmao Diao, Taposh Banerjee, _Member, IEEE_ , and Vahid Tarokh, _Fellow, IEEE_
Enmao Diao and Vahid Tarokh are with the Department of Electrical and Computer
Engineering, Duke University, Durham, NC 27708, USA. Taposh Banerjee is with
the Department of Industrial Engineering, University of Pittsburgh, PA 15213,
USA. Vahid Tarokh was supported in part by Air Force Research Lab grant number
FA-8750-20-2-0504. Taposh Banerjee was supported in part by the U.S. Army
Research Lab under grant W911NF2120295.

###### Abstract

Score-based statistical models play an important role in modern machine
learning, statistics, and signal processing. For hypothesis testing, a score-
based hypothesis test is proposed in [1]. We analyze the performance of this
score-based hypothesis testing procedure and derive upper bounds on the
probabilities of its Type I and II errors. We prove that the exponents of our
error bounds are asymptotically (in the number of samples) tight for the case
of simple null and alternative hypotheses. We calculate these error exponents
explicitly in specific cases and provide numerical studies for various other
scenarios of interest.

##  I Introduction

Score matching [2, 3] is a procedure that has emerged due to its ability to
outperform likelihood-based benchmarks in image generations [4, 5]. It has
been shown that it is possible to efficiently model $\nabla\log p(x)$ using
deep neural networks (DNNs), where $p(x)$ denotes the data generating density.
However, it is computationally non-trivial to calculate even un-normalized
versions of $p(x)$ from these trained DNNs. For log-likelihood-based
hypothesis testing and change detection, there exists a well-established
literature. Recently change detection and hypothesis testing for score-based
systems have been studied. In particular, score-based hypothesis testing [1]
was proposed as a procedure to decide whether or not to reject a hypothesis.
We note that the likelihood ratio test (LRT) is a standard method commonly
used for hypothesis testing, and the celebrated Neyman-Pearson lemma gives the
uniformly most powerful optimality property of LRT for simple hypotheses
testing. This implies that score-based hypothesis testing cannot outperform
the LRT test for simple hypotheses testing when the densities of data under
the null and alternative hypotheses are exactly known. However, the evaluation
of exact likelihoods may be computationally cumbersome (or even intractable)
for various modern statistical models including graphical models [6], energy-
based models [7] and deep generative models [8]. This gives some importance to
score-based hypothesis testing [1] and the analysis of its performance. This
motivates our work in this paper, where we derive upper bounds on the
probabilities of its Type I and II errors for score-based hypothesis testing,
and establish that the exponents of our error bounds are asymptotically (in
the number of samples) tight for the case of simple null and alternative
hypotheses.

The outline of our paper is given next. In Section II, we will motivate the
problem to be considered in this paper, review the related work, and present
the score-based test for binary hypothesis testing [1]. We then derive upper
bounds on the probabilities of its Type I and II errors for finite sample size
in Section III. In Section IV, we show that these bounds are asymptotically
tight using large deviation theory [9, 10]. In Section V, we discuss
simulation methods for estimating these error exponents. In Section VI, we
provide explicit calculation for the error exponent in the multivariate
Gaussian case. In Section VII, we show the accuracy of these error bounds for
experiments performed on real and simulated data. We will make our concluding
remarks in Section VIII.

##  II Score-based Hypothesis Testing for Unnormalized Statistical Models

###  II-A Binary Hypothesis Testing

Let ${x}\in\mathcal{X}\subseteq\mathbb{R}^{d}$ be the realization of a random
vector $X$. We denote
$\mathbf{X}_{n}\overset{\Delta}{=}\\{{X}_{1},\dots,{X}_{n}\\}$ as independent
and identically distributed (i.i.d.) observations from an unknown distribution
$P$ with probability density $p$. For two probability measures $P_{\infty}$
and $P_{1}$ with respective densities $p_{\infty}$ and $p_{1}$, we investigate
the simple binary hypotheses testing problem:

$\displaystyle\mathcal{H}_{0}:p=p_{\infty}\quad\text{against}\quad\mathcal{H}_{1}:p=p_{1}.$
(1)

It is well-known that the optimal test in Bayesian, minimax, and variational
settings is the likelihood ratio test (LRT):

$\delta_{L}(X)=\begin{cases}1,&\text{ if
}S_{L}(\mathbf{X}_{n},p_{\infty})-S_{L}(\mathbf{X}_{n},p_{1})>T\\\
0&\quad\text{ otherwise}.\end{cases}$ (2)

Here,

$S_{L}(\mathbf{X}_{n},p)=-\log\prod_{i=1}^{n}p({X}_{i}).$

is the log score of the sample $(X_{1},\dots,X_{n})$. The choice of the
threshold $T$ depends on the problem formulation. In the variational setting
of Neymann and Pearson, the threshold is chosen to set the type-I error. The
performance and optimality properties of this test can be found in most
standard texts; see, for example, [11, 12]. When the densities $p_{\infty}$
and $p_{1}$ are not known, the LRT statistic is replaced by a suitable
statistic depending on the assumptions on the data. Such families of tests
include the generalized likelihood ratio (GLR) test, mixture tests, and
nonparametric tests [12, 13].

###  II-B Limitations of the LRT for Score-Based Models

In modern machine learning applications, two new classes of models have
emerged:

  1. 1.

Unnormalized statistical models: In these models, we know the densities
$p_{\infty}$ and $p_{1}$ within normalizing constants. Specifically, we have

$\begin{split}p_{\infty}(x)=\frac{\tilde{p}_{\infty}(x)}{Z_{\infty}},\quad\text{
and }\quad p_{1}(x)=\frac{\tilde{p}_{1}(x)}{Z_{1}}.\end{split}$ (3)

Here

$\begin{split}Z_{\infty}=\int_{x}\tilde{p}_{\infty}(x)dx,\quad\text{ and
}\quad Z_{1}=\int_{x}\tilde{p}_{1}(x)dx.\end{split}$ (4)

are normalizing constants that are hard (or even impossible) to calculate by
numerical integration. The unnormalized models $\tilde{p}_{\infty}(x)$ and
$\tilde{p}_{1}(x)$ are known in precise functional forms. Examples include
continuous-valued Markov random fields or undirected graphical models which
are used for image modeling [2, 14]. In this case, the LRT can still be
implemented. This is because the LRT can be written as

$\begin{split}\delta_{L}(X)&=\begin{cases}1,&\text{ if
}\prod_{i=1}^{n}\frac{p_{1}({X}_{i})}{p_{\infty}({X}_{i})}>T\\\ 0&\quad\text{
otherwise},\end{cases}\\\ &=\begin{cases}1,&\text{ if
}\prod_{i=1}^{n}\frac{\tilde{p}_{1}({X}_{i})}{\tilde{p}_{\infty}({X}_{i})}>T_{1}\\\
0&\quad\text{ otherwise}.\end{cases}\end{split}$ (5)

Here $T_{1}=\left(\frac{Z_{1}}{Z_{\infty}}\right)^{n}T$. If $T_{1}$ is chosen
to satisfy a constraint on the type-I error (which can be done by sampling
from $p_{\infty}$ using its unnormalized version $\tilde{p}_{\infty}$), then
the LRT test can still be implemented.

  2. 2.

Score-based models: In many modern machine learning applications, even
$\tilde{p}_{\infty}(x)$ and $\tilde{p}_{1}(x)$ are unknown. But, we may learn
the scores,

$\nabla_{x}\log p_{\infty}(x)\quad\text{ and }\quad\nabla_{x}\log p_{1}(x),$

from data. This is possible using the idea of score-matching. Specifically,
these scores can be learned using a deep neural network [2, 4, 14, 15, 16,
17]. We note that a score-based model is also unnormalized where the exact
form of the unnormalized function is hard to estimate. In this case, the LRT
cannot be implemented and a fresh approach is needed.

###  II-C Score-Based Hypothesis Testing

A score-based approach to hypothesis testing was taken in [1]. This test was
based on the concept of Hyvärinen Score [2]. To define this score and the
corresponding test, let $p$ and $q$ be two probability densities defined on
$\mathcal{X}$. The Fisher divergence between $p$ and $q$ is defined as

$\displaystyle\mathbb{D}_{\texttt{F}}(p\;||\;q)\overset{\Delta}{=}\mathbb{E}_{{X}\sim
p}\left[\frac{1}{2}\left\|\nabla_{x}\log p({X})-\nabla_{{x}}\log
q({X})\right\|_{2}^{2}\right],$ (6)

whenever the integral is well defined. Under some mild regularity conditions
on $p$ and $q$ [2], the Fisher divergence can be rewritten as

$\displaystyle\mathbb{D}_{\texttt{F}}(p\;||\;q)=\mathbb{E}_{{X}\sim
p}\left[\frac{1}{2}\left\|\nabla_{{x}}\log
p({X})\right\|_{2}^{2}+S_{\texttt{H}}({X},q)\right],$ (7)

with

$\displaystyle
S_{\texttt{H}}({X},q)\overset{\Delta}{=}\frac{1}{2}\left\|\nabla_{{x}}\log
q({X})\right\|_{2}^{2}+\Delta_{{x}}\log q({X}),$ (8)

where $\Delta_{{x}}=\sum_{i=1}^{d}\frac{\partial^{2}}{\partial x_{i}^{2}}$
denotes the Laplacian operator with respect to
${x}=(x_{1},\cdots,x_{d})^{\top}$. The score $S_{\texttt{H}}({X},q)$ is the
Hyvärinen score [2] and is a counterpart of the logarithmic score function
$S_{L}(X,p)=-\log q({X})$ that corresponds to the widely used negative log-
likelihood. We further use the notation ${S}_{\texttt{H}}(\mathbf{X}_{n},p)$
to denote

$S_{\texttt{H}}(\mathbf{X}_{n},p)\overset{\Delta}{=}\frac{1}{n}\sum_{i=1}^{n}S_{\texttt{H}}({X}_{i},p).$
(9)

The score-based test for binary hypothesis testing proposed in [1] is given by

$\delta(\mathbf{X}_{n})=\begin{cases}1,&\text{ if
}S_{\texttt{H}}(\mathbf{X}_{n},p_{\infty})-S_{\texttt{H}}(\mathbf{X}_{n},p_{1})>T\\\
0&\quad\text{ otherwise}.\end{cases}$ (10)

In this paper, we analyze the type-I and type-II error probabilities of this
test as $n\to\infty$ using large deviation theory [9, 10]. In the test of the
paper, we drop the subscript $H$ from the $S_{\texttt{H}}(\mathbf{X}_{n},p)$
and simply refer to the Hyvärinen score by $S(\mathbf{X},p)$ for a collection
of $n$ data points and by $S({X},p)$ for a single data point.

Recently, there have been some attempts to perform out-of-distribution (OOD)
detection based on score-matching estimates [18, 19]. Note that in OOD
detection, the in-distribution samples are assumed to follow a distribution
$p(\cdot)$, while the OOD samples are distributed according to another
distribution $q(\cdot)$. For an observation $x$, hypothesis
$\mathcal{H}_{0}:x\sim p$ is then tested against $\mathcal{H}_{1}:x\not\sim
p$. In [18], the authors utilized the norm of the gradient of logarithmic
likelihood at multiple noise scales for anomaly detection. The authors of [19]
estimated the Fisher divergence between the null and alternative distributions
for feature shift detection assuming that the data is generated according to
the null distribution. In [1], the scores of null distribution and alternative
distributions were used to propose a score-based hypothesis test. As we will
be analyzing the performance of this method, we will review its required
background and subsequently discuss it below.

##  III Fixed-Sample Error Analysis of Score-Based Test for Binary Hypothesis
Testing

Consider the single-sample score-based simple hypotheses test between
densities $p_{\infty}$ and $p_{1}$ using $n$ data points
$\mathbf{X}_{n}=(X_{1},\dots,X_{n})$ (after supressing the subscript $H$ in
(10)):

$\delta_{n}\overset{\Delta}{=}\delta(\mathbf{X}_{n})=\begin{cases}1,&\text{ if
}S(\mathbf{X}_{n},p_{\infty})-S(\mathbf{X}_{n},p_{1})>T\\\ 0&\quad\text{
otherwise},\end{cases}$ (11)

where recall that $S(\mathbf{X}_{n},p)=\frac{1}{n}\sum_{i=1}^{n}S({X}_{i},p)$.
Thus, choosing $\delta_{n}=1$ is interpreted as the selection of the
alternative hypothesis $\mathcal{H}_{1}:p=p_{1}$. The type I error
(probability of false alarm) is then given by

$\alpha_{n}(\delta_{n})\overset{\Delta}{=}\mathsf{P}_{\infty}(\delta_{n}=1)=\mathsf{P}_{\infty}\left(S(\mathbf{X}_{n},p_{\infty})-S(\mathbf{X}_{n},p_{1})>T\right).$
(12)

This probability can be estimated using Chernoff’s bound: for $\theta>0$

$\begin{split}\alpha_{n}(\delta_{n})&=\mathsf{P}_{\infty}\left(S(\mathbf{X}_{n},p_{\infty})-S(\mathbf{X}_{n},p_{1})>T\right)\\\
&=\mathsf{P}_{\infty}\left(\frac{1}{n}\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1})>T\right)\\\
&=\mathsf{P}_{\infty}\left(\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1})>nT\right)\\\
&=\mathsf{P}_{\infty}\left(e^{\theta(\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1})))}>e^{n\theta
T}\right)\\\ &\leq e^{-n\theta
T}\mathsf{E}_{\infty}\left[e^{\theta(\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1})))}\right]\\\
&=\left[e^{-\theta
T}e^{\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]}\right]^{n}\\\
&=\exp\left[n\left(\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]-\theta
T\right)\right],\end{split}$ (13)

where $\mathsf{E}_{\infty}(\cdot)$ denotes expectation with respect to
$\mathsf{P}_{\infty}$.

Since this inequality is true for every $\theta>0$, we can take an infimum
over $\theta$ on the right-hand side.

$\begin{split}\alpha_{n}(\delta_{n})&\leq\inf_{\theta\geq
0}\exp\left(n\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]-\theta
T\right)\\\ &=\exp\left(n\inf_{\theta\geq
0}\left(\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]-\theta
T\right)\right)\\\ &=\exp\left(-n\sup_{\theta\geq 0}\left(\theta
T-\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]\right)\right)\\\
\end{split}$ (14)

Define

$\phi(\theta)=\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,P_{\infty})-S(X,P_{1}))}\right].$

It is well-known that $\phi(\theta)$ is a convex function. Define the Legendre
transformation of the convex function $\phi(\theta)$ as

$\phi^{*}(T)=\sup_{\theta\geq 0}\left[\theta T-\phi(\theta)\right].$

Thus,

$\alpha_{n}(\delta_{n})\leq e^{-n\phi^{*}(T)},$ (15)

for all $n\geq 1$. Note that $\phi^{*}(T)$ is strictly positive if $T$ is
larger than the the slope of $\phi(\theta)$ at $\theta=0$, which is

$\frac{d}{d\theta}\phi(\theta)\bigg{|}_{\theta=0}=\frac{d}{d\theta}\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]\bigg{|}_{\theta=0}=-\mathbb{D}_{\texttt{F}}(p_{\infty}\;\|\;p_{1})<0.$

Thus, for $T>-\mathbb{D}_{\texttt{F}}(p_{\infty}\;\|\;p_{1})$,

$\alpha_{n}(\delta_{n})\leq e^{-n\phi^{*}(T)}<1$ (16)

for all $n\geq 1$. Thus,

$\frac{\log\alpha_{n}(\delta_{n})}{n}\leq-\phi^{*}(T).$ (17)

In the next section, we show that this bound $-\phi^{*}(T)$ is asymptotically
tight using large deviation theory.

We note that a similar argument can be made for the type II error by changing
the roles of $P_{1}$ and $P_{\infty}$. Specifically, we have

$\beta_{n}(\delta_{n})=\mathsf{P}_{1}\left(S(\mathbf{X}_{n},p_{1})-S(\mathbf{X}_{n},p_{\infty})>-T\right)\leq\exp\left(-n\sup_{\theta\geq
0}\left(-\theta
T-\log\mathsf{E}_{1}\left[e^{\theta(S(X,p_{1})-S(X,p_{\infty}))}\right]\right)\right),$

for all $n\geq 1$, where $\mathsf{E}_{1}(\cdot)$ denotes expectation with
respect to $\mathsf{P}_{1}$. We conclude that

$\frac{\log\beta_{n}(\delta_{n})}{n}\leq-\sup_{\theta\geq 0}\left(-\theta
T-\log\mathsf{E}_{1}\left[e^{\theta(S(X,P_{1})-S(X,P_{\infty})}\right]\right).$

##  IV Large Deviation Analysis of the Score-Based Test

By Cramer’s theorem [9], the type-I error has a positive error exponent given
by

$\begin{split}\frac{1}{n}\log\mathsf{P}_{\infty}(\delta_{n}=1)&=\frac{1}{n}\log\mathsf{P}_{\infty}\left(\frac{1}{n}\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1}))>T\right)\\\
&\to\inf_{\theta\geq
0}\left(\log\mathsf{E}_{\infty}\left[e^{\theta(S(X_{1},p_{\infty})-S(X_{1},p_{1}))}\right]-\theta
T\right)\\\ &=-\sup_{\theta\geq 0}\left(\theta
T-\log\mathsf{E}_{\infty}\left[e^{\theta(S(X_{1},p_{\infty})-S(X_{1},p_{1}))}\right]\right)\\\
&=-\phi^{*}(T).\end{split}$

Thus, the exponent we calculated in (17) for a sample of size $n$ (using the
Chernoff bound) turns out to be tight in the asymptotic regime as
$n\rightarrow\infty$.

Similarly, by Cramer’s theorem again, the type-II error has a positive error
exponent given by

$\begin{split}\frac{1}{n}\log\mathsf{P}_{1}(\delta_{n}=0)&=\frac{1}{n}\log\mathsf{P}_{1}\left(\frac{1}{n}\sum_{i=1}^{n}(S(X_{i},p_{\infty})-S(X_{i},p_{1}))<T\right)\\\
&=\frac{1}{n}\log\mathsf{P}_{1}\left(\frac{1}{n}\sum_{i=1}^{n}(S(X_{i},p_{1})-S(X_{i},p_{\infty}))>-T\right)\\\
&\to\inf_{\theta\geq
0}\left(\log\mathsf{E}_{1}\left[e^{\theta(S(X_{1},p_{1})-S(X_{1},p_{\infty}))}\right]+\theta
T\right)\\\ &=-\sup_{\theta\geq 0}\left(-\theta
T-\log\mathsf{E}_{1}\left[e^{\theta(S(X_{1},p_{1})-S(X_{1},p_{\infty}))}\right]\right).\end{split}$

Let

$-\mathbb{D}_{\texttt{F}}(p_{\infty}\;\|\;p_{1})<T<\mathbb{D}_{\texttt{F}}(p_{1}\;\|\;p_{\infty}).$

Then,

$-T>-\mathbb{D}_{\texttt{F}}(p_{1}\;\|\;p_{\infty}),$

and this error exponent is also positive.

##  V Estimating the Error Exponent

The error exponent $\phi^{*}(T)$ can be estimated empirically by using
numerical simulation to design the hypothesis test. Since we have access to
the score $\nabla\log p_{\infty}(x)$, we can sample from it using Markov Chain
Monte Carlo (MCMC) techniques such as Metropolis Adjusted Langevin Algorithm
(MALA), Hamilitonian Monte Carlo (HMC), etc. Let $X_{1},X_{2},\dots,X_{m}$ be
$m$ samples generated using such an MCMC Algorithm. Then the function of
$\theta$ given by

$\theta
T-\log\mathsf{E}_{\infty}\left[e^{\theta(S(X_{1},p_{\infty})-S(X_{1},p_{1}))}\right]$

can be empirically estimated by

$\theta
T-\log\frac{1}{m}\sum_{k=1}^{m}\left[e^{\theta(S(X_{k},p_{\infty})-S(X_{k},p_{1}))}\right].$

The above expression can be optimized over $\theta$ using standard gradient
descent algorithms in order to estimate $\phi^{*}(T)$. The error exponent for
type II error can be numerically calculated in an analogous manner.

##  VI Explicit Calculation For Specific Cases

We note that for some specific distributions, the above error exponents can be
calculated in closed form. As an example, we consider the standard simple
Gaussian hypothesis testing with where the null and alternative hypothesis
have the same co-variance matrices. Without loss of generality, we can assume
that $P_{\infty}$ is multidimensional Gaussian $N(0,\Sigma)$ with mean zero
and co-variance matrix $\Sigma$ and $P_{1}$ is $N(\mu,\Sigma)$. Then by direct
calculation:

$S(X,p_{\infty})-S(X,p_{1})=-\frac{\mu^{T}}{2}\Sigma^{-2}(X-\frac{\mu}{2}).$

Another straightforward calculation gives

$\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]=\mathsf{E}_{\infty}\left[e^{-\theta\frac{\mu^{T}}{2}\Sigma^{-2}(X-\frac{\mu}{2})}\right]=e^{\theta\frac{\mu^{T}\Sigma^{-2}\mu}{4}}e^{\theta^{2}\frac{\mu^{T}\Sigma^{-3}\mu}{8}}.$

We thus have

$\theta
T-\log\mathsf{E}_{\infty}\left[e^{\theta(S(X,p_{\infty})-S(X,p_{1}))}\right]=\theta
T-\theta\frac{\mu^{T}\Sigma^{-2}\mu}{4}-\theta^{2}\frac{\mu^{T}\Sigma^{-3}\mu}{8}.$

If $T<\frac{\mu^{T}\Sigma^{-2}\mu}{4}$ then the maxima is achieved at
$\theta=0$ and the type I error exponent is given by $\phi^{*}(T)=0$.
Otherwise, the maxima is achieved at
$\theta=\frac{4T-\mu^{T}\Sigma^{-2}\mu}{\mu^{T}\Sigma^{-3}\mu}$ and the type I
error exponent is given by

$\phi^{*}(T)=\frac{(4T-\mu^{T}\Sigma^{-2}\mu)^{2}}{8\mu^{T}\Sigma^{-3}\mu},$

which is quadratic in $T$, which is the same behavior as the log-likelihood
ratio test.

##  VII Numerical Experiments

In this section, numerical experiments are conducted using various datasets in
order to numerically demonstrate the effectiveness of our large deviation
analysis of score-based hypothesis testing. Our experiments utilize both
synthetic and real-world data.

  1. 1.

Synethetic datasets: For synthetic datasets, we consider samples generated
from multivariate normal distributions, exponential family, and Gauss-
Bernoulli Restricted Boltzmann Machines (RBMs) [1]. When the exact likelihood
is available, as in the case of multivariate normal distributions and the
exponential family, we compare the performance of the Hyvarinen Score Test
(HST) with that of the Likelihood Ratio Test (LRT). For each distribution
type, we perform $10$ runs, with each run producing a datasets of size
$10,000$ created by sampling from both the null and alternative distributions.

  2. 2.

Real datasets: For real-world datasets, we utilize the KDD Cup’99 dataset
[20]. This dataset contains various types of simulated adversarial network
attacks and is typically used for bench-marking network intrusion detection
algorithms.

For all experiments, we evaluate the Positive Error Exponent (for Type I
error) and Negative Error Exponent (for Type II error) for various thresholds.
Using numerical calculations, we compare the estimates of the theoretical
results obtained from the large deviation analysis given above, with the
empirical error exponent obtained from hypothesis tests of size $n$ by
sampling with replacement (ranging from $n=1$ to $n=128$). The threshold for
calculating the theoretical limit is obtained empirically with $n=1$.

###  VII-A Details of Experimental Setup

For the multivariate normal distribution data, we consider the bivariate
normal distribution $\mathcal{N}(\mu,\Sigma)$, performing tests on $\mu$ with
known $\Sigma$, and vice versa. The null hypothesis $P_{\infty}$ is: The mean
$\mu_{0}$ (respectively the covariance $\Sigma_{0}$) is equal to its true
value $\mu^{\star}$ (or $\Sigma^{\star}$). The corresponding values of the
alternative hypothesis $P_{1}$ are assigned by adding i.i.d. samples of a
noise term $\mathcal{N}(0,\sigma_{ptb}^{2})$ to coordinates of $\mu_{\infty}$
(or by adding a i.i.d. samples of a log-Normal noise to diagonal elements of
$\Sigma_{\infty}$, where the logarithm of each noise term is distributed
according to a Normal distribution $\mathcal{N}(0,\sigma_{ptb}^{2})$). In our
simulations, we let $\mu_{\infty}=(0,0)^{T}$ (respectively
$\Sigma_{\infty}=\left[{\begin{array}[]{cc}1&0.7\\\ 0.7&1\\\
\end{array}}\right]$) for the null hypotheses in each trial, testing against
alternatives $\mu_{1}$ (respectively $\Sigma_{1}$) as described above with
$\sigma_{ptb}=0.01$.

For the exponential family dataset, we use the random variable
$\mathbf{x}\in\mathbb{R}^{d}$ distributed according to the un-normalized pdf

$\displaystyle
p_{\tau}(\mathbf{x})\propto\exp\left\\{-\tau\left(\sum_{i=1}^{d}x_{i}^{4}+\sum_{1\leq
i\leq d,i\leq j\leq d}x_{i}^{2}x_{j}^{2}\right)\right\\},$ (18)

where $\tau\in\mathcal{T}\subset\mathbb{R}^{+}$ is the model parameter. This
subfamily of the exponential family correspond to pairwise interaction
graphical models [21]. We first consider the hypothesis test
$\mathcal{H}_{\infty}:\tau=\tau_{\infty}=1$ versus
$\mathcal{H}_{1}:\tau=\tau_{1}=\tau_{0}+\tau_{ptb}$, with $\tau_{ptb}=0.01$.
Hamiltonian Monte Carlo (HMC) is used to generate samples from the un-
normalized density functions. In order to perform LRT, we compute the
normalizing constant by numerical integration.

We also consider the RBM [7] model which is defined based on an undirected bi-
partite graphical model consisting of hidden and visible variables. The Gauss-
Bernoulli RBM has binary-valued hidden variables
$\mathbf{h}\in\\{0,1\\}^{d_{h}}$ and real-valued visible variables
$\mathbf{x}\in R^{d_{x}}$ with joint distribution

$\displaystyle p(\mathbf{x},\mathbf{h})=$
$\displaystyle\frac{1}{Z_{\theta}}\text{exp}\left\\{-\left(\frac{1}{2}\sum_{i=1}^{d_{x}}\sum_{j=1}^{d_{h}}\frac{x_{i}}{\sigma_{i}}W_{ij}h_{j}+\sum_{i=1}^{d_{x}}b_{i}x_{i}+\sum_{j=1}^{d_{h}}c_{j}h_{j}-\frac{1}{2}\sum_{i=1}^{d_{x}}\frac{x_{i}^{2}}{\sigma_{i}^{2}}\right)\right\\},$
(19)

where $\theta=(\mathbf{W},\mathbf{b},\mathbf{c})$ are model parameters and
$Z_{\theta}$ is the normalizing constant. We set $\sigma_{i}=1$ for all
$i=1,\dots,d_{x}$ in our experiments. The probability of the visible variable
$\mathbf{x}$ written as
$p(\mathbf{x})=\sum_{h\in\\{0,1\\}^{d_{h}}}p(\mathbf{x},\mathbf{h})=\frac{1}{Z_{\theta}}\exp\\{-F_{\theta}(\mathbf{x})\\}$,
where $F_{\theta}(\mathbf{x})$ is the free energy given by

$\displaystyle F_{\theta}(\mathbf{x})=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{d_{x}}(x_{i}-b_{i})^{2}-\sum_{j=1}^{d_{h}}\operatorname{Softplus}\left(\sum_{i=1}^{d_{x}}W_{ij}x_{i}+b_{j}\right).$
(20)

The $\operatorname{Softplus}$ function is defined as
$\operatorname{Softplus}(t)\overset{\Delta}{=}\log(1+\exp(t))$ with a default
scale parameter $\beta=1$. The corresponding Hyvärinen score
$S_{\texttt{H}}(\mathbf{X}_{n},\theta)$ is given by

$\displaystyle
S_{\texttt{H}}(\mathbf{X}_{n},\theta)=\sum_{n=1}^{n}\sum_{i=1}^{d_{x}}\left[\frac{1}{2}\left(x_{in}-b_{i}+\sum_{j=1}^{d_{h}}W_{ij}\delta_{jn}\right)^{2}+\sum_{j=1}^{d_{h}}W_{ij}^{2}\delta_{jn}\left(1-\delta_{jn}\right)-1\right],$
(21)

where
$\delta_{jn}\overset{\Delta}{=}\operatorname{Sigmoid}(\sum_{i=1}^{d_{x}}W_{ij}x_{in}+b_{j})$.
The $\operatorname{Sigmoid}$ function is defined as
$\operatorname{Sigmoid}(t)\overset{\Delta}{=}(1+\exp(-t))^{-1}$.

We randomly draw the weight matrix
$\mathbf{W}_{0}\in\mathbb{R}^{d_{x}}\times\mathbb{R}^{d_{h}}$ from the
standard Normal distribution. In our experiments, we set the dimension of
visible variables to $d_{x}=50$ and hidden variables to $d_{h}=40$. The weight
matrix of the alternative hypothesis is constructed by perturbing elements of
$\mathbf{W}_{0}$ with i.i.d samples of Normal distribution
$\mathcal{N}(0,\sigma_{ptb}^{2})$ with $\sigma_{ptb}=0.01$. Samples of RBMs
are drawn using Gibbs sampling with $1000$ RBM iterations in order to ensure
convergence.

For the KDD Cup’99 dataset, the ‘normal’ traffic network is treated as the
null hypothesis, against various types of adversarial network attacks as
alternative hypotheses. A Gauss-Bernoulli RBM is trained with all available
training data of the ‘normal’ traffic network to estimate an unnormalized
version of the null distribution $P_{\infty}$. Hypothesis tests are conducted
against various alternative hypotheses, including ‘back’, ‘ipsweep’,
‘neptune’, and others. The ‘unknown’ adversarial network attack combines all
types of attacks with a dataset size of $\leq 100$. To provide a comprehensive
overview, the number of available data points of the network traffic in the
KDD Cup’99 dataset are depicted in Table I.

TABLE I: Statistics of various network traffic in KDD Cup’99 dataset.

normal | neptune | back | teardrop | satan | warezclient | ipsweep | smurf | portsweep | pod | nmap | unknown
---|---|---|---|---|---|---|---|---|---|---|---
87832 | 51820 | 968 | 918 | 906 | 893 | 651 | 641 | 416 | 206 | 158 | 177

We conducted two sets of experiments. In the first set, we explicitly use the
un-normalized null ($P_{\infty}$) and alternative ($P_{1}$) distributions from
true data distributions of synthetic data. The synthetic data experiments
include normal distributions, exponential family, and Gauss-Bernoulli RBM. In
the second set of simulations (the Gauss-Bernoulli RBM and KDD Cup’99
datasets), we fit the alternative distribution $P_{1}$ using training data of
size $N$, with cold-start from the null distribution. For the null
distribution, we use the explicit un-normalized data distribution for the
synthetic Gauss-Bernoulli RBM experiment. For the KDD Cup’99 dataset
experiment, we fit the unnormalized null distribution with a Gauss-Bernoulli
RBM using all the available training data of the ‘normal’ traffic network. In
order to demonstrate the robustness of our results to the sampling of
alternative training data, we produced $100$ alternative distributions, each
trained with $N$ alternative data samples. For each test, one of the
alternative distributions is randomly selected. Furthermore, we conducted
ablation studies on the training data size $N$ in order to demonstrate the
influence of the underfitting/overfitting of the alternative distribution on
our large deviation analysis results. This is important because in many
practical applications (such as detecting adversarial network attacks),
alternative data are much more scarce than the null data.

###  VII-B Experimental Results

In our first set of experiments (as described above), we evaluate the
performance of the Hyvarinen Score Test (HST) and Likelihood Ratio Test (LRT)
for explicit data distributions. These distributions include multivariate
normal distributions, the exponential family, and Gauss-Bernoulli Restricted
Boltzmann Machines (RBMs). Results from these distributions are presented in
Figures 1 to 2. These results are consistent with our theoretical analysis
given in Section II. Specifically, as the sample size $n$ in the composite
hypothesis tests is increased, both the empirical positive and negative error
exponents converge to their corresponding theoretical limits.

In the second set of experiments (as described above), we focus on scenarios
where the alternative distribution $P_{1}$ is fit using the training data.
These results are depicted in Figures 5 through 16. In these experiments, we
observed that the empirical positive and negative error exponents occasionally
exceed the theoretical limits as $n$ grows large. However, this trend
disappears as the number of samples $N$ (used to fit the alternative
distribution) increases. As $N$ increases (for instance from $10$ to $100$ in
Figure 5 and from $10$ to $All$ in Figure 6), the empirical results more
frequently adhere to the theoretical limits, aligning with our expectation.
Furthermore, we also observe that when $N$ increases, the empirical error
exponent becomes more convex. These results suggest that a more accurately
fitted alternative distribution facilitates the convergence of the empirical
error exponents towards the theoretical limits. Further research could
investigate the influence of the fitting accuracy of the alternative
distribution on the applicability and effectiveness of the large deviation
analysis.

Figure 1: Large deviation analysis of likelihood-based and sore-based
hypothesis testing for multivariate normal distribution with perturbation on
$\mu$ and $\sigma_{ptb}=0.01$. Figure 2: Large deviation analysis of
likelihood-based and sore-based hypothesis testing for multivariate normal
distribution with perturbation on $\sigma$ and $\sigma_{ptb}=0.01$. Figure 3:
Large deviation analysis of likelihood-based and sore-based hypothesis testing
for exponential family with perturbation on $\tau$ and $\sigma_{ptb}=0.01$.
Figure 4: Large deviation analysis of sore-based hypothesis testing for Gauss-
Bernoulli RBM with perturbation on $W$ and $\sigma_{ptb}=0.01$. Figure 5:
Large deviation analysis of sore-based hypothesis testing for Gauss-Bernoulli
RBM ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$ samples) with
perturbation on $W$ and $\sigma_{ptb}=0.01$. Figure 6: Large deviation
analysis of sore-based hypothesis testing for ‘back’ attack on KDD Cup’99
dataset ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$ samples). Figure 7:
Large deviation analysis of sore-based hypothesis testing for ‘ipsweep’ attack
on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$
samples). Figure 8: Large deviation analysis of sore-based hypothesis testing
for ‘neptune’ attack on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-
Bernoulli RBM and $N$ samples). Figure 9: Large deviation analysis of sore-
based hypothesis testing for ‘nmap’ attack on KDD Cup’99 dataset ($P_{1}$
fitted with Gauss-Bernoulli RBM and $N$ samples). Figure 10: Large deviation
analysis of sore-based hypothesis testing for ‘pod’ attack on KDD Cup’99
dataset ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$ samples). Figure 11:
Large deviation analysis of sore-based hypothesis testing for ‘portsweep’
attack on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$
samples). Figure 12: Large deviation analysis of sore-based hypothesis testing
for ‘satan’ attack on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-Bernoulli
RBM and $N$ samples). Figure 13: Large deviation analysis of sore-based
hypothesis testing for ‘smurf’ attack on KDD Cup’99 dataset ($P_{1}$ fitted
with Gauss-Bernoulli RBM and $N$ samples). Figure 14: Large deviation analysis
of sore-based hypothesis testing for ‘teardrop’ attack on KDD Cup’99 dataset
($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$ samples). Figure 15: Large
deviation analysis of sore-based hypothesis testing for ‘warezclient’ attack
on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-Bernoulli RBM and $N$
samples). Figure 16: Large deviation analysis of sore-based hypothesis testing
for ‘unknown’ attack on KDD Cup’99 dataset ($P_{1}$ fitted with Gauss-
Bernoulli RBM and $N$ samples).

##  VIII Conclusion

In this work, we analyzed the performance of score-based hypothesis testing
method [1]. We derived upper bounds on Type I and II error probabilities, and
proved that the exponents of our error bounds become precise in large sample
size regimes. We calculated these error exponents numerically for a variety of
scenarios of interest. Our experiments using both synthetic and real-world
data demonstrated that the empirical error exponents follow our theoretical
analysis. Future research may focus on the impact of lack of accuracy in
fitting the alternative distribution on the applicability and effectiveness of
large deviation analysis.

## References

  * Wu et al. [2022] S. Wu, E. Diao, K. Elkhalil, J. Ding, and V. Tarokh, “Score-based hypothesis testing for unnormalized models,” _IEEE Access_ , vol. 10, pp. 71 936–71 950, 2022.
  * Hyvärinen [2005] A. Hyvärinen, “Estimation of non-normalized statistical models by score matching.” _J. Mach. Learn. Res._ , vol. 6, no. 4, 2005.
  * Hyvärinen [2007] ——, “Some extensions of score matching,” _Comput. Stat. Data Anal._ , vol. 51, no. 5, pp. 2499–2512, 2007.
  * Song et al. [2020] Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Poole, “Score-based generative modeling through stochastic differential equations,” _arXiv preprint arXiv:2011.13456_ , 2020.
  * Vahdat et al. [2021] A. Vahdat, K. Kreis, and J. Kautz, “Score-based generative modeling in latent space,” _Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 34, pp. 11 287–11 302, 2021.
  * Koller and Friedman [2009] D. Koller and N. Friedman, _Probabilistic graphical models: principles and techniques_. The MIT Press, 2009.
  * LeCun et al. [2006] Y. LeCun, S. Chopra, R. Hadsell, M. Ranzato, and F. Huang, “A tutorial on energy-based learning,” in _Predicting structured data_. The MIT Press, 2006, vol. 1.
  * Papamakarios et al. [2021] G. Papamakarios, E. T. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, “Normalizing flows for probabilistic modeling and inference,” _J. Mach. Learn. Res._ , vol. 22, pp. 57:1–57:64, 2021.
  * Cerf and Petit [2011] R. Cerf and P. Petit, “A short proof of cram´er’s theorem in $\mathbb{R}$,” _The American Mathematical Monthly_ , vol. 118, no. 10, pp. 925–931, 2011.
  * Dembo and Zeitouni [2009] A. Dembo and O. Zeitouni, _Large deviations techniques and applications_. Springer Science & Business Media, 2009, vol. 38.
  * Moulin and Veeravalli [2018] P. Moulin and V. V. Veeravalli, _Statistical Inference for Engineers and Data Scientists_. Cambridge University Press, 2018.
  * Lehmann et al. [1986] E. L. Lehmann, J. P. Romano, and G. Casella, _Testing statistical hypotheses_. Springer, 1986, vol. 3.
  * Wasserman [2006] L. Wasserman, _All of nonparametric statistics_. Springer Science & Business Media, 2006.
  * Wu et al. [2023a] S. Wu, E. Diao, T. Banerjee, J. Ding, and V. Tarokh, “Score-based change point detection for unnormalized models,” _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , 2023.
  * Song and Ermon [2019] Y. Song and S. Ermon, “Generative modeling by estimating gradients of the data distribution,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * Vincent [2011] P. Vincent, “A connection between score matching and denoising autoencoders,” _Neural computation_ , vol. 23, no. 7, pp. 1661–1674, 2011.
  * Wu et al. [2023b] S. Wu, E. Diao, T. Banerjee, J. Ding, and V. Tarokh, “Robust quickest change detection for unnormalized models,” _Conference on Uncertainty in Artificial Intelligence (UAI)_ , 2023.
  * Mahmood et al. [2020] A. Mahmood, J. Oliva, and M. A. Styner, “Multiscale score matching for out-of-distribution detection,” in _International Conference on Learning Representations_ , 2020.
  * Kulinski et al. [2020] S. Kulinski, S. Bagchi, and D. I. Inouye, “Feature shift detection: Localizing which features have shifted via conditional distribution tests,” _Advances in neural information processing systems_ , vol. 33, pp. 19 523–19 533, 2020.
  * Lippmann et al. [2000] R. Lippmann, J. W. Haines, D. J. Fried, J. Korba, and K. Das, “Analysis and results of the 1999 darpa off-line intrusion detection evaluation,” in _International Workshop on Recent Advances in Intrusion Detection_. Springer, 2000, pp. 162–182.
  * Yu et al. [2016] M. Yu, M. Kolar, and V. Gupta, “Statistical inference for pairwise graphical models using score matching,” _Advances in Neural Information Processing Systems (NeurIPS)_ , vol. 29, 2016.