mmd short=MMD, long=maximum mean discrepancy, long-plural-form=maximum mean
discrepancies, ksd short=KSD, long=kernel Stein discrepancy, long-plural-
form=kernel Stein discrepancies, ipm short=IPM, long=integral probability
metric, rkhs short=RKHS, long=reproducing kernel Hilbert space, pd short=PD,
long=positive definite, dtmc short=DTMC, long=discrete time Markov chain, ctmc
short=CTMC, long=continuous time Markov chain, mc short=MC, long=Markov chain,
lr short=LR, long=likelihood ratio, gof short=GoF, long=goodness-of-fit,

A Kernel Stein Test of Goodness of Fit for Sequential Models

Jerome Baum Heishiro Kanagawa Arthur Gretton

University College London Gatsby Computational Neuroscience Unit University
College London

###### Abstract

We propose a goodness-of-fit measure for probability densities modelling
observations with varying dimensionality, such as text documents of differing
lengths or variable-length sequences. The proposed measure is an instance of
the kernel Stein discrepancy (KSD), which has been used to construct goodness-
of-fit tests for unnormalised densities. Existing KSDs require the model to be
defined on a fixed-dimension space. As our major contributions, we extend the
KSD to the variable dimension setting by identifying appropriate Stein
operators, and propose a novel KSD goodness-of-fit test. As with the previous
variants, the proposed KSD does not require the density to be normalised,
allowing the evaluation of a large class of models. Our test is shown to
perform well in practice on discrete sequential data benchmarks.

##  1 INTRODUCTION

This paper addresses the problem of evaluating a probabilistic model for
observations with variable dimensionality. This problem commonly arises in
sequential data analysis, where sequences of different lengths are observed,
and generative models (e.g., Markov models) are used to draw inferences. Our
problem setting also concerns a scenario where a data point is a collection of
observations that may not be sequentially ordered, as in the topic modelling
of text documents. Our task is formalized as follows: given a sample
$\\{X_{i}\\}_{i=1}^{n}$ from a distribution $Q$ on a sample space
$\mathcal{X}$ (e.g., the set of all possible sequences), we aim to quantify
the discrepancy of our model distribution $P$ with respect to $Q$ .

One approach to this problem is generating samples from the model and
computing a sample-based discrepancy measure, such as the maximum mean
discrepancy [11, 13]. A disadvantage of this approach is that it potentially
requires many samples to see the departure of the model from the data. For
example, if a sequence model misspecifies a particular transition from a given
history, we would need to generate sequences sharing this history to observe
the disparity. Unfortunately, observing such sequences becomes more
challenging as the state space size increases or the history length gets
longer, necessitating repeated sampling. This observation motivates using a
measure that exploits information provided by the model, such as dependence
relations between different time points.

A model-dependent measure may be derived based on Stein’s method [4], a
technique from probability theory developed originally to obtain explicit
rates of convergence to normality. The key construct from Stein’s method is a
distribution-specific operator called a Stein operator that modifies a
function to have zero expectation under the distribution. A model-specific
Stein operator $\mathcal{A}_{P}$ may be defined to construct a zero-
expectation function $\mathcal{A}_{P}f$; its expectation $\mathbb{E}_{X\sim
Q}[\mathcal{A}_{P}f(X)]$ under the sample distribution serves as a discrepancy
measure, since a non-zero expectation indicates the deviation from the model.
One may generalize this idea to a family of functions $\mathcal{F}$, and the
resulting worst-case summary $\sup_{f\in\mathcal{F}}\lvert\mathbb{E}_{X\sim
Q}[\mathcal{A}_{P}f(X)]\rvert$ is called a Stein discrepancy, as proposed by
[12]. An appropriate choice of the operator and the function class yields a
computable Stein discrepancy. This paper focuses on an instance called the
kernel Stein discrepancy (KSD) [16, 14, 15, 17], where functions from a
reproducing kernel Hilbert space (RKHS) are used.

The KSD is a versatile framework for designing goodness-of-fit measures. Given
a Stein operator and a reproducing kernel, the KSD is expressed by an
expectation of a Stein-modified kernel, leading to tractable estimators only
involving sample evaluations of the modified kernel. This feature allows us to
use a wealth of kernels from the literature: by designing a Stein operator,
one can adapt the kernel to the model and define a KSD appropriate for the
problem. Based on the Stein kernel theory by [16], [14] and [15] originated
this line of work, where they combined the Langevin Stein operator proposed by
[12] with an RKHS; remarkably, the resulting (Langevin) KSD is computable for
densities with unknown normalizing constants. There have been numerous
extensions to models in other data domains: categorical data [18], point-
pattern data [24], censored data [26], directional data [28], and functional
data [31]. We refer the reader to [29] for a review on the KSD’s applications
out goodness-of-fit testing.

A limitation of the preceding Stein works is that they require the model
distribution to be defined on a fixed-dimension space. For example, [21]
evaluate latent Dirichlet allocation models [7] for text documents assuming a
fixed document length, an assumption unlikely to hold in practice. Relevant
work has been accomplished by [31], where they propose a KSD for distributions
on an (infinite-dimensional) Hilbert space, and treat continuous-time models,
such as Brownian motions and stochastic differential equation models. While
their test applies to sequential models, our target setting is different, as
we focus on discrete-time models such as Markov chains and do not require a
density with respect to a Gaussian measure.

In the present work, we extend the KSD framework to the variable-dimension
setting. Specifically, we identify a Stein operator for distributions on the
set $\mathcal{X}=\bigsqcup_{\ell=1}^{\ell_{\mathrm{max}}}\mathcal{S}^{(l)}$ of
all univariate sequences of finite alphabets; i.e., a disjoint union of finite
sets $\mathcal{S}^{(l)}$. Our Stein operator builds on the Zanella Stein
operator recently proposed by [27], which is derived from a Markov jump
process having the target distribution as invariant measure; we review this in
Section 2. Next, in Section 3, we derive a Stein operator using a Markov
process admitting transitions between sets of different sizes. Based on the
new Stein operator and the associated KSD, we propose a novel goodness-of-fit
test for sequential models. As in previous variants, the proposed KSD does not
require the model to be normalized, allowing the evaluation of models such as
Markov random fields, as demonstrated in Section 4. As the proposed operator
involves a number of tunable parameters, as our second contribution, we offer
guidance on how we select them based on empirical power analysis.

##  2 BACKGROUND

Given an independent identically distributed sample $\\{X_{i}\\}_{i=1}^{n}$
from a distribution $Q$ on a discrete set $\mathcal{X}$, we wish to test the
goodness of fit of a model $P$. In this section, we recall the definition of
the Zanella Stein operator and the associated kernel Stein discrepancy (KSD),
a goodness-of-fit measure used by our test. In the following, we assume that
the model’s probability mass function $p(x)$ is positive everywhere in
$\mathcal{X}$.

##### Zanella Stein operator

The Zanella Stein operator is a Stein operator for distributions on a
countable set $\mathcal{X}$, proposed by [27] following the generator method
of [5]. For a real-valued function on $\mathcal{X}$, the Zanella Stein
operator is defined by

$\displaystyle\left(\mathcal{A}_{P}f\right)(x)$
$\displaystyle\coloneqq\sum_{y\in\partial
x}g\left(\frac{p(y)}{p(x)}\right)\left(f(y)-f(x)\right),$ (1)

where $\partial x\subset\mathcal{X}$ is a set of _neighborhood_ points, and
$g:(0,\infty)\to(0,\infty)$ is a _balancing function_ , which must satisfy
$g\left(t\right)=tg\left(1/t\right)$ for $t>0$. Several choices of $g$ are
known: the minimum probability flow operator discussed in [19] uses
$g\left(t\right)=\sqrt{t}$; the choice $g\left(t\right)=t/(1+t)$ is known as
the Barker Stein operator [3, 30].

The Zanella Stein operator is the infinitesimal generator of a Markov jump
process called the Zanella process [25, 22] designed to have the target
distribution $P$ as an invariant distribution. The process is based on the
idea of _locally informed proposals_ : these jump from a given point $x$ to
any of its neighbors $y$, at a rate $g\left(p(y)/p(x)\right)$ so that detailed
balance is satisfied for $P$, ensuring that $P$ is an invariant distribution.
The Zanella process is not the only allowable choice: other Markov chains or
processes have also been considered, including birth-death processes [30], the
Glauber dynamics Markov chain [23, 20]. We refer the reader to [30] for a
review. The Zanella process is attractive since it applies to any discrete set
given an appropriate neighborhood structure. We make use of this property to
construct our discrepancy measure in the next section.

For the operator to characterize $P$, the neighborhood must be chosen
accordingly. Let $(\mathcal{V},\mathcal{E})$ be the directed graph with
vertices $\mathcal{V}=\mathcal{X}$ and edges
$\mathcal{E}=\\{(x,y)|y\in\partial x\\}$ induced by the chosen neighborhood
structure. The graph $(\mathcal{V},\mathcal{E})$ must have the following
properties:

  1. 1.

Symmetry: $(x,y)\in\mathcal{E}\Leftrightarrow(y,x)\in\mathcal{E}$.

  2. 2.

Aperiodicity: there is no integer $j>1$ that divides the length of every cycle
of the graph.

  3. 3.

Strong connectivity : the transitive closure of $\mathcal{E}$ is
$\mathcal{X}\times\mathcal{X}$, so that there is a path from every point to
every other point.

In fact, symmetry alone implies that $P$ is in the set of invariant
distributions; the other two properties are required to ensure that $P$ is the
unique invariant distribution, as otherwise the corresponding Stein
discrepancy cannot distinguish $P$ from other distributions even for a rich
test function class. Note that we can conveniently choose a sparse
neighborhood to accelerate the computation, although this may result in the
reduced sensitivity of the discrepancy. In Section 4, we discuss in more
detail the tradeoff between computational cost and test power.

##### Kernel Stein Discrepancy

A computable discrepancy measure may be defined based on the Stein operator
and a reproducing kernel Hilbert space (RKHS) [2, 8]. Following [16, 14, 15,
17], we define the kernel Stein discrepancy (KSD) associated with the Zanella
Stein operator as follows:

$\displaystyle\operatorname{KSD}\bigl{[}Q\|P\bigr{]}\coloneqq\sup_{\lVert
f\rVert_{\mathcal{H}}\leq 1}\bigl{|}\mathbb{E}_{X\sim
Q}[\mathcal{A}_{P}(f)(X)]\bigr{|},$ (2)

where $\mathcal{H}$ is the RKHS of real-valued functions on $\mathcal{X}$
corresponding to a positive definite kernel
$k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ with the natural norm $\lVert
f\rVert_{\mathcal{H}}=\sqrt{\langle f,f\rangle_{\mathcal{H}}}$ defined by the
inner product $\langle\cdot,\cdot\rangle_{\mathcal{H}}$. By the vanishing
property $\mathbb{E}_{Y\sim P}[\mathcal{A}_{P}f(Y)]=0$, the KSD is zero if
$Q=P$. The other direction requires further assumptions on the operator and
the RKHS. According to [27, Proposition 4], when the corresponding Zanella
process is exponentially ergodic and the RKHS is $C_{0}$-universal, we have
$\operatorname{KSD}\bigl{[}Q\|P\bigr{]}=0$ only if $P=Q$. For a finite state
space $\mathcal{X}$, exponential ergodicity is satisfied by any irredicible
Markov jump process. The $C_{0}$-universality is equivalent to the integrally
strict positive definiteness [10]; in the particular case
$\lvert\mathcal{X}\rvert<\infty$, this condition states that the Gram matrix
defined over all the configurations on $\mathcal{X}$ is strictly positive
definite.

The use of an RKHS allows us to obtain an explicit solution to the
optimisation problem in the definition, and yields a closed-form expression
[27],

$\displaystyle\operatorname{KSD}^{2}\bigl{[}Q\|P\bigr{]}=\mathbb{E}_{X,X^{\prime}\sim
Q\otimes Q}[h_{p}(X,X^{\prime})],$ (3)

where $X$, $X^{\prime}$ are i.i.d. random variables with law $Q$, provided
$\mathbb{E}_{X\sim Q}[h_{p}(X,X)^{1/2}]<\infty$. The proof of this result
follows as in [17], noting that the inner product between $f$ and
$\mathcal{A}_{P}k(x,\cdot)$ _reproduces_ $\mathcal{A}_{P}f(x)$ for any
$f\in\mathcal{H}$. The function $h_{p}$ is called a _Stein kernel_ , defined
by

$\displaystyle\begin{aligned} &h_{p}(x,y)\\\ &=\sum_{\nu\in{\cal
N}_{x}}\sum_{\tilde{\nu}\in{\cal
N}_{y}}g_{\nu}(x)g_{\tilde{\nu}}(y)\left\\{k(\nu(x),\tilde{\nu}(y))+k(x,y)\right.\\\
&\hphantom{=\sum_{\nu\in{\cal N}_{x}}\sum_{\tilde{\nu}\in{\cal
N}_{y}}g_{\nu}(x)g_{\tilde{\nu}}(y)}\left.\quad-k(x,\tilde{\nu}(y))-k(\nu(x),y)\right\\}\end{aligned},$

where we identify each point in the neighborhood $\partial x$ as a mapping
from $x$ to itself, and $g_{\nu}(x)=g\\{p(\nu(x))/p(x)\\}$.

Given a sample $\\{X_{i}\\}_{i=1}^{n}\overset{\text{i.i.d.}}{\sim}Q$, the
squared KSD expression (3) admits a simple unbiased estimator

$U_{n,P}\left(\\{X_{i}\\}_{i=1}^{n}\right)\coloneqq\frac{1}{n(n-1)}\sum_{i\neq
j}^{n}h_{p}(X_{i},X_{j}),$ (4)

which is a U-statistic [1]. The U-statistic is degenerate if $Q=P$ [14, 15].
In this case, the scaled statistic
$nU_{n,P}\left(\\{X_{i}\\}_{i=1}^{n}\right)$ asymptotically follows the law of
$\sum_{j=1}^{\infty}\lambda_{j}(Z_{j}^{2}-1)$, where $Z_{j}$ are i.i.d.
standard normal variables, and $\lambda_{j}$ are eigenvalues given by the
eigenvalue problem
$\sum_{x^{\prime}\in\mathcal{X}}h_{p}(x,x^{\prime})a(x^{\prime})p(x^{\prime})=\lambda
a(x)$ with $\sum_{x\in\mathcal{X}}a(x)^{2}p(x)<\infty$ [9, see, e.g.,]. One of
our proposed tests in Section 3 simulates this distribution using a wild
bootstrap procedure.

##  3 TESTING SEQUENTIAL MODELS

We now address the variable-dimension setting. We begin by formally defining
the sample space $\mathcal{X}$. Let $\mathcal{S}$ be a nonempty finite set of
symbols. For an integer $\ell\geq 1$, we denote by $\mathcal{S}^{(\ell)}$ the
set of all length-$\ell$ sequences formed by symbols in $\mathcal{S}$; i.e.,
$\mathcal{S}^{(\ell)}=\prod_{j=1}^{\ell}\mathcal{S}$. Then, the sample space
$\mathcal{X}$ is given by the set of all sequences
$\bigsqcup_{\ell=1}^{\ell_{\mathrm{max}}}\mathcal{S}^{(\ell)}$, where
$\bigsqcup$ denotes the disjoint union, and $\ell_{\mathrm{max}}$ is $\infty$
or a finite positive integer. In this setting, a sample $X$ from a
distribution on $\mathcal{X}$ corresponds to a sequence whose length is
randomly determined.

Algorithm 1 Parametric bootstrap test

1:$P$; a target distribution on $\mathcal{X}$.

2:$D_{n}=\\{X_{1},\dots,X_{n}\\}$; data.

3:$\alpha$; the desired level of the test.

4:for $b=1,\ldots,B$ do

5: Sample
$D_{n}^{(b)}=\\{X^{(b)}_{1},\dots,X^{(b)}_{n}\\}\overset{\text{i.i.d.}}{\sim}P$

6: $\Delta_{b}\leftarrow U_{n,P}\left(D_{n}^{(b)}\right)$ using Equation 4

7:end for

8:$\hat{t}_{\alpha}\leftarrow(1-\alpha)$-quantile of
$\\{\Delta_{1},\dots,\Delta_{B}\\}$

9:if $U_{n,P}\left(D_{n}\right)>\hat{t}_{\alpha}$ then

10: Reject $H_{0}$

11:end if

Algorithm 2 Wild bootstrap test

1:$P$; a target distribution on $\mathcal{X}$.

2:$D_{n}=\\{X_{1},\dots,X_{n}\\}$; data.

3:$\alpha$; the desired level of the test.

4:for $b=1,\ldots,B$ do

5: Draw $W_{b}\sim\mathrm{Multi}(n;1/n,\dots,1/n)$

6: $\Delta_{b}\leftarrow U_{n,P}^{W_{b}}\left(D_{n}\right)$ using Equation 6

7:end for

8:$\hat{t}_{\alpha}\leftarrow(1-\alpha)$-quantile of
$\\{\Delta_{1},\dots,\Delta_{B}\\}$

9:if $U_{n,P}\left(D_{n}\right)>\hat{t}_{\alpha}$ then

10: Reject $H_{0}$

11:end if

##### Neighborhood Choice

To obtain a KSD in this setting, we need to specify a neighborhood structure
for the the Zanella Stein operator. We first introduce a simple structure for
sequences, which edits only the end of the string, to illustrate the required
connectivity principles. We then describe a a more advanced neighbourhood
structure, allowing for greater modifications, which will be the foundation of
our tests.

As noted in Section 2, for the Stein operator to uniquely characterise the
model, we require the strong connectivity of the induced graph. We achieve
this by introducing _inter- and intra-length state transitions_. Specifically,
for a length-$\ell$ sequence
$x_{(1:\ell)}=(x_{1},\dots,x_{\ell})\in\mathcal{S}^{(l)}$, we propose the
following neighborhood:

$\displaystyle\partial
x_{(1:\ell)}=\mathcal{I}_{x}\cup\\{x_{(1:l-1)}\\}\cup\mathcal{R}_{x}$ (5)

with

$\displaystyle\mathcal{I}_{x}=\\{(x_{1},\dots,x_{l},s):s\in\mathcal{S}\\},$
$\displaystyle\mathcal{R}_{x}=\\{(x_{1},\dots,s):s\in\mathcal{N}_{\mathrm{rep}}(x_{\ell})\\},$

where $\mathcal{N}_{\mathrm{rep}}(s)$ denotes a neighborhood chosen for a
symbol $s\in\mathcal{S}$; this is chosen such that the induced graph is
strongly connected. The first two sets in (5) represent the _inter-length
trantisions_ : Each point of $\mathcal{I}_{x}$ corresponds to adding a symbol
to the end of the given sequence $x_{(1:\ell)}$, whereas $\\{x_{(1:l-1)}\\}$
stands for the deletion of the last symbol and shortening the sequence. The
third set $\mathcal{R}_{x}$ represents _intra-length transitions_ , where the
last symbol of the sequence is replaced. A trivial choice for satisfying the
strong connectivity is using the whole alphabet set $\mathcal{S}$ for
$\mathcal{N}_{\mathrm{rep}}(s)$ for each $s\in\mathcal{S}$. This approach may
be permitted when the alphabet size is small. One might otherwise want to
consider using smaller sets to reduce the computational cost. For example, we
can use a smaller $\mathcal{N}_{\mathrm{rep}}(s)$, provided that the graph
(with vertices $\mathcal{S}$) induced by the choice satisfies the three
conditions required in Section 2; e.g., as in [18], we may introduce a cyclic
order $\\{0,\dots,\lvert\mathcal{S}\rvert-1\\}$ to $\mathcal{S}$, and take two
adjacent points $\mathcal{N}_{\mathrm{rep}}(s)=\\{s+1,s-1\\}$, where $s\pm 1$
denotes $(s\pm 1)\ \mathrm{mod}\ \lvert\mathcal{S}\rvert-1$.

The structure proposed above is a minimal neighborhood choice that guarantees
strong connectivity. As we will see in Section 4, however, this minimal choice
may not produce a powerful test, as it only uses information of neighboring
length sequences. For example, the tail insertion operation by
$\mathcal{I}_{x}$ only encodes the structure of the target distribution in
terms of the end of the sequence. For sequential models such as Markov chain
models, the longer a sequence is, the more evidence we should have in favor
of, or against, a given model. Editing only the final element of a sequence
would therefore not make use of such accumulated evidence. Indeed, this
approach performs relatively weakly in our benchmark experiments (Figure 2(b),
case of $J=1$).

##### Proposed Neighborhood

Based on the above observation, we generalize the above neighborhood set.
Specifically, we consider expanding the neighborhood set by modifying sequence
elements at different locations: we propose the $J$-location modification
neighborhood

$\partial_{x_{(1:\ell)}}^{J}=\cup_{j=1}^{J}\mathcal{I}_{x_{(1:\ell)},j}\cup\mathcal{D}_{x_{(1:\ell)},j}\cup\mathcal{R}_{x_{(1:\ell)},j},$

where

$\displaystyle\mathcal{I}_{x_{(1:\ell)},j}\subset\mathcal{S}^{(\ell+1)}=\\{\operatorname{ins}(x_{(1:\ell)},j,s):s\in{\mathcal{N}}_{\mathrm{\operatorname{ins}}}\\},$
$\displaystyle\mathcal{D}_{x_{(1:\ell)},j}\subset\mathcal{S}^{(\ell-1)}=\\{\operatorname{del}(x_{(1:\ell)},j)\text{
if }x_{\ell-j+1}\in\mathcal{N}_{\operatorname{ins}}\\},$
$\displaystyle\mathcal{R}_{x_{(1:\ell)},j}\subset\mathcal{S}^{(\ell)}=\\{\operatorname{rep}(x_{(1:\ell)},j,s):s\in\mathcal{N}_{\operatorname{rep}}(x_{\ell})\\}.$

Here, $\operatorname{ins}(x_{(1:\ell)},j,s)$ denotes the sequence extending
$x_{(1:\ell)}$ by inserting $s$ at the $(j-1)$-th location counting back from
the end of the string, using symbols from a fixed set
$\mathcal{N}_{\operatorname{ins}}\subset\mathcal{S}$;
$\operatorname{del}(x_{(1:\ell)},j)$ deletes the $(j-1)$-th element counting
back from the end of $x_{(1:\ell)}$; and
$\operatorname{rep}(x_{(1:\ell)},j,s)$ replaces by $s$ the $(j-1)$-th element
from the end of $x_{(1:\ell)}$. Note that the deletion and insertion
operations must be paired to ensure the symmetry of the graph.

##### Test Procedure

Using one of the proposed neighborhoods for the Zanella Stein operator, we can
define a KSD goodness-of-fit test for sequential models. We consider testing
the null $H_{0}:\operatorname{KSD}\bigl{[}Q\|P\bigr{]}=0$ against the
alternative $H_{1}:\operatorname{KSD}\bigl{[}Q\|P\bigr{]}\neq 0$. Note that
when the KSD distinguishes any distributions, the null and alternative become
$H_{0}:P=Q$ and $H_{1}:P\neq Q$, respectively. Below, we present two tests,
which differ in how they compute test thresholds.

The first test is the parametric bootstrap test in Algorithm 1. The parametric
bootstrap repeatedly draws samples from the model $P$ and simulates the
distribution of the statistic (4). Because of this feature, this procedure
does not apply if sampling from the model is infeasible.

The second test is the wild bootstrap test described in Algorithm 2, which is
analogous to the existing KSD tests [14, 15, 18]. As we have seen in Section
2, the asymptotic distribution of the (scaled) KSD estimate (4) is known. The
wild bootstrap test simulates this distribution by repeatedly computing the
following quantity:

$U_{n,P}^{W}\left(\\{X_{i}\\}_{i=1}^{n}\right)\coloneqq\frac{\sum_{i\neq
j}^{n}(W_{i}-1)(W_{j}-1)h_{p}(X_{i},X_{j})}{n(n-1)}\ $ (6)

where $W=(W_{1},\dots,W_{n})$ is a sample from a multinomial distribution
$\mathrm{Multi}(n;1/n,\dots,1/n)$ with $n$ trials and $n$ events and
independent of the data $\\{X_{1},\dots,X_{n}\\}$. The test has an only
asymptotically correct size [6], and may not be suitable for a small sample
size $n$. In contrast to the parametric bootstrap test, this test does not
require sampling from the model.

##  4 EXPERIMENTS

We investigate the performance of the proposed test through synthetic
experiments.111Code is available at https:/github.com/test-for-sequential-
models/code. In particular, we aim to answer the following questions: (a) how
leveraging model structure using a Stein kernel improves test performance over
the maximum mean discrepancy (MMD) [11]; (b) the effect of neighborhood choice
on the KSD test power.

Throughout our experiments, we will make use of two kernels. The first is the
exponentiated Hamming kernel,
$k_{\mathrm{H}}(x_{(1:\ell)},y_{(1:\ell)})\coloneqq\exp\\{-\frac{1}{\ell}\sum_{i=1}^{\ell}\delta_{x_{i}}\left(y_{i}\right)\\}$,
where $\delta_{x}\left(y\right)=1$ if $x=y$ and $0$ otherwise. Additionally,
we define $k_{\mathrm{H}}(x_{(1:\ell)},y_{(1:\ell^{\prime})})=0$ whenever
$\ell\neq\ell^{\prime}$. The second is a string kernel which we refer to as
the _contiguous subsequence kernel_ (CSK). This kernel counts the number of
contiguous subsequences of a particular length, present in the two input
sequences,

$\displaystyle k_{\mathrm{CSK}}\left(x,y\right)$
$\displaystyle\coloneqq\frac{k_{\text{u}}(x,y)}{\sqrt{k_{\text{u}}(x,x)\cdot
k_{\text{u}}(y,y)}},$

where

$\displaystyle k_{\text{u}}(x_{(1:\ell)},y_{(1:\ell^{\prime})})$
$\displaystyle\coloneqq\sum_{i=1}^{\ell-t+1}\sum_{j=1}^{\ell^{\prime}-t+1}\delta_{x_{(i:i+t-1)}}\left(y_{(j:j+t-1)}\right),$

with parameter $t$ for the subsequence length. We use the normalized version
so that the kernel does not overly depend on the sequence lengths. The Hamming
kernel does not take into account the sequence structure and can be considered
as a naive choice, whereas the CSK kernel does.

###  4.1 Demonstration of Test Power

We first show that the KSD indeed benefits from information supplied by the
model. For this purpose, we construct a problem where the alternative is given
by adding a slight perturbation to the model. We choose
$\mathcal{S}=\\{0,\ldots,n-1\\}$ for our alphabet, with $n=8$. The model
distribution is a random walk through $\mathcal{S}$, starting from a random
uniform element of $\mathcal{S}$, with increments drawn uniformly from
$\\{-1,+1\\}$. We use a cyclic structure as described earlier, which
identifies $0$ with $n$ such that the random walk wraps around. The chain
terminates with probability $1/20$ after each step. We construct an
alternative distribution as follows: At every step, if and only the current
element is $0\in\mathcal{S}$, we randomly replace the increment with $0$ with
a probability $p_{\text{hold}}\in[0,1/4]$. That is, we introduce a random
holding step that occurs on average in $1/(4n)=1/32$ of steps, when the
perturbation parameter is at its maximum. The perturbation is unlikely to be
observed in any given sample, making the problem challenging.

For the KSD test we use the Zanella-Stein operator with the $J$-location
modification neighbourhood $\partial^{J}_{x_{(1:\ell)}}$, with $J=\infty$;
i.e., we allow edits anywhere in the sequence. As baselines we use an MMD test
that draws 100 samples from the model, a likelihood-ratio test where the
alternative consists of all first-order Markov chains on $\mathcal{S}$, and an
oracle likelihood-ratio test where the alternative is the (generally unknown)
ground truth distribution. Both kernel tests are built using the CSK kernel.
We use parametric bootstrap for all four tests. Tests are evaluated on $n=30$
i.i.d. sample points from the perturbed distribution.

$0$$0.1$$0.2$$0$$0.2$$0.4$$0.6$$0.8$$1$Holding probability in the ground
truth, $p_{\text{hold}}$Estimated rejection rateKSDMMDLR oracleLR MCLevel
$\alpha=0.05$ Figure 1: Estimated rejection rate (test power) of the proposed
test when the model distribution and ground truth differ by a perturbation
that is unlikely to be observed in any given sample. The model distribution is
a random walk through $\mathbb{Z}/n\mathbb{Z}$ with $n=8$. The ground truth
introduces a random holding period in the state $0$. The baselines are an MMD
test that samples from the model distribution, and a likelihood-ratio test
where the alternative consists of first-order Markov chains.

In Figure 1 we can see that the KSD test outperforms the two non-oracle
baselines. At low degrees of perturbation, the problem is challenging because
perturbations are rarely observed. The KSD test is able to overcome this
difficulty by exploiting information provided by the model.

###  4.2 Choice of Neighbourhood Graph

The test that we chose for the preceding experiment performed well, in part
due to a good choice of neighbourhood. We now further investigate the choice
of neighbourhood for the Zanella-Stein operator. In order to compare the
performance of different choices of Stein kernel, we constructed 12 synthetic
testing scenarios with various properties, ensuring our conclusions are
broadly applicable, and not specific to properties of a particular model. The
experiments in this section will evaluate the average performance of each test
across these scenarios. We use discrete-time Markov chains for the model $P$
and ground truth $Q$. The testing scenarios are fully specified in the
Appendix.

More densely connected neighbourhood graphs should give more powerful tests,
but such tests are slower to evaluate. To understand the tradeoff between
compute budget and test power, we construct a family of KSD tests
parameterized by the size of the neighbourhood of each point. We consider
three classes of neighbourhood: The $J$-location modification neighbourhood
$\partial^{J}_{x_{(1:\ell)}}$; the substitution neighbourhood
$\cup_{j=1}^{J}\mathcal{R}_{x_{(1:\ell)},j}$; and the neighbourhood of
insertions and deletions
$\cup_{j=1}^{J}\mathcal{I}_{x_{(1:\ell)},j}\cup\mathcal{D}_{x_{(1:\ell)},j}$.
In each case, the symbol neighbourhoods $\mathcal{N}_{\text{ins}}$ and
$\mathcal{N}_{\operatorname{rep}}$ are set to the entire alphabet
$\mathcal{S}$. We construct a Zanella-Stein operator using the Barker
balancing function $g\left(t\right)=t/(1+t)$. For the resulting operators, we
write $\mathrm{ZS}$, $\mathrm{ZS}^{\prime}$, and $\mathrm{ZS}^{\prime\prime}$,
respectively. Using these Stein operators, we construct a KSD by applying the
operators on the two kernels specified earlier: the Hamming kernel and the CSK
kernel. To estimate the critical value, we use parametric bootstrap. We also
include a pair of MMD tests, one for each kernel, as a baseline.

$1$$3$$5$$7$$9$$0$$0.2$$0.4$$0.6$$\cdots$$\infty$Neighbourhood size parameter
$J$Estimated rejection rateHamming $\mathrm{ZS}$Hamming
$\mathrm{ZS}^{\prime}$Hamming $\mathrm{ZS}^{\prime\prime}$Hamming MMDCSK
$\mathrm{ZS}$CSK $\mathrm{ZS}^{\prime}$CSK $\mathrm{ZS}^{\prime\prime}$CSK
MMDLevel $\alpha=0.05$ (a)

$0$$5$$10$$15$$20$$0$$0.2$$0.4$$0.6$Median time in secondsEstimated rejection
rateHamming $\mathrm{ZS}$Hamming $\mathrm{ZS}^{\prime}$Hamming
$\mathrm{ZS}^{\prime\prime}$CSK $\mathrm{ZS}$CSK $\mathrm{ZS}^{\prime}$CSK
$\mathrm{ZS}^{\prime\prime}$Level $\alpha=0.05$ (b)

Figure 2: Estimated rejection rate of three different families of Zanella-
Stein test with varying neighbourhood size and underlying kernel. We evaluate
the tests on 12 different synthetic testing scenarios, evaluating each test
multiple times on independently sampled synthetic datasets. To estimate the
overall rejection rate, we report a simple average across all test
evaluations. Figure (a) shows the estimated rejection rate as a function of
the neighbourhood size parameter $J$. Figure (b) shows the same as a function
of the median time taken to bootstrap and evaluate the test once, in seconds.
Test power increases as we increase the size of the neighbourhood. We see
diminishing returns particularly when the Hamming kernel is used to construct
the test.

Results are in Figure 2(a). Larger neighbourhoods lead to more powerful tests,
as expected. However, Figure 2(b) shows that we face diminishing returns when
trading off compute budget against test power. The CSK kernel always does
better than the Hamming kernel. Interestingly, for the CSK kernel, the simpler
ZS’ test performs almost as well as the ZS test. Both of these Stein tests
outperform the CSK MMD for $J=\infty$.

Another approach to controlling the size and connectivity of the neighbourhood
is to control the symbol neighbourhoods $\mathcal{N}_{\text{ins}}$ and
$\mathcal{N}_{\text{rep}}$. We use the operator $\mathrm{ZS}^{\prime}$ which
uses a neighbourhood consisting of substitutions, choose a cyclic structure
over the alphabet $\mathcal{S}$, and limit the symbol neighbourhood
$\mathcal{N}_{\text{rep}}\left(s\right)$ using the distance in the cyclic
alphabet.

$1$$2$$3$$0$$0.2$$0.4$$0.6$$\cdots$$\infty$Symbol neighbourhood size parameter
$d$Estimated rejection rateHamming $\mathrm{ZS}^{\prime},J=1$Hamming
$\mathrm{ZS}^{\prime},J=3$CSK $\mathrm{ZS}^{\prime},J=1$CSK
$\mathrm{ZS}^{\prime},J=3$Level $\alpha=0.05$ Figure 3: Estimated rejection
rate of several Zanella-Stein tests based on the substitution neighbourhood.
We control the size of the symbol neighbourhood $\mathcal{N}_{\text{rep}}$ and
estimate rejection rate as in Figure 2.

In Figure 3 we can see no notable impact of symbol neighbourhood on test
power. This shows that we can construct tests using a sparse neighbourhood. An
infinite alphabet could be addressed with such an approach, as described
earlier. The experiment also shows that changes to the size of the
neighbourhood are not sufficient to fully explain our previous results. Test
power is evidently also impacted by the location of edits rather than just the
number.

###  4.3 Use of Evidence in Long Sequences

Edits across the entire sequence lead to higher test power. One explanation
for this effect is that editing the entire sequence, instead of just editing
near the end, allows the test to make better use of the evidence available
from longer sequences. Where the data are distributed according to a Markov
chain or some other model that mixes quickly, longer sequences provide for
more information about the distribution. We now show that the proposed test
can indeed make use of such information.

In order to illustrate this, we construct the following problem: The model
distribution and ground truth are second-order Markov chains over a finite
alphabet with $|\mathcal{S}|=8$. We select two transition kernels randomly
from a Dirichlet distribution with concentration parameter $\alpha=1$. The
model distribution uses one of the two kernels. The ground truth uses an equal
mixture of the two kernels, so that at each step either of the kernels is used
with equal probability. The two chains both stop with a fixed probability
$1/\lambda$ after each step. We fix the number of samples at 8 and vary the
stopping probability in order to control the length of the samples.

The KSD test for this experiment is constructed using a large, strongly
connected, neighbourhood: specifically, the $J$-location modification
neighbourhood with $J=\infty$, which allows insertion, deletion, and
substitution anywhere in the sequence. As a baseline, we compare against an
MMD test using 100 samples from the model distribution, as well as two
likelihood-ratio tests: an oracle test, where $H_{1}:\\{P,Q\\}$, and a test
where the alternative consists of all second-order Markov chains over
$\mathcal{S}$. The two kernel tests are constructed using the CSK kernel with
subsequence length $n=3$, and we use parametric bootstrap for all tests.

The results are shown in Figure 4, and we can see that the test is able to
reject more often on longer sequences, without the need for more samples. The
test outperforms the MMD baseline.

$5$$10$$15$$20$$0$$0.2$$0.4$$0.6$$0.8$$1$Expected length
$\mathbb{E}\left[\operatorname{len}\left(X\right)\right]=1/\lambda$Estimated
rejection rateKSDMMDLR oracleLR MCLevel $\alpha=0.05$ Figure 4: Estimated
rejection rate of the proposed test, applied to sequences of increasing
length. The model distribution and ground truth are second order Markov
models. We control the expect length of the samples, while holding the sample
count fixed at 8. As the expected length increases, test power grows.

###  4.4 Sensitivity to Changes in Length

There is a caveat with the family of neighbourhoods we have described: these
neighbourhoods are not particularly sensitive to changes in the distribution
over sequence length
$\mathcal{L}\left(\operatorname{len}\left(X\right)\right)$. We demonstrate
this by choosing a randomly generated Markov chain over a finite alphabet
$|\mathcal{S}|=10$ as the model distribution, with the transition kernel again
chosen from a Dirichlet distribution as above. The ground truth and model
distribution use the same transition kernel, but we vary the stopping
probability. The stopping probability is selected such that the expected
sequence length under the model distribution is 8, and we vary the sequence
length under the ground truth between 8 and 20.

We use the same tests as in the previous experiment, except that we use a
subsequence length of 2 for the CSK kernel and the alternative for the LR
baseline now consists of first- rather than second-order Markov chains.

$0$$0.2$$0.4$$0.6$$0.8$$1$$0$$0.2$$0.4$$0.6$$0.8$$1$Degree of perturbation in
the length distribution
$\mathcal{L}\left(\operatorname{len}\left(X\right)\right)$Estimated rejection
rateKSDMMDLR oracleLR MCLevel $\alpha=0.05$ Figure 5:  Estimated rejection
rate of the proposed KSD test, compared with three baselines: MMD, a suitable
likelihood-ratio test, and an oracle LR test. We use a randomly generated
Markov chain as the model distribution and perturb the distribution of
lengths, $\mathcal{L}\left(\operatorname{len}\left(X\right)\right)$, for the
ground truth. The proposed KSD test fails to reject in this setting.

In Figure 5, we can see that the KSD test fails to reject entirely. Both types
of baseline significantly outperform.

This experimental result highlights the need for bespoke tests, when we desire
that a test is sensitive to a specific aspect of the distribution. A test user
who is particularly interested in the length distribution would likely benefit
from directly testing the distribution over lengths using a classical test.

$0.7$$0.8$$0.9$$1$$1.1$$1.2$$1.3$$0$$0.2$$0.4$$0.6$$0.8$$1$Concentration
parameter $\alpha$ of ground truthEstimated rejection rateKSDMMDLR oracleLR
MCLevel $\alpha=0.05$ Figure 6:  Experiment showing performance of the
proposed test on a simple MRF model. Our toy MRF model permits a normalized
density and we can sample effectively from it. We use this to construct three
baselines that would not be available in the general case. Perturbation is
controlled via a concentration parameter $\alpha$, which is set fo $1$ for the
model distribution and ranges from $0.75$ to $1.25$ for the ground truth. The
concentration parameter controls the correlation of successive elements in a
sequence, so that repetitions become more (less) common as $\alpha$ grows
(shrinks). Our proposed test outperforms both non-oracle baselines in one
region and outperforms MMD in the other region.

###  4.5 Markov Random Fields

A major benefit of the proposed test is that it does not require us to sample
from the model distribution, nor to compute a normalized density. It is
nonetheless worthwhile to compare against baselines where normalised densities
are known, and/or for which sampling is possible. We investigate this setting
using a Markov Random Field (MRF) which allows us both to compute a normalized
density and to sample from the distribution. This permits us to draw i.i.d.
samples for the experiment as well as evaluate the MMD and LR baselines.
However, the KSD tests will make no use of the normalized density nor of the
ability to sample; the MMD tests will only make use of the ability to sample
in order to perform the required two-sample test. We no longer perform
parametric bootstrap for the kernel tests but rather wild bootstrap, as the
latter does not require sampling.

For the MRF model we specify a single potential across the entire sample space
$\mathcal{X}$:

$\displaystyle h\left(x_{1:l}\right)$
$\displaystyle\coloneqq\begin{cases}C\cdot
l+\alpha\sum_{i=1}^{l-1}\delta_{x_{i}}\left(x_{i+1}\right),&\text{if }l\leq
M\\\ -\infty,&\text{else.}\end{cases}$ (7)

In our experiment we hold the length parameters $C$ and $M$ fixed, while
perturbing the distribution by varying $\alpha$. The parameter $\alpha$ is a
concentration parameter and controls the correlation between successive
elements of a sequence. We specify $\alpha=1$ for the model distribution and
vary $\alpha\in[0.75,1.25]$ for the ground truth. We specify $M=20$ in all
cases, and vary $C$ with $\alpha$ so that the mean length is fixed at 10 and
the length distribution does not depend on $\alpha$. The alphabet size is 3.

We compare the proposed test against the MMD and LR baselines. For the LR
baseline we use an oracle test, as well as a test where $H_{1}$ consists of
first-order Markov chains. This baseline does not contain the toy model MRF
family in the alternative. Note that a general likelihood-ratio test for MRF
models is intractable. The MMD and KSD tests use the CSK kernel with
subsequence length 2.

The proposed test performs competitively with the non-oracle baselines. In
Figure 6 we can see that the proposed test outperforms the MMD baseline across
the full range of perturbations, and outperforms the non-oracle LR baseline in
one region but not the other. The LR baseline fails as $\alpha$ grows, due to
lack of consistency. We emphasize that our proposed test applies more
generally than the LR test, since we do not require normalized densities.
Moreover, using the wild Bootstrap approach, we do not need to sample from the
model, unlike the MMD test.

$200$$300$$400$$500$$600$$700$$0$$0.2$$0.4$$0.6$$0.8$$1$Number of
samplesEstimated rejection rateKSDMMDLR oracleLR MCLevel $\alpha=0.05$ Figure
7:  Performance of the proposed KSD test on a simple MRF model. We perform the
same experiment as shown in Figure 6, in this case holding fixed $\alpha=0.9$
and varying the number of samples. In this region, our proposed test
outperforms the MMD baseline but not the LR baseline.

We also show how test power varies with the data size in Figure 7. We choose
the region where our test outperforms the MMD baseline, but not the LR
baseline. The concentration parameter is fixed to $\alpha=0.9$, while we vary
the sample size $n$. We can see that the MMD baseline has a low sample
efficiency for this problem, confirming again that the MMD is unable to use
the model structure. Both the LR baseline and our proposed KSD test make use
of model information. Only the KSD test can work in the absence of a
normalized density, however.

#### Acknowledgements

Heishiro Kanagawa and Arthur Gretton thank the Gatsby Charitable Foundation
for the financial support.

## References

  * [1] Wassily Hoeffding  “A Class of Statistics with Asymptotically Normal Distribution”  In _Ann. Math. Statist._ 19.3, 1948, pp. 293–325
  * [2] Nachman Aronszajn  “Theory of reproducing kernels”  In _Transactions of the American mathematical society_ 68.3, 1950, pp. 337–404
  * [3] AA Barker  “Monte Carlo Calculations of the Radial Distribution Functions for a Proton?Electron Plasma”  In _Australian Journal of Physics_ 18.2 CSIRO Publishing, 1965, pp. 119  DOI: 10.1071/ph650119
  * [4] Charles M. Stein  “A bound for the error in the normal approximation to the distribution of a sum of dependent random variables”, 1972
  * [5] Andrew D Barbour  “Stein’s method and Poisson process convergence” Publisher: Cambridge University Press  In _Journal of Applied Probability_ 25, 1988, pp. 175–184
  * [6] Herold Dehling and Thomas Mikosch  “Random quadratic forms and the bootstrap for U-statistics” Publisher: Elsevier  In _Journal of Multivariate Analysis_ 51.2, 1994, pp. 392–413
  * [7] D. Blei, A. Ng and M. Jordan  “Latent Dirichlet Allocation”  In _Journal of Machine Learning Research_ 3, 2003, pp. 993–1022
  * [8] Ingo Steinwart and Andreas Christmann  “Support Vector Machines”, Information Science and Statistics  Springer, 2008
  * [9] Robert J. Serfling  “Approximation theorems of mathematical statistics”  John Wiley & Sons, 2009
  * [10] Bharath K. Sriperumbudur, Kenji Fukumizu and Gert R.. Lanckriet  “Universality, characteristic kernels and RKHS embedding of measures”  In _Journal of Machine Learning Research_ 12, 2011, pp. 2389–2410
  * [11] Arthur Gretton et al.  “A kernel two-sample test” Publisher: JMLR. org  In _The Journal of Machine Learning Research_ 13.1, 2012, pp. 723–773
  * [12] Jackson Gorham and Lester Mackey  “Measuring sample quality with Stein’s method”  In _Advances in Neural Information Processing Systems_ 28, 2015
  * [13] J. Lloyd and Z. Ghahramani  “Statistical model criticism using kernel two sample tests”  In _Advances in Neural Information Processing Systems_ , 2015, pp. 829–837
  * [14] Kacper Chwialkowski, Heiko Strathmann and Arthur Gretton  “A kernel test of goodness of fit”  In _International conference on machine learning_ PMLR, 2016, pp. 2606–2615
  * [15] Qiang Liu, Jason Lee and Michael Jordan  “A kernelized Stein discrepancy for goodness-of-fit tests”  In _International conference on machine learning_ PMLR, 2016, pp. 276–284
  * [16] Chris J. Oates, Mark Girolami and Nicolas Chopin  “Control functionals for Monte Carlo integration”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ 79.3 Wiley, 2016, pp. 695–718  DOI: 10.1111/rssb.12185
  * [17] Jackson Gorham and Lester Mackey  “Measuring Sample Quality with Kernels”  In _Proceedings of The 34th International Conference on Machine Learning_ , 2017, pp. 1292–1301
  * [18] Jiasen Yang, Qiang Liu, Vinayak A. Rao and Jennifer Neville  “Goodness-of-fit Testing for Discrete Distributions via Stein Discrepancy”  In _ICML_ , 2018
  * [19] Alessandro Barp et al.  “Minimum Stein Discrepancy Estimators”  In _Advances in Neural Information Processing Systems_ 32, 2019
  * [20] Guy Bresler and Dheeraj Nagaraj  “Stein’s Method for Stationary Distributions of Markov Chains and Application to Ising Models”  In _Annals of Applied Probability_ 29, 2019
  * [21] Heishiro Kanagawa et al.  “A kernel stein test for comparing latent variable models”  In _arXiv preprint arXiv:1907.00586_ , 2019
  * [22] Samuel Power and Jacob Vorstrup Goldman  “Accelerated sampling on discrete spaces with non-reversible Markov processes”  In _arXiv preprint arXiv:1912.04681_ , 2019
  * [23] Gesine Reinert and Nathan Ross  “Approximating stationary distributions of fast mixing Glauber dynamics, with applications to exponential random graphs”  In _Annals of Applied Probability_ 29, 2019
  * [24] Jiasen Yang, Vinayak Rao and Jennifer Neville  “A Stein–Papangelou Goodness-of-Fit Test for Point Processes”  In _Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics_ , 2019, pp. 226–235
  * [25] Giacomo Zanella  “Informed Proposals for Local MCMC in Discrete Spaces”  In _Journal of the American Statistical Association_ 115, 2019, pp. 852–865
  * [26] Tamara Fernandez, Nicolas Rivera, Wenkai Xu and Arthur Gretton  “Kernelized Stein Discrepancy Tests of Goodness-of-fit for Time-to-Event Data”  In _Proceedings of the 37th International Conference on Machine Learning_ , 2020, pp. 3112–3122
  * [27] Liam Hodgkinson, Robert Salomone and Fred Roosta  “The reproducing Stein kernel approach for post-hoc corrected sampling”  In _arXiv preprint arXiv:2001.09266_ , 2020
  * [28] Wenkai Xu and Takeru Matsuda  “A Stein Goodness-of-fit Test for Directional Distributions”  In _Proceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics_ , 2020, pp. 320–330
  * [29] Andreas Anastasiou et al.  “Stein’s Method Meets Statistics: A Review of Some Recent Developments”  In _arXiv preprint arXiv:2105.03481_ , 2021
  * [30] Jiaxin Shi et al.  “Gradient Estimation with Discrete Stein Operators”  In _ArXiv_ abs/2202.09497, 2022
  * [31] George Wynne, Mikołaj Kasprzak and Andrew B. Duncan  “A Spectral Representation of Kernel Stein Discrepancy with Application to Goodness-of-Fit Tests for Measures on Infinite Dimensional Hilbert Spaces”  arXiv, 2022  DOI: https://doi.org/10.48550/arXiv.2206.04552

##  Appendix A PROGRAM CODE

The program code used to produce our experimental results is available online
at: https://github.com/test-for-sequential-models/code

##  Appendix B EXPERIMENT DETAILS

We provide further details for the experiments that we conducted.

###  B.1 Likelihood Ratio Test

The LR baselines use the likelihood ratio statistic:

$\displaystyle
T^{\text{LR}}_{m}\left(\left\\{Y^{(l)}\right\\}_{l=1}^{m}\right)$
$\displaystyle\coloneqq 2\left(\sup_{p\in H_{0}\cup H_{1}}\log
p\left(\left\\{Y^{(l)}\right\\}_{l=1}^{m}\right)-\sup_{p\in H_{0}}\log
p\left(\left\\{Y^{(l)}\right\\}_{l=1}^{m}\right)\right),$ (8)

where $H_{0}$ and $H_{1}$ are the null and alternative hypotheses,
respectively. In our experiments, we use a simple null hypothesis for the
likelihood ratio statistic. The alternative used in the oracle baseline is a
simple hypothesis containing only the (normally unknown) ground truth. The
alternative used in the MC baseline is a composite hypothesis consisting of
all Markov chains of some order, with the order depending on the experiment.
In experiments where we use an alternative consisting of first-order Markov
chains, the alternative consists of all initial distributions over the
alphabet together with all transition kernels through the alphabet, with an
additional stopping point:

$\displaystyle H_{1}$
$\displaystyle\coloneqq\operatorname{Sim}_{\mathcal{S}}\times\left(\operatorname{Sim}_{\mathcal{S}\cup\\{\text{stop}\\}}\right)^{\mathcal{S}},$
(9) $\displaystyle\operatorname{Sim}_{\mathcal{S}}$
$\displaystyle\coloneqq\left\\{(\theta_{i})_{i\in\mathcal{S}}:\sum_{i\in\mathcal{S}}\theta_{i}=1\text{
and }\theta_{i}\geq 0\,\text{ for any }i\in\mathcal{S}\right\\}.$ (10)

In experiments where we use an alternative consisting of second-order Markov
chains, the alternative consists of the analogous collection of second-order
kernels:

$\displaystyle H_{1}$
$\displaystyle\coloneqq\operatorname{Sim}_{\mathcal{S}}\times\left(\operatorname{Sim}_{\mathcal{S}\cup\\{\text{stop}\\}}\right)^{\mathcal{S}}\times\left(\operatorname{Sim}_{\mathcal{S}\cup\\{\text{stop}\\}}\right)^{\mathcal{S}\times\mathcal{S}}.$
(11)

###  B.2 Details of Testing Scenarios

In several experiments, we referred to a collection of 12 synthetic testing
scenarios to evaluate the test power. We provide the details of these
scenarios. The scenarios are intended to be dissimilar to each other in form
and difficulty, although all model distributions and ground truths are Markov
chains of some order.

Some of the distributions described here do not have support on the entire
sequence space. In our earlier discussion, we required support everywhere as a
technical condition relating to the connectivity of the Zanella process. To
ensure that the probability of all sequences is positive, we introduce low-
probability “restart” events into all of our Markov chains that resample
uniformly from the alphabet with low probability, rather than sampling from
the specified transition kernel. The probability of these events is small,
$\varepsilon\approx 10^{-3}$. We do not mention this further, in order to
simplify our discussion.

##### Binary i.i.d. sequences

We draw sequences of i.i.d. Bernoulli variables. The model and ground truth
differ in the parameter of the Bernoulli distribution.

Alphabet: $\\{0,1\\}$.

Model distribution: Length follows a Poisson distribution with mean 20.
Contents of the sequence are i.i.d. Bernoulli with $p=0.6$.

Ground truth: Length follows a Poisson distribution with mean 20. Contents of
the sequence are i.i.d. Bernoulli with $p=0.4$.

Sample size: 10.

##### Binary sequences: misspecified Markov order

We simulate the case where we are mistaken about the Markov order of the data.
The model consists of i.i.d. Bernoulli variables, while the ground truth
follows a simple first-order Markov chain.

Alphabet: $\\{0,1\\}$.

Model distribution: Length follows a geometric distribution with mean 20.
Contents of the sequence are i.i.d. Bernoulli with $p=0.6$.

Ground truth: First-order Markov chain. Initial distribution uniform on {0,
1}. The transition kernel is

$\displaystyle\operatorname{Pr}\left(\text{stop}\middle|X_{t}=x\right)$
$\displaystyle=1/\lambda,$ (12)
$\displaystyle\operatorname{Pr}\left(X_{t+1}=y\middle|X_{t}=x\right)$
$\displaystyle=\left(1-1/\lambda\right)\delta_{x}(y),$ (13)

with $\lambda=20$. This gives sequences alternating between $\\{0,1\\}$.

Sample size: 30.

##### Random walk I

We take a random walk through a cyclic alphabet and perturb it by letting the
ground truth hold, or “skip a step”, with low probability.

Alphabet: $\mathcal{S}=\\{1,\ldots,8\\}$ with a cyclic structure.

Model distribution: Random walk through $\mathcal{S}$ according to the cyclic
structure. Start uniformly on $\mathcal{S}$. After each step, stop with
probability $1/\lambda$, with $\lambda=8$.

Ground truth: Random walk with holding. Start uniformly on $\mathcal{S}$. At
each step, hold or “skip a step”, with probability $p=0.2$, letting
$X_{t+1}=X_{t}$. Otherwise, sample a step from $\\{-1,+1\\}$ uniformly. After
each step, stop with probability $1/\lambda$, with $\lambda=8$.

Sample size: 30.

##### Random walk II

The same scenario, but with few long sequences instead of many short
sequences. We only let the ground truth hold in a subset of states, which
reduces the degree of perturbation and makes the problem more challenging.

Alphabet: $\mathcal{S}=\\{1,\ldots,30\\}$ with a cyclic structure.

Model distribution: Same as in Random walk I with $\lambda=30$.

Ground truth: Same as in Random walk I but with $\lambda=30$ and sampling the
step from $\\{-1,0,+1\\}$ according to:

$\displaystyle\operatorname{Pr}\left(X_{t+1}-X_{t}=\Delta\middle|X_{t}=x\right)$
$\displaystyle=\begin{cases}p,&\text{if }\Delta=0,x\leq i_{\text{max}},\\\
\frac{1-p}{2},&\text{if }\Delta\in\\{-1,+1\\},x\leq i_{\text{max}},\\\
\frac{1}{2},&\text{if }\Delta\in\\{-1,+1\\},x>i_{\text{max}},\\\
0,&\text{otherwise,}\end{cases}$ (14)

with $i_{\text{max}}=8$ and $p=0.2$.

Sample size: 8.

##### Random walk III

We introduce memory into the random walk by correlating or anti-correlating
successive steps, which gives us a simple second-order MC.

Alphabet: $\mathcal{S}=\\{1,\ldots,10\\}$ with a cyclic structure.

Model distribution: Same as in Random walk I, but we modify the transition
kernel to correlate successive steps. Specifically, letting $\delta=0.95$, we
sample the step from $\\{-1,0,+1\\}$ according to:

$\displaystyle\operatorname{Pr}\left(X_{t+2}-X_{t+1}=d\middle|X_{t+1},X_{t}\right)$
$\displaystyle=\begin{cases}\delta,&\text{if }d=X_{t+1}-X_{t},|d|=1,\\\
1-\delta,&\text{if }d=X_{t}-X_{t+1},|d|=1,\\\ \frac{1}{2},&\text{if
}|d|=1<|X_{t}-X_{t+1}|,\\\ 0,&\text{otherwise.}\end{cases}$ (15)

We have to treat the case where $|X_{t}-X_{t+1}|>1$ due to the low probability
restart event mentioned earlier.

Ground truth: Same as model distribution but with $\delta=0.05$.

Sample size: 30.

##### Random walk IV

Same as in Random walk III but with $\lambda=30$ and sample size 8. That is,
few long sequences instead of many short ones.

##### Random second-order MC I

We define a distribution over second-order MC and sample two random ones. The
stopping probability is constant so that we can control the length
independently of the perturbation of the contents. We sample many short
sequences.

Alphabet: $\mathcal{S}=\\{1,\ldots,10\\}$.

Model distribution: We sample a second-order Markov chain by sampling the
initial distribution and transition kernels from Dirichlet distributions with
concentration parameter $\alpha=1$. At each step, the chain stops with
probability $1/\lambda$ where $\lambda=8$.

Ground truth: Same as model distribution, but the Markov chain is sampled from
a different random number generator (RNG) seed in order to give a different
MC.

Sample size: 30.

##### Random second-order MC II

Same as in Random second-order MC I but with $\lambda=20$ and sample size 8.
That is, we sample a small number of long sequences.

##### Random second-order MC III

Same as in Random second-order MC I but with $\lambda=8$ and sample size 8.
That is, we sample a small number of short sequences.

##### Random MC with varied initial distribution I

We sample a single random first-order MC for the model and ground truth, and
set a fixed initial distribution which we perturb. This poses a hard problem,
as most of the observed data relates to transitions that come from an
unperturbed transition kernel.

Alphabet: $\mathcal{S}=\\{1,\ldots,10\\}$.

Model distribution: We sample a first-order Markov chain by sampling the
transition kernel from a Dirichlet distribution with concentration parameter
$\alpha=1$. The initial distribution is uniform over $\mathcal{S}$. At each
step, the chain stops with probability $1/\lambda$ where $\lambda=8$.

Ground truth: Same as model distribution, with the same RNG seed. in order to
give the same transition kernel. The initial distribution is an equal mixture
of $\operatorname{Unif}\left(\mathcal{S}\right)$ and
$\operatorname{Unif}\left(\\{1,2\\}\right)$.

Sample size: 30.

##### Random MC with varied initial distribution II

Same as in Random MC with varied initial distribution I but with $\lambda=20$
and sample size 8. This problem is significantly harder, as we observe 8
rather than 30 instances of the initial distribution, and the longer sequences
do not provide any additional evidence.

##### Random MC with varied length distribution

We vary the distribution over lengths while keeping the content distribution
the same.

Alphabet: $\mathcal{S}=\\{1,\ldots,10\\}$.

Model distribution: We sample a first-order Markov chain by sampling the
transition kernel from a Dirichlet distribution with concentration parameter
$\alpha=1$. The initial distribution is uniform over $\mathcal{S}$. At each
step, the chain stops with probability $1/\lambda$ where $\lambda=8$.

Ground truth: We use the same initial distribution and transition kernel as in
the model. At each step, the chain stops with probability $1/\lambda$ where
$\lambda=20$.

Sample size: 30.