# Data thinning for convolution-closed distributions

Anna Neufeld  Department of Statistics, University of Washington  Ameer
Dharamshi  Department of Biostatistics, University of Washington  Lucy L. Gao
Department of Statistics, University of British Columbia  Daniela Witten
Department of Statistics, University of Washington  Department of
Biostatistics, University of Washington

###### Abstract

We propose _data thinning_ , a new approach for splitting an observation into
two or more independent parts that sum to the original observation, _and that
follow the same distribution as the original observation_ , up to a (known)
scaling of a parameter. This proposal is very general, and can be applied to
any observation drawn from a “convolution closed” distribution, a class that
includes the Gaussian, Poisson, negative binomial, Gamma, and binomial
distributions, among others. It is similar in spirit to — but distinct from,
and more easily applicable than — a recent proposal known as _data fission_.
Data thinning has a number of applications to model selection, evaluation, and
inference. For instance, cross-validation via data thinning provides an
attractive alternative to the “usual” approach of cross-validation via sample
splitting, especially in unsupervised settings in which the latter is not
applicable. In simulations and in an application to single-cell RNA-sequencing
data, we show that data thinning can be used to validate the results of
unsupervised learning approaches, such as $k$-means clustering and principal
components analysis.

##  1 Introduction

As datasets continue to grow in size, scientists fit increasingly complex
models to their data. Thus, there is an ever-growing need for out-of-box
methods that can be used to validate the outputs of complex statistical
learning models. In many settings, the simplest and most natural way to do
this is via _sample splitting_ , in which the observations are split into a
training set, used to fit a model, and a test set, used to validate it [Hastie
et al., 2009]. Sample splitting can also be applied to conduct inference after
model selection, as in Fithian et al. [2014]. It is an invaluable tool in the
arsenal of any practicing data scientist, and can be applied to virtually any
supervised learning problem.

In recent work, Leiner et al. [2022] proposed an attractive alternative to
sample splitting, which they call _data fission_. Data fission involves
decomposing a _single observation_ $X$ into two parts, $X^{(1)}$ and
$X^{(2)}$, such that (i) $X^{(1)}$ and $X^{(2)}$ can together be used to
reconstruct $X$, and (ii) the joint distribution of
$\left(X^{(1)},X^{(2)}\right)$ is tractable. In the special cases where $X$ is
drawn from either a Gaussian or a Poisson distribution, they propose a way to
obtain $X^{(1)}$ and $X^{(2)}$ that are independent; see also independent
proposals by Rasines and Young [2021], Oliveira et al. [2021], Chen et al.
[2021], and Neufeld et al. [2022]. However, outside of these two special
cases, Leiner et al. [2022] do not propose an approach for generating
$X^{(1)}$ and $X^{(2)}$ that are independent. Instead, they propose a
strategy, using conjugate priors, to obtain $X^{(1)}$ and $X^{(2)}$ such that
the distributions of both $X^{(1)}$ and $X^{(2)}\mid X^{(1)}$ are known.
Leiner et al. [2022] argue that we can use knowledge of the distributions of
$X^{(1)}$ and $X^{(2)}\mid X^{(1)}$ in settings where sample splitting would
typically be used: e.g. we can fit a lasso-penalized model using $X^{(1)}$,
and can then conduct inference on the selected coefficients using $X^{(2)}\mid
X^{(1)}$. In principle, though not a primary focus of their paper, to evaluate
the suitability of a model for $X$, we can fit the model using $X^{(1)}$, and
can validate it using $X^{(2)}\mid X^{(1)}$ — however, as we will see, this
approach is hard to apply in practice.

In many cases, the strategy of Leiner et al. [2022] leads to decompositions
where the distributions of $X^{(1)}$ and $X^{(2)}\mid X^{(1)}$ are quite
different from each other, and from the distribution of $X$. We can see this
in the following example, in which a single Gamma random variable is split
into a $B$-vector of geometric random variables, and a variable that
conditionally follows a Gamma distribution.

######  Example 1.1 (Gamma decomposition, data fission approach).

Suppose $X\sim\mathrm{Gamma}(\alpha,\beta)$, where
$\operatorname{\text{E}}[X]=\alpha/\beta$ and
$\operatorname{\text{Var}}(X)=\alpha/(\beta^{2})$. For a tuning parameter
$B\in\\{1,2,\ldots\\}$, Leiner et al. [2022] propose drawing
$Z=(Z_{1},\ldots,Z_{B})$, where
$Z_{i}\overset{\mathrm{ind.}}{\sim}\mathrm{Poisson}(X)$, and thus marginally
$Z_{i}\overset{\mathrm{ind.}}{\sim}\mathrm{NegativeBinomial}(\alpha,\beta/(\beta+1))$.
Take $X^{(1)}=Z$, and $X^{(2)}=X$. Then, the conditional distribution of
$X^{(2)}\mid X^{(1)}$ is $\mathrm{Gamma}(\alpha+\sum_{i=1}^{B}Z_{i},\beta+B)$.

In this example, the distributions of $X$, $X^{(1)}$, and $X^{(2)}\mid
X^{(1)}$ are all quite different. This poses challenges in using this
decomposition in practice: for instance, in evaluating the suitability of a
model for $X$ by fitting it to $X^{(1)}$ and validating it using $X^{(2)}\mid
X^{(1)}$.

In this paper, we propose _data thinning_ , a recipe for decomposing an
observation $X$ into two parts, $X^{(1)}$ and $X^{(2)}$, such that (i)
$X=X^{(1)}+X^{(2)}$, (ii) $X^{(1)}$ and $X^{(2)}$ are independent, and (iii)
$X^{(1)}$ and $X^{(2)}$ follow the same distribution as $X$, up to a (known)
scaling of a parameter. Critically, properties (ii) and (iii) guarantee that
once we have obtained $X^{(1)}$ and $X^{(2)}$, _we know what to do with them_
, since they are independent draws from the distribution used to generate $X$.
For instance, to evaluate the suitability of a model for $X$, we can fit it to
$X^{(1)}$, and can validate it using $X^{(2)}$. Our recipe can be applied to
any distribution that is _convolution-closed_ , in the sense of Joe [1996]:
this includes the Gaussian, Poisson, negative binomial, Gamma, exponential,
and binomial distributions, among others. Our proposal coincides with that of
Leiner et al. [2022] only in the special case of the Poisson distribution [see
also Neufeld et al., 2022].

To see the benefit of our data thinning approach over the data fission
proposal of Leiner et al. [2022], we return to the setting of Example 1.1.

######  Example 1.2 (Gamma decomposition, data thinning approach).

Suppose $X\sim\mathrm{Gamma}(\alpha,\beta)$, with $\beta$ unknown. For a
tuning parameter $\epsilon\in(0,1)$, we take $X^{(1)}=XZ$ and
$X^{(2)}=X-X^{(1)}$, where
$Z\sim\mathrm{Beta}(\epsilon\alpha,(1-\epsilon)\alpha)$. Then $X^{(1)}$ and
$X^{(2)}$ are independent, and
$X^{(1)}\sim\mathrm{Gamma}(\epsilon\alpha,\beta)$ and
$X^{(2)}\sim\mathrm{Gamma}((1-\epsilon)\alpha,\beta)$.

In other words, data thinning allows us to decompose a
$\mathrm{Gamma}(\alpha,\beta)$ random variable, for which the parameter
$\beta$ is unknown, into two independent Gamma random variables, $X^{(1)}$ and
$X^{(2)}$. Therefore, fitting a model to $X^{(1)}$ and validating it using
$X^{(2)}$ is straightforward.

######  Remark 1.

Comparing Examples 1.1 and 1.2, we see that data fission does not require
knowledge of any parameters, whereas data thinning requires knowledge of the
shape parameter $\alpha$. In the special case of the exponential distribution,
this is not an issue, as we know that $\alpha=1$ (because
$\mathrm{Exponential}(\lambda)$ is the same as $\mathrm{Gamma}(1,\lambda)$).
Moreover, note that $\alpha$ appears in the conditional distribution of
$X^{(2)}\mid X^{(1)}$ in Example 1.1. Thus, if the goal is to do inference on
$\beta$, then fruitful application of the data fission framework will also
require a good estimate of $\alpha$. We explore the impact of using the
incorrect value of $\alpha$ in data thinning in Section 2.4.

Furthermore, unlike the data fission proposal of Leiner et al. [2022], data
thinning allows us to _decompose an observation into an arbitrary number of
independent components_. To see this, we continue Example 1.2.

######  Example 1.3 (Gamma decomposition into $M$ components, data thinning
approach).

Once again, we suppose that $X\sim\mathrm{Gamma}(\alpha,\beta)$, where $\beta$
is unknown. We take $(X^{(1)},\ldots,X^{(M)})=XZ$, where
$Z\sim\mathrm{Dirichlet}(\alpha/M,\ldots,\alpha/M)$. Then
$X^{(1)},\ldots,X^{(M)}$ are independent, they sum to $X$, and each is
marginally drawn from a $\mathrm{Gamma}(\alpha/M,\beta)$ distribution.

Thus, for $m=1,\ldots,M$, it is straightforward to fit a model to
$\sum_{m^{\prime}\neq m}X^{(m^{\prime})}$ and validate it using $X^{(m)}$.
This enables us to develop an analog to $K$-fold cross-validation using sample
splitting that is based on data thinning.

Data thinning has a number of potential applications, including the post-
selection inference examples that were the primary focus of Leiner et al.
[2022]. In this paper, we instead focus on its use in model assessment and
validation. Specifically, we focus on the unsupervised learning setting, in
which the “usual” cross-validation via sample splitting approach cannot be
directly applied [see, e.g. Owen and Perry, 2009]. We show that _cross-
validation via data thinning_ provides an attractive alternative in this
setting.

The rest of this paper is organized as follows. In Section 2, we briefly
review the class of convolution-closed distributions, and introduce a thinning
procedure to split a single random variable into two independent random
variables, each of which follows the same distribution as the original random
variable (up to a scaling of the parameter(s)). We further extend this
thinning procedure to split a single random variable into an arbitrary number
of independent random variables in Section 3. In Section 4, we compare our
decompositions to those of Leiner et al. [2022]. Section 5 contains a
simulation study, in which we show that cross-validation via data thinning
provides an alternative to cross-validation via sample splitting in
unsupervised settings. An application to single-cell RNA-sequencing data is in
Section 6. We close with a discussion in Section 7. All proofs are in the
supplementary materials.

##  2 The data thinning proposal

###  2.1 A review of convolution-closed distributions

We begin by defining a convolution-closed distribution; see Joe [1996] and
Jørgensen and Song [1998] for additional discussion.

######  Definition 1 (Convolution-closed).

Let $F_{\lambda}$ denote a distribution indexed by a parameter $\lambda$ in
parameter space $\Lambda$. Let ${X}^{\prime}\sim F_{\lambda_{1}}$ and
${X}^{\prime\prime}\sim F_{\lambda_{2}}$ with
${X}^{\prime}\perp\\!\\!\\!\perp{X}^{\prime\prime}$. If
${X}^{\prime}+{X}^{\prime\prime}\sim F_{\lambda_{1}+\lambda_{2}}$ whenever
$\lambda_{1}+\lambda_{2}\in\Lambda$, then $F_{\lambda}$ is _convolution-
closed_ in the parameter $\lambda$.

Many well-known distributions are convolution-closed. While the
$\mathrm{Poisson}(\lambda)$ distribution is convolution-closed in its single
parameter $\lambda$ and the $\mathrm{N}(\mu,\sigma^{2})$ distribution is
convolution-closed in the two-dimensional parameter $(\mu,\sigma^{2})$, other
distributions, such as the Gamma, are convolution-closed in just one parameter
with the other parameter(s) held fixed. Table 1 provides details about some
well-known convolution-closed distributions.

######  Remark 2 (Expectation is linear in $\lambda$).

Let $F_{\lambda}$ for $\lambda\in\Lambda$ be a convolution-closed family. By
definition, if $X^{\prime}\sim F_{\lambda_{1}}$ and $X^{\prime\prime}\sim
F_{\lambda_{2}}$ and $\lambda_{1}+\lambda_{2}\in\Lambda$, then
$X^{\prime}+X^{\prime\prime}\sim F_{\lambda_{1}+\lambda_{2}}$. If these
distributions have first moments, then
$\operatorname{\text{E}}[X^{\prime}+X^{\prime\prime}]=\operatorname{\text{E}}[X^{\prime}]+\operatorname{\text{E}}[X^{\prime\prime}]$.
Thus, outside of contrived counterexamples, $\operatorname{\text{E}}[X]$ is a
linear function of the parameter $\lambda$. We say that a distribution whose
expectation $\operatorname{\text{E}}[X]$ is linear in the convolution-closed
parameter $\lambda$ satisfies the _linear expectation property_. All
distributions in Table 1 satisfy this property.

Distribution | Notes
---|---
$X\sim\mathrm{Poisson}(\lambda)$, where $\operatorname{\text{E}}[X]=\lambda$ and $\operatorname{\text{Var}}(X)=\lambda$. |  Convolution-closed in $\lambda$.
$X\sim\mathrm{N}(\mu,\sigma^{2})$, where $\operatorname{\text{E}}[X]=\mu$ and $\operatorname{\text{Var}}[X]=\sigma^{2}$. |  Convolution-closed in $(\mu,\sigma^{2})$.
$X\sim\mathrm{NegativeBinomial}(r,p)$, where $\operatorname{\text{E}}[X]=r\frac{1-p}{p}$ and $\operatorname{\text{Var}}[X]=r\frac{1-p}{p^{2}}$. |  Convolution-closed in $r$ if $p$ is fixed.
$X\sim\mathrm{Gamma}(\alpha,\beta)$, where $\operatorname{\text{E}}[X]=\frac{\alpha}{\beta}$ and $\operatorname{\text{Var}}(X)=\frac{\alpha}{\beta^{2}}$. |  Convolution-closed in $\alpha$ if $\beta$ is fixed.
$X\sim\mathrm{Binomial}(r,p)$, where $\operatorname{\text{E}}[X]=rp$ and $\operatorname{\text{Var}}(X)=rp(1-p)$. |  Convolution-closed in $r$ if $p$ is fixed.
$X\sim\mathrm{InverseGaussian}(\mu w,\lambda w^{2})$ with $\operatorname{\text{E}}[X]=\mu w$ and $\operatorname{\text{Var}}(X)=\frac{w^{3}\mu^{3}}{w^{2}\lambda}=\frac{w\mu^{3}}{\lambda}$. |  Convolution-closed in $w$ if $\mu$ and $\lambda$ are fixed.
$X\sim\mathrm{GeneralizedPoisson}(\lambda,\theta)$, see Jørgensen [1992] for parameterization. |  Convolution-closed in $\lambda$ if $\theta$ is fixed.
$X\sim\mathrm{Tweedie}_{p}(\lambda,\theta),$ see Jørgensen [1992] for parameterization. |  Convolution-closed in $\lambda$ if $\theta$ and $p$ are fixed.
$\mathbb{X}\sim\mathrm{N}_{k}\left(\boldsymbol{\mu},\mathbb{\Sigma}\right)$, with $\operatorname{\text{E}}[\mathbb{X}]=\boldsymbol{\mu}$ and $\operatorname{\text{Var}}(\mathbb{X})=\mathbb{\Sigma}$. |  Convolution-closed in $(\boldsymbol{\mu},\mathbb{\Sigma})$.
$\mathbb{X}\sim\mathrm{Multinomial}_{k}\left(r,\mathbb{p}\right)$, with $\operatorname{\text{E}}[\mathbb{X}]=r\mathbb{p}$ and $\operatorname{\text{Var}}(\mathbb{X})=r\left(\mathrm{diag}(\mathbb{p})-\mathbb{p}\mathbb{p}^{T}\right)$. |  Convolution-closed in $r$ if $\mathbb{p}$ is fixed.
Table 1: A partial list of convolution-closed distributions. The last two rows
contain multivariate distributions. The results in each row are easily
verifiable. The generalized Poisson and Tweedie distributions are written in
their additive exponential dispersion family parameterization; see Jørgensen
[1992] or Jørgensen and Song [1998] for details of how these and other
distributions can be expressed in this way.

For a convolution-closed distribution $F_{\lambda}$, suppose that
$X^{\prime}\sim F_{\lambda_{1}}$ and $X^{\prime\prime}\sim F_{\lambda_{2}}$
with ${X}^{\prime}\perp\\!\\!\\!\perp{X}^{\prime\prime}$. Let
$G_{\lambda_{1},\lambda_{2},x}$ denote the conditional distribution of
$X^{\prime}\mid X^{\prime}+X^{\prime\prime}=x$. The density of the
distribution $G_{\lambda_{1},\lambda_{2},x}$ can be written down for any
$F_{\lambda}$ with a known density function [Jørgensen, 1992]. Furthermore, it
turns out that $G_{\lambda_{1},\lambda_{2},x}$ has a simple closed form for
several of the well-known distributions from Table 1; see Table 2. For
example, if $F_{\lambda}$ is the $\mathrm{Poisson}(\lambda)$ distribution,
then $G_{\lambda_{1},\lambda_{2},x}$ is the
$\mathrm{Binomial}\left(x,\frac{\lambda_{1}}{\lambda_{1}+\lambda_{2}}\right)$
distribution.

###  2.2 Data thinning

Recall from Section 2.1 that $G_{\lambda_{1},\lambda_{2},x}$ is the
conditional distribution of $X^{\prime}\mid X^{\prime}+X^{\prime\prime}=x$,
where ${X}^{\prime}\sim F_{\lambda_{1}}$ and ${X}^{\prime\prime}\sim
F_{\lambda_{2}}$ with ${X}^{\prime}\perp\\!\\!\\!\perp{X}^{\prime\prime}$. We
now introduce our proposal.

######  Algorithm 1 (Data thinning).

Suppose that we observe a realization $x$ of $X\sim F_{\lambda}$, where the
distribution $F_{\lambda}$ is convolution-closed in $\lambda$ with parameter
space $\Lambda$. For any value of $\epsilon\in(0,1)$ such that
$\epsilon\lambda\in\Lambda$ and $(1-\epsilon)\lambda\in\Lambda$:

  1. 1.

Draw $X^{(1)}\mid{X}=x\sim G_{\epsilon\lambda,(1-\epsilon)\lambda,x}$.

  2. 2.

Let $X^{(2)}=X-X^{(1)}$.

We now introduce our main theorem, which is closely related to a proposal by
Joe [1996] to construct autoregressive time series processes with known
marginal distributions.

######  Theorem 1.

Suppose that we apply Algorithm 1 to a realization $x$ of $X\sim F_{\lambda}$.
Then, the following results hold:

  1. 1.

$X^{(1)}\sim F_{\epsilon\lambda}$.

  2. 2.

$X^{(2)}\sim F_{(1-\epsilon)\lambda}$.

  3. 3.

$X^{(1)}\perp\\!\\!\\!\perp X^{(2)}$.

  4. 4.

If $F_{\lambda}$ satisfies the linear expectation property from Remark 2, then
$\operatorname{\text{E}}[X^{(1)}]=\epsilon\operatorname{\text{E}}[X]$ and
$\operatorname{\text{E}}[X^{(2)}]=(1-\epsilon)\operatorname{\text{E}}[X]$.

Theorem 1 is proven in Appendix A.1. The intuition for parts 1-3 is as
follows: if $X\sim F_{\lambda}$, then $X$ could have arisen as the sum of two
independent random variables $X^{\prime}\sim F_{\lambda_{1}}$ and
$X^{\prime\prime}\sim F_{\lambda_{2}}$, with
$\lambda_{1}+\lambda_{2}=\lambda$. Algorithm 1 works backwards to “undo” this
sum by generating random variables $X^{(1)}$ and $X^{(2)}$ that follow the
same distribution as $X^{\prime}$ and $X^{\prime\prime}$. Part 4 follows
directly from Remark 2. As we will see in Section 2.3, $\epsilon\in(0,1)$ is a
tuning parameter that governs a tradeoff between how much information is in
$X^{(1)}$ as opposed to $X^{(2)}$.

Theorem 1 guarantees that the decomposition provided by Algorithm 1 satisfies
the goals given in Section 1: namely (i) $X=X^{(1)}+X^{(2)}$; (ii)
$X^{(1)}\perp\\!\\!\\!\perp X^{(2)}$; and (iii) $X^{(1)}$ and $X^{(2)}$ follow
the same distribution as $X$, up to a (known) scaling of a parameter.

Table 2 summarizes the data thinning proposal for several well-known
distributions. The proposal in this paper extends well beyond this set of
distributions, although in some cases the conditional distribution
$G_{\lambda_{1},\lambda_{2},x}$ may not have a recognizable form.

######  Remark 3.

Some of the decompositions presented in Table 2 require knowledge of an
additional parameter that is not of primary interest. For example, thinning
the $\mathrm{N}(\mu,\sigma^{2})$ distribution requires knowledge of
$\sigma^{2}$. In fact, the same is true of the Gaussian decomposition given in
Leiner et al. [2022]; see Appendix C. In Section 2.4, we explore the
implications of performing data thinning in the presence of an unknown
nuisance parameter.

######  Remark 4.

Table 2 indicates that thinning the $\mathrm{Binomial}(r,p)$ distribution
requires that both $\epsilon r$ and $(1-\epsilon)r$ take on integer values; to
use the terminology of Joe [1996] and Jørgensen and Song [1998], the binomial
distribution is not “infinitely divisible”. This restriction becomes more
limiting in the extension to multiple folds given in Section 3, and prevents
us from thinning the Bernoulli distribution.

We now give an example of an application where data thinning is useful in
practice.

######  Application 1 (Model evaluation for unsupervised learning using data
thinning).

Suppose we observe ${X}_{ij}$ for $i=1,\ldots,n$ and $j=1,\ldots,d$, where
either each ${X}_{ij}$ is drawn independently from a univariate convolution-
closed distribution that satisfies the linear expectation property from Remark
2, or where each row $(X_{i1},\ldots,X_{id})^{T}$ is drawn independently from
a multivariate convolution-closed distribution that satisfies the linear
expectation property.

We wish to evaluate $\hat{\mu}({X})$ obtained from unsupervised learning (e.g.
clustering) on $X$, as an estimator for $\operatorname{\text{E}}[{X}]$.
Computing a loss function between $\hat{\mu}({X})$ and ${X}$ is
unsatisfactory, since the loss function will take on a small value if we
overfit the mean. Instead, we apply Algorithm 1 with $\epsilon\in(0,1)$ to
thin each element ${X}_{ij}$ into ${X}_{ij}^{(1)}$ and ${X}_{ij}^{(2)}$, and
then:

  1. 1.

Compute $\hat{\mu}({X}^{(1)})$ by applying unsupervised learning to $X^{(1)}$.
This provides an estimator of
$\operatorname{\text{E}}[{X}^{(1)}]=\epsilon\operatorname{\text{E}}[{X}]$
(Theorem 1, part 4).

  2. 2.

Compute a loss function between $\hat{\mu}({X}^{(1)})$ and ${X}^{(2)}$.

In Application 1, if $\epsilon=0.5$, then
$\operatorname{\text{E}}[{X}^{(1)}]=\operatorname{\text{E}}[{X}^{(2)}]=0.5\operatorname{\text{E}}[{X}]$.
Thus, $\hat{\mu}({X}^{(1)})$ is an estimator of
$\operatorname{\text{E}}[{X}^{(2)}]$, and so devising a suitable loss function
in Step 2 is straightforward. If $\epsilon\neq 0.5$, then
$\hat{\mu}({X}^{(1)})$ is a plug-in estimator of
$\frac{\epsilon}{1-\epsilon}\operatorname{\text{E}}[{X}^{(2)}]$ (Theorem 1,
part 4); we discuss this further in Section 2.3. The independence between
$X^{(1)}$ and $X^{(2)}$ prevents the loss function from taking on a small
value due to overfitting.

Distribution of ${X}$ |  Generate $X^{(1)}\mid{X}=x$ as: |  Dist. of $X^{(1)}$ |  Dist. of $X^{(2)}$ | Notes
---|---|---|---|---
$\text{Poisson}(\lambda)$ |  Draw $X^{(1)}\mid{X}=x\sim\text{Binomial}(x,\epsilon).$ | $\text{Poisson}(\epsilon\lambda)$ | $\text{Poisson}((1-\epsilon)\lambda)$ |
$\text{N}(\mu,\sigma^{2})$ |  Draw $X^{(1)}\mid{X}=x\sim\text{N}(\epsilon x,\epsilon(1-\epsilon)\sigma^{2}).$ | $\text{N}(\epsilon\mu,\epsilon\sigma^{2})$ | $\text{N}((1-\epsilon)\mu,(1-\epsilon)\sigma^{2})$ |  $\sigma^{2}$ must be known.
$\mathrm{NegBin}(r,p)$ |  Draw $X^{(1)}\mid{X}=x\sim\text{BetaBin}(x,\epsilon r,(1-\epsilon)r).$ | $\text{NegBin}(\epsilon r,p)$ | $\text{NegBin}((1-\epsilon)r,p)$ |  $r$ must be known.
$\text{Gamma}(\alpha,\beta)$ |  Draw $Z\sim\mathrm{Beta}\left(\epsilon\alpha,(1-\epsilon)\alpha\right)$, and let $X^{(1)}=x\cdot Z$. | $\text{Gamma}(\epsilon\alpha,\beta)$ | $\text{Gamma}((1-\epsilon)\alpha,\beta)$ |  $\alpha$ must be known.
$\text{Exp}(\lambda)=\text{Gamma}(1,\lambda)$ |  Draw $Z\sim\mathrm{Beta}(\epsilon,(1-\epsilon))$, and let $X^{(1)}=x\cdot Z$. | $\text{Gamma}(\epsilon,\lambda)$ | $\text{Gamma}((1-\epsilon),\lambda)$ |
$\text{Binomial}(r,p)$ |  Draw $X^{(1)}\mid{X}=x\sim\text{Hypergeometric}(\epsilon r,(1-\epsilon)r,x).$ | $\text{Binomial}(\epsilon r,p)$ | $\text{Binomial}((1-\epsilon)r,p)$ |  $r$ must be known. $\epsilon r$, $(1-\epsilon)r$ must be integers.
$\text{N}_{k}(\boldsymbol{\mu},\mathbb{\Sigma})$ |  Draw $\mathbb{X}^{(1)}\mid\mathbb{{X}}=\mathbb{x}\sim\text{N}(\epsilon\mathbb{x},\epsilon(1-\epsilon)\mathbb{\Sigma}).$ | $\text{N}_{k}(\epsilon\boldsymbol{\mu},\epsilon\mathbb{\Sigma})$ | $\text{N}_{k}((1-\epsilon)\boldsymbol{\mu},(1-\epsilon)\mathbb{\Sigma})$ |  $\Sigma$ must be known.
$\text{Multinom}_{k}(r,\mathbb{p})$ |  Draw $\mathbb{X}^{(1)}\mid\mathbb{{X}}=\mathbb{x}\sim$ | $\text{Multinom}_{k}(\epsilon r,\mathbb{p})$ | $\text{Multinom}_{k}((1-\epsilon)r,\mathbb{p})$ |  $r$ must be known.
|  $\text{MultivariateHypergeometric}(\mathbb{x}_{1},\mathbb{x}_{2},\ldots,\mathbb{x}_{k},\epsilon r)$. |  |  |  $\epsilon r$, $(1-\epsilon)r$ must be integers.
Table 2: Details of data thinning for several well-known distributions, using
the parameterizations given in Table 1. Note that in all cases,
$X^{(1)}\perp\\!\\!\\!\perp X^{(2)}$. While the exponential distribution
itself is not convolution-closed, writing the exponential distribution as a
special case of a Gamma distribution yields a useful decomposition. In the
last row, the parameterization of the multivariate hypergeometric distribution
is such that $\operatorname{\text{E}}\left[\mathbb{x}\right]=\epsilon
r\mathbb{x}/\sum_{j=1}^{k}\mathbb{x}_{j}$.

###  2.3 Role of the parameter $\epsilon$

In Algorithm 1, the parameter $\epsilon$ governs a trade-off between the
information in $X^{(1)}$ versus in $X^{(2)}$.

######  Example 2.1 (Fisher information in thinned Poisson distribution).

Let $X\sim\mathrm{Poisson}(\lambda)$. We thin $X$ to obtain
$X^{(1)}\sim\mathrm{Poisson}(\epsilon\lambda)$ and
$X^{(2)}\sim\mathrm{Poisson}((1-\epsilon)\lambda)$. Let $I_{X}(\lambda)$
denote the Fisher information contained in $X$ about the parameter $\lambda$.
Then $I_{X}(\lambda)=\frac{1}{\lambda}$, $I_{X^{(1)}}(\lambda)=\epsilon
I_{X}(\lambda)$, and $I_{X^{(2)}}(\lambda)=(1-\epsilon)I_{X}(\lambda)$.

######  Example 2.2 (Fisher information in thinned binomial distribution).

Let $X\sim\mathrm{Binomial}(r,p)$. We thin $X$ to obtain
$X^{(1)}\sim\mathrm{Binomial}(\epsilon r,p)$ and
$X^{(2)}\sim\mathrm{Binomial}((1-\epsilon)r,p)$. Let $I_{X}(p)$ denote the
Fisher information contained in $X$ about the parameter $p$. Then
$I_{X}(p)=\frac{r}{p(1-p)}$, $I_{X^{(1)}}(p)=\epsilon I_{X}(p)$, and
$I_{X^{(2)}}(p)=(1-\epsilon)I_{X}(p)$.

Similar results hold for other distributions in Table 2. Intuitively, as
$\epsilon$ increases, the amount of information in $X^{(1)}$ about the
parameter of interest increases, and the amount of information in $X^{(2)}$
decreases. This has implications for Application 1: as $\epsilon$ increases,
the quality of the estimator of the expected value increases (Step 1), but the
information available for computing the loss between this estimator and
$X^{(2)}$ decreases (Step 2).

In Section 2.2, we mentioned that the loss function in Step 2 of Application 1
must be chosen with care when $\epsilon\neq 0.5$. This can be seen in the
following example.

######  Example 2.3 (Application 1 with mean squared error loss).

Consider Step 2 of Application 1 with mean squared error (MSE) loss. Since
$\operatorname{\text{E}}[{X}^{(2)}]=\frac{1-\epsilon}{\epsilon}\operatorname{\text{E}}[{X}^{(1)}]$,
we define

$MSE=\frac{1}{nd}\left\|{X}^{(2)}-\frac{1-\epsilon}{\epsilon}\hat{\mu}({X}^{(1)})\right\|_{F}^{2},$

where the factor of $\frac{1-\epsilon}{\epsilon}$ turns an estimate of
$\operatorname{\text{E}}[{X}^{(1)}]$ into an estimate of
$\operatorname{\text{E}}[{X}^{(2)}]$.

As we will see in Section 5, we can also use alternative loss functions, such
as negative log likelihood, when considering Application 1.

Application 1 focuses on using data thinning to evaluate an estimator, and
Sections 5 and 6 of this paper focus on model selection and assessment.
However, data thinning could also be applied in a variety of other settings,
such as inference after variable selection in regression [Leiner et al.,
2022]. In these settings, the parameter $\epsilon$ still governs a trade-off
between the information available in $X^{(1)}$ (used for selection) and
$X^{(2)}$ (used for inference).

###  2.4 Effect of unknown nuisance parameters

For several of the distributions in Table 2, data thinning requires knowledge
of a nuisance parameter. For example, thinning a $\mathrm{N}(\mu,\sigma^{2})$
distribution requires knowledge of $\sigma^{2}$.

We now consider what happens when we perform data thinning on normally-
distributed data using an incorrect value of the variance. We refer to this
incorrect value as $\tilde{\sigma}^{2}$.

######  Proposition 1.

Suppose that we observe $x$ from $X\sim\mathrm{N}(\mu,\sigma^{2})$. We draw
$X^{(1)}\mid X=x\sim\mathrm{N}\left(\epsilon
x,\epsilon(1-\epsilon)\tilde{\sigma}^{2}\right)$, for some $\tilde{\sigma}$
that is not a function of $x$, and let $X^{(2)}=X-X^{(1)}$. Then:

  * (i) $X^{(1)}\sim\mathrm{N}\left(\epsilon\mu,\epsilon^{2}\sigma^{2}+\epsilon(1-\epsilon)\tilde{\sigma}^{2}\right)$,

  * (ii) $X^{(2)}\sim\mathrm{N}\left((1-\epsilon)\mu,(1-\epsilon)^{2}\sigma^{2}+\epsilon(1-\epsilon)\tilde{\sigma}^{2}\right)$, and

  * (iii) $\mathrm{Cov}\left(X^{(1)},X^{(2)}\right)=\epsilon(1-\epsilon)\left(\sigma^{2}-\tilde{\sigma}^{2}\right)$.

Proposition 2(iii) indicates that if we apply data thinning with _not enough_
noise ($\tilde{\sigma}^{2}<\sigma^{2}$), then $X^{(1)}$ and $X^{(2)}$ are
positively correlated. On the other hand, if we apply data thinning with _too
much_ noise ($\tilde{\sigma}^{2}>\sigma^{2}$), then $X^{(1)}$ and $X^{(2)}$
are negatively correlated. Similar results hold for the negative binomial
distribution and the Gamma distribution.

######  Proposition 2.

Suppose that we observe $x$ from $X\sim\mathrm{NegBin}(r,p)$. We draw
$X^{(1)}\mid
X=x\sim\mathrm{BetaBin}\left(x,\epsilon\tilde{r},(1-\epsilon)\tilde{r}\right)$
for some $\tilde{r}$ that is not a function of $x$, and let
$X^{(2)}=X-X^{(1)}$. Then
$\mathrm{Cov}\left(X^{(1)},X^{(2)}\right)=\epsilon(1-\epsilon)r\left(\frac{1-p}{p}\right)^{2}\left(1-\frac{r+1}{\tilde{r}+1}\right).$

######  Proposition 3.

Suppose that we observe $x$ from $X\sim\mathrm{Gamma}(\alpha,\beta)$. We let
$X^{(1)}=x\times Z$, where
$Z\sim\mathrm{Beta}\left(\epsilon\tilde{\alpha},(1-\epsilon)\tilde{\alpha}\right)$
for some $\tilde{\alpha}$ that is not a function of $x$. We let
$X^{(2)}=X-X^{(1)}$. Then
$\mathrm{Cov}\left(X^{(1)},X^{(2)}\right)=\epsilon(1-\epsilon)\frac{\alpha}{\beta^{2}}\left(1-\frac{\alpha+1}{\tilde{\alpha}+1}\right).$

Propositions 1–3 are proven in Appendix A.2. Figure 1 verifies these results
empirically.

The results in this section assume that $\tilde{\sigma}$, $\tilde{b}$, and
$\tilde{\alpha}$ are not a function of $X$. In practice, one might estimate
the unknown parameters $\sigma$, $b$, and $\alpha$ using additional data.

Figure 1: _Left:_ We generate 100,000 realizations of $X\sim\mathrm{N}(7,5)$.
For 50 values of $\tilde{\sigma}^{2}$, we thin $X$ into $X^{(1)}$ and
$X^{(2)}$ using $\tilde{\sigma}^{2}$ instead of $\sigma^{2}=5$. _Center:_ We
generate 100,000 realizations of $X\sim\mathrm{NB}(7,0.7)$. For 50 values of
$\tilde{r}$, we thin $X$ into $X^{(1)}$ and $X^{(2)}$ using $\tilde{r}$
instead of $r=7$. _Right:_ We generate 100,000 realizations of
$X\sim\mathrm{Gamma}(7,5)$. For 50 values of $\tilde{\alpha}$, we thin $X$
into $X^{(1)}$ and $X^{(2)}$ using $\tilde{\alpha}$ instead of $\alpha=7$.
_All:_ In each panel, for each value of the nuisance parameter, we display the
empirical correlation between $X^{(1)}$ and $X^{(2)}$, along with the
theoretical correlation suggested by Propositions 1–3.

##  3 Multi-fold data thinning

Data thinning involves decomposing $X$ into $X^{(1)}$ and $X^{(2)}$, which
each have the same distribution as $X$ (up to a known parameter scaling). It
can be applied recursively to create $M$ independent data folds,
$X^{(1)},\ldots,X^{(M)}$, that sum to $X$, as in the following example.

######  Example 3.1 (Recursive thinning of the normal distribution).

Let $x$ denote a realization of $X\sim N(\mu,\sigma^{2})$. Given
$\epsilon_{1},\epsilon_{2},\epsilon_{3}\in(0,1)$ with
$\epsilon_{1}+\epsilon_{2}+\epsilon_{3}=1$, we proceed as follows:

  1. 1.

Draw $X^{(1)}\mid
X\sim\mathrm{N}\left(\epsilon_{1}X,\epsilon_{1}(1-\epsilon_{1})\sigma^{2}\right)$.
Let $X^{(2,3)}=X-X^{(1)}$. By Theorem 1, $X^{(1)}\sim
N\left(\epsilon_{1}\mu,\epsilon_{1}\sigma^{2}\right)$, $X^{(2,3)}\sim
N\left((1-\epsilon_{1})\mu,(1-\epsilon_{1})\sigma^{2}\right)$, and
$X^{(1)}\perp\\!\\!\\!\perp X^{(2,3)}$.

  2. 2.

Draw $X^{(2)}\mid
X^{(2,3)}\sim\mathrm{N}\left(\frac{\epsilon_{2}}{1-\epsilon_{1}}X^{(2,3)},\frac{\epsilon_{2}}{1-\epsilon_{1}}(1-\frac{\epsilon_{2}}{1-\epsilon_{1}})(1-\epsilon_{1})\sigma^{2}\right)$.
Let $X^{(3)}=X-X^{(1)}-X^{(2)}$. By Theorem 1, $X^{(2)}\sim
N(\epsilon_{2}\mu,\epsilon_{2}\sigma^{2})$ and $X^{(3)}\sim
N(\epsilon_{3}\mu,\epsilon_{3}\sigma^{2})$ with $X^{(2)}\perp\\!\\!\\!\perp
X^{(3)}$. Furthermore, since $\left(X^{(2)},X^{(3)}\right)$ is a function of
$X^{(2,3)}$, each of $X^{(2)}$ and $X^{(3)}$ remain independent of $X^{(1)}$.

Thus, $X^{(m)}\sim N(\epsilon_{m}\mu,\epsilon_{m}\sigma^{2})$ for
$m=1,\ldots,3$, $X^{(1)}+X^{(2)}+X^{(3)}=X$, and $X^{(m)}\perp\\!\\!\\!\perp
X^{(m^{\prime})}$ for all $m\neq m^{\prime}$.

While Example 3.1 can be extended to create $M>3$ folds, this recursive
approach can be cumbersome. In Example 1.3 of Section 1, we saw that in the
case of the Gamma distribution, there is a simpler way to create multiple
folds without recursion. We will now provide a general form of this result.
This requires additional notation: we let
$G_{\lambda_{1},\lambda_{2},\ldots,\lambda_{M},x}$ denote the joint
distribution of
$\left(X_{1},\ldots,X_{M}\right)\mid{X_{1}+X_{2}+\ldots+X_{M}}=x$, where
$X_{m}\overset{\mathrm{ind.}}{\sim}F_{\lambda_{m}}$, for $m=1,\ldots,M$, and
where $F_{\lambda}$ is a convolution-closed distribution. The following
algorithm and theorem mimic Algorithm 1 and Theorem 1.

######  Algorithm 2 (Multi-fold data thinning).

Suppose that we observe a realization $x$ of $X\sim F_{\lambda}$, where
$F_{\lambda}$ is a convolution-closed distribution with parameter space
$\Lambda$. Choose $\epsilon_{1},\ldots,\epsilon_{M}\in(0,1)$ such that
$\sum_{m=1}^{M}\epsilon_{m}=1$ and $\epsilon_{m}\lambda\in\Lambda$ for
$m=1,\ldots M$. Then, draw $\left(X^{(1)},\ldots,X^{(M)}\right)\sim
G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\ldots,\epsilon_{M}\lambda,x}$.

######  Theorem 2.

Suppose that we apply Algorithm 2 to a realization $x$ of $X\sim F_{\lambda}$,
for a convolution-closed distribution $F_{\lambda}$. Then, the following
results hold:

  1. 1.

$X^{(m)}\sim F_{\epsilon_{m}\lambda}$ for $m=1,\ldots,M$.

  2. 2.

$X^{(m)}\perp\\!\\!\\!\perp X^{(m^{\prime})}$ for all $m\neq m^{\prime}$.

  3. 3.

$X^{(1)}+X^{(2)}+\ldots+X^{(M)}=X$.

  4. 4.

If $F_{\lambda}$ satisfies the linear expectation property from Remark 2, then
$\operatorname{\text{E}}[X^{(m)}]=\epsilon_{m}\operatorname{\text{E}}[X]$ for
$m=1,\ldots,M$.

The proof of Theorem 2 is included in Appendix A.1, and is a straightforward
extension of that of Theorem 1. The intuition for parts 1-3 is as follows: we
know that $X\sim F_{\lambda}$ could have arisen as the sum of $M$ independent
random variables $X_{1},\ldots,X_{M}$ such that $X_{m}\sim
F_{\epsilon_{m}\lambda}$. If we draw
$\left(X^{(1)},\ldots,X^{(M)}\right)|X=x\sim
G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\ldots,\epsilon_{M}\lambda,x}$,
then the joint distribution of $\left(X^{(1)},\ldots,X^{(M)}\right)$ equals
the joint distribution of $\left(X_{1},\ldots,X_{M}\right)$, i.e. it is the
joint distribution of $M$ independent random variables with distributions
$F_{\epsilon_{1}\lambda},\ldots,F_{\epsilon_{M}\lambda}$. Part 4 follows
directly from Remark 2.

We now revisit the case of the normal distribution from Example 3.1.

######  Example 3.2 (Multi-fold thinning of the normal distribution).

Let $X\sim N(\mu,\sigma^{2})$. To generate $M=3$ independent folds that
contain proportions $\epsilon_{1},\epsilon_{2}$, and
$\epsilon_{3}=1-\epsilon_{1}-\epsilon_{2}$ of the data, we draw

$\small\left[\begin{matrix}X^{(1)}\\\ X^{(2)}\\\
X^{(3)}\end{matrix}\right]\mid X=x\sim
N\left(\left[\begin{matrix}\epsilon_{1}x\\\ \epsilon_{2}x\\\
\epsilon_{3}x\end{matrix}\right],\left[\begin{matrix}\epsilon_{1}(1-\epsilon_{1})\sigma^{2}&-\epsilon_{1}\epsilon_{2}\sigma^{2}&-\epsilon_{1}\epsilon_{3}\sigma^{2}\\\
-\epsilon_{1}\epsilon_{2}\sigma^{2}&\epsilon_{2}(1-\epsilon_{2})\sigma^{2}&-\epsilon_{2}\epsilon_{3}\sigma^{2}\\\
-\epsilon_{1}\epsilon_{3}\sigma^{2}&-\epsilon_{2}\epsilon_{3}\sigma^{2}&\epsilon_{3}(1-\epsilon_{3})\sigma^{2}\end{matrix}\right]\right).$

One can verify that this multivariate normal corresponds to
$G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\epsilon_{3}\lambda,x}$. By
Theorem 2, $X^{(1)},X^{(2)},$ and $X^{(3)}$ are independent and $X^{(m)}\sim
N(\epsilon_{m}\mu,\epsilon_{m}\sigma^{2})$ for $m=1,2,3$. We note that this
distribution
$G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\epsilon_{3}\lambda,x}$ is a
degenerate multivariate normal distribution, which enforces the constraint
that the realized values of $X^{(1)},X^{(2)}$, and $X^{(3)}$ sum to $x$.

Table 3 reveals that
$G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\ldots,\epsilon_{M}\lambda,x}$ in
Algorithm 2 has a very simple form for every univariate distribution in Table
2. We omit the multivariate distributions to avoid cumbersome notation.

Distribution of ${X}$ |  Generate $\left(X^{(1)},\ldots,X^{(M)}\right)\mid{X}=X$ as: |  Dist. of $X^{(m)}$ | Notes
---|---|---|---
$\text{Poisson}(\lambda)$ |  Draw $\left(X^{(1)},\ldots,X^{(M)}\right)\mid{X}=x\sim\text{Multinomial}(x,\epsilon_{1},\ldots,\epsilon_{M})$. | $\text{Poisson}(\epsilon_{m}\lambda)$ |
$\text{N}(\mu,\sigma^{2})$ |  Draw $\left(X^{(1)},\ldots,X^{(M)}\right)\mid{X}=x\sim\text{N}\left(\mu\boldsymbol{\epsilon},\sigma^{2}\mathrm{diag}(\boldsymbol{\epsilon})-\sigma^{2}\boldsymbol{\epsilon}\boldsymbol{\epsilon}^{T}\right)$, |  $\text{N}(\epsilon_{m}\mu,\epsilon_{m}\sigma^{2})$, |  $\sigma^{2}$ must be known.
|  where $\boldsymbol{\epsilon}=(\epsilon_{1},\ldots,\epsilon_{M})^{T}$. Note that this multivariate normal distribution |  |
| is degenerate (the covariance matrix is not invertible). |  |
$\mathrm{NegBin}(r,p)$ |  Draw $\left(X^{(1)},\ldots,X^{(M)}\right)\mid{X}=x\sim\text{DirichletMultinomial}(X,\epsilon_{1}r,\ldots,\epsilon_{M}r)$. | $\text{NegBin}(\epsilon_{m}r,p)$ |  $r$ must be known.
$\text{Gamma}(\alpha,\beta)$ |  Draw $Z\sim\text{Dirichlet}\left(\epsilon_{1}\alpha,\ldots,\epsilon_{M}\alpha\right)$, and and let $\left(X^{(1)},\ldots,X^{(M)}\right)=x\cdot Z$ | $\text{Gamma}(\epsilon_{m}\alpha,\beta)$ |  $\alpha$ must be known.
$\text{Exp}(\lambda)=\text{Gamma}(1,\lambda)$ |  Draw $Z\sim\mathrm{Dirichlet}(\epsilon_{1},\ldots,\epsilon_{M})$, and let $\left(X^{(1)},\ldots,X^{(M)}\right)=x\cdot Z$. | $\text{Gamma}(\epsilon_{m},\lambda)$ |
$\text{Binomial}(r,p)$ |  Draw $\left(X^{(1)},\ldots,X^{(M)}\right)\mid{X}=x\sim\text{MultivariateHypergeometric}(\epsilon_{1}r,\ldots,\epsilon_{M}r,x)$. | $\text{Binomial}(\epsilon_{m}r,p)$ |  $r$ must be known.
|  |  |  each $\epsilon_{m}r$ must be an integer.
Table 3: Details of how to perform multi-fold data thinning using Algorithm 2
for several common univariate distributions. Each row can be verified using
standard properties of these distributions.

The parameters $\epsilon_{1},\ldots,\epsilon_{M}$ in Algorithm 2 govern the
same information tradeoff that was described in Section 2.3. Folds with the
largest values of $\epsilon_{m}$ contain the most information about the
parameter of interest. We now consider the following extension of Application
1.

######  Application 2 (Cross validation for unsupervised learning using multi-
fold thinning).

We observe ${X}_{ij}$ for $i=1,\ldots,n$ and $j=1,\ldots,d$, where the
elements ${X}_{ij}$ are drawn independently from a univariate convolution-
closed family that satisfies the linear expectation property from Remark 2.
Alternatively, we can allow the rows $(X_{i1},\ldots,X_{id})^{T}$ to be drawn
independently from a multivariate convolution-closed family that satisfies the
linear expectation property. We wish to evaluate $\hat{\mu}({X})$ obtained
through unsupervised learning (e.g. clustering) as an estimator for
$\operatorname{\text{E}}[{X}]$. We apply Algorithm 2 with parameters
$\epsilon_{1},\ldots,\epsilon_{M}$ to thin each element ${X}_{ij}$ into
${X}_{ij}^{(1)},\ldots,{X}_{ij}^{(M)}$, and proceed as follows:

  1. 1.

For $m=1,\ldots,M$:

    1. (a)

Define $X^{(-m)}:=X-X^{(m)}$. Apply unsupervised learning to ${X}^{(-m)}$ to
obtain $\hat{\mu}\left({X}^{(-m)}\right)$, which is an estimator of
$\operatorname{\text{E}}[{X}^{(-m)}]=(1-\epsilon_{m})\operatorname{\text{E}}[{X}]$.

    2. (b)

Compute a loss function between $\hat{\mu}\left({X}^{(-m)}\right)$ and
${X}^{(m)}$. For example, as in Example 2.3, compute the mean squared error
between $\frac{\epsilon_{m}}{1-\epsilon_{m}}\hat{\mu}({X}^{(-m)})$ and
$X^{(m)}$.

  2. 2.

Evaluate the estimator $\hat{\mu}(\cdot)$ by aggregating (e.g. averaging) the
loss across folds.

For simplicity, we suggest setting $\epsilon_{1}=\ldots=\epsilon_{M}=1/M$. The
advantage of Application 2 over Application 1 with $\epsilon=1/M$ is reduction
of the variance of the loss function via averaging. We will demonstrate the
practical advantages of multi-fold thinning in Section 5.

##  4 Comparison of data thinning and data fission

As mentioned in Section 1, Leiner et al. [2022]’s data fission proposal
provides an alternate set of strategies to decompose $X$ into $X^{(1)}$ and
$X^{(2)}$. In general, the data fission decompositions provided by Leiner et
al. [2022] face the following drawbacks:

  1. 1.

_With the exception of the normal and Poisson decompositions, $X^{(1)}$ and
$X^{(2)}$ are not independent. _ In the Poisson case, the data thinning and
data fission proposals are identical [and have also been explored by a number
of authors previously; see Neufeld et al., 2022, Oliveira et al., 2022, Chen
et al., 2021, Gerard, 2020, among others]. In Appendix C, we compare the data
thinning and data fission decompositions of the normal distribution. We note
that the data fission decomposition of the normal distribution has also been
explored by a number of authors [Tian, 2020, Oliveira et al., 2021, Rasines
and Young, 2021].

  2. 2.

_The distribution of $X^{(1)}$, and the conditional distribution of
$X^{(2)}|X^{(1)}$, need not resemble the distribution of $X$._ For instance,
in the Gamma decomposition in Example 1.1, $X^{(1)}$ is a set of $B$ negative
binomial random variables.

  3. 3.

_The parameters of interest are entangled in the conditional distribution of
$X^{(2)}\mid X^{(1)}$_. We see this in the Gamma decomposition in Example 1.1,
and will additionally see it in the negative binomial example below (Example
4.1). This entanglement poses challenges in applying data fission.

  4. 4.

_The tuning parameter that governs the information trade-off between $X^{(1)}$
and $X^{(2)}$ can be hard to interpret_. For instance, in the case of the
Gamma decomposition in Example 1.1, the tuning parameter is
$B\in\\{1,2,\ldots\\}$. In the case of the normal decomposition (described in
Appendix C), the tuning parameter is $\tau>0$. And in the case of the binomial
decomposition in Example 4.1 below, it is $\epsilon\in(0,1)$.

  5. 5.

_The proposal has no clear extension to the case of multiple folds._

By contrast, data thinning provides a simple recipe to decompose convolution-
closed distributions — a large class that includes examples not considered by
Leiner et al. [2022], such as the multinomial, the Tweedie, or the Inverse
Gaussian (Table 1) — into two or more _independent_ random variables _that
follow the same distribution as $X$, up to scaling by a known constant_.
Furthermore, in this recipe, the tuning parameter $\epsilon$ always lies in
$(0,1)$, and can always be interpreted as governing a simple trade-off between
the information in the thinned random variables.

We now revisit Application 1 to see a concrete example in which data thinning
is straightforward but the proposal of Leiner et al. [2022] is difficult to
use in practice.

######  Example 4.1 (A comparison of negative binomial decompositions).

We observe
${X}_{ij}\sim\mathrm{NegativeBinomial}\left({r}_{ij},{p}_{ij}\right)$ for
$i=1,\ldots,n$ and $j=1,\ldots,p$.

From Table 2, data thinning requires ${r}_{ij}$ to be known, and yields
${X}^{(1)}_{ij}\sim\mathrm{NegativeBinomial}\left(\epsilon{r}_{ij},{p}_{ij}\right)$
and
${X}^{(2)}_{ij}\sim\mathrm{NegativeBinomial}\left((1-\epsilon){r}_{ij},{p}_{ij}\right)$,
with ${X}^{(1)}\perp\\!\\!\\!\perp{X}^{(2)}$,
$\operatorname{\text{E}}[{X}^{(1)}]=\epsilon\operatorname{\text{E}}[{X}]$, and
$\operatorname{\text{E}}[{X}^{(2)}]=(1-\epsilon)\operatorname{\text{E}}[{X}]$.
Thus, as in Application 1 and Example 2.3, $\hat{\mu}\left({X}^{(1)}\right)$
is an estimator for $\epsilon\operatorname{\text{E}}[{X}]$, and we can
evaluate $\hat{\mu}({X}^{(1)})$ by computing the mean squared error between
${X}^{(2)}$ and $\frac{1-\epsilon}{\epsilon}\hat{\mu}\left({X}^{(1)}\right)$.

For $\epsilon\in(0,1)$, the data fission proposal of Leiner et al. [2022]
draws ${X}^{(1)}_{ij}\mid{X}_{ij}\sim\mathrm{Binomial}({X}_{ij},\epsilon)$ and
sets ${X}^{(2)}={X}-{X}^{(1)}$. Under this decomposition,
${X}^{(1)}_{ij}\sim\mathrm{NegativeBinomial}\left({r}_{ij},\frac{{p}_{ij}}{{p}_{ij}+\epsilon(1-{p}_{ij})}\right)$,
and so
$\operatorname{\text{E}}[{X}^{(1)}]=\epsilon\operatorname{\text{E}}[{X}]$.
Moreover,
${X}^{(2)}_{ij}\mid{X}^{(1)}_{ij}\sim\mathrm{NegativeBinomial}\left({r}_{ij}+{X}^{(1)}_{ij},{p}_{ij}+\epsilon-{p}_{ij}\epsilon\right)$,
and so $\operatorname{\text{E}}[{X}^{(2)}_{ij}\mid{X}^{(1)}_{ij}]=$
$\left({r}_{ij}+{X}^{(1)}_{ij}\right)\left(\frac{1-{p}_{ij}-\epsilon+\epsilon{p}_{ij}}{{p}_{ij}+\epsilon-\epsilon{p}_{ij}}\right)$.
As ${X}^{(1)}$ and ${X}^{(2)}$ are not independent, we cannot simply use mean
squared error loss between ${X}^{(2)}$ and
$\frac{1-\epsilon}{\epsilon}\hat{\mu}\left({X}^{(1)}\right)$ to evaluate the
estimator.

At a glance, an advantage of data fission is that it does not require
knowledge of ${r}_{ij}$ to obtain ${X}_{ij}^{(1)}$ and ${X}_{ij}^{(2)}$.
However, looking at the conditional distribution of ${X}^{(2)}\mid{X}^{(1)}$,
we see that $E[{X}^{(2)}]$ is a function of $\operatorname{\text{E}}[{X}]$
that is tractable only if ${r}_{ij}$ is known. Thus, in practice, using the
data fission approach to do inference on $\operatorname{\text{E}}[{X}]$ will
require knowing or accurately estimating the nuisance parameters ${r}_{ij}$.

Similar issues arise for other decompositions given in Leiner et al. [2022].

##  5 Simulation Study

The data thinning proposal can be applied in a variety of settings, including
the model selection and validation problems considered in Leiner et al.
[2022]. Here, we focus instead on the application of data thinning to cross-
validation for unsupervised learning. Data thinning is particularly attractive
in this setting, as a traditional sample splitting approach is not suitable;
see Owen and Perry [2009] or Fu and Perry [2020] for further discussion. In
what follows, we apply data thinning and multi-fold thinning to unsupervised
learning problems, and contrast their performance against naive approaches
that use the same data to both fit and validate the unsupervised models.
Specifically, we consider the problems of selecting the number of principal
components (PCs) in binomial distributed data, and selecting the number of
clusters in $K$-means clustering on Gamma distributed data.

###  5.1 Simulation setup

####  5.1.1 Choosing the number of PCs on binomial data

We generate binomial-distributed data with $n=250$ observations and $d=100$
dimensions. Specifically, for $i=1,\dots,n$ and $j=1,\dots,d$, we generate
$X_{ij}\overset{\mathrm{ind.}}{\sim}\text{Binomial}(r,p_{ij})$ where $r=100$
and $p$ is an unknown $n\times d$ matrix of probabilities. We construct
$\text{logit}(p)$ as a rank-$K^{*}=10$ matrix with singular values
$5,6,\dots,14$. Additional details are provided in Appendix D.1. In this
simulation, our goal is to estimate $K^{*}$.

####  5.1.2 Choosing the number of clusters on Gamma data

We generate datasets $X\in\mathbb{R}^{n\times d}$ such that there are 100
observations in each of $K^{*}$ clusters, for a total of $n=100K^{*}$
observations. Our objective is to estimate $K^{*}$.

We let
$X_{ij}\overset{\mathrm{ind.}}{\sim}\text{Gamma}(\lambda,\theta_{c_{i},j})$,
for $i=1,\ldots,n$ and $j=1,\ldots,d$, where $c_{i}\in\\{1,2,\dots,K^{*}\\}$
indexes the true cluster membership of the $i$th observation. The shape
parameter $\lambda$ is a known constant common across all clusters and all
dimensions, whereas the rate parameter $\theta$ is an unknown $K^{*}\times d$
matrix such that each cluster has its own $d$-dimensional rate parameter. We
generate data under two regimes:

  * _Small $d$, small $K^{*}$:_ We take $d=2$ and $K^{*}=4$.

  * _Large $d$, large $K^{*}$:_ We take $d=100$ and $K^{*}=10$.

The values of $\lambda$ and $\theta$ are provided in Appendix D.1. A sample
“small $d$, small $K^{*}$” dataset is presented in Figure 2, alongside the
output of data thinning with $\epsilon=0.5$.

Figure 2: Left: A simulated dataset in the “small $d$, small $K^{*}$” setting
described in Section 5.1.2. Center/Right: The result of data thinning with
$\epsilon=0.5$. Note that the rate parameter has not changed, but the shape
parameter has decreased by a factor of two.

###  5.2 Methods

####  5.2.1 Selecting the number of principal components in binomial data

The goal of the following algorithm is to quantify the negative log-likelihood
(NLL) associated with approximating a binomial matrix, $X^{(\mathrm{test})}$,
with the first $K$ PCs of a binomial matrix, $X^{(\mathrm{train})}$, in the
setting described in Section 5.1.1.

######  Algorithm 3 (Evaluating binomial PCA with negative log-likelihood
loss).

.
Input: A positive integer $K$, and matrices
$X^{(\mathrm{train})}\in\mathbb{R}^{n\times d}$ and
$X^{(\mathrm{test})}\in\mathbb{R}^{n\times d}$, where
$X^{(\mathrm{train})}_{ij}\overset{\textup{ind}}{\sim}\textup{Binomial}\left(a^{\textup{(train)}}r,p_{ij}\right)$,
$X^{(\mathrm{test})}_{ij}\overset{\textup{ind}}{\sim}\textup{Binomial}\left(a^{\textup{(test)}}r,p_{ij}\right)$,
$r$ is the known number of trials, and $a^{\textup{(train)}}$ and
$a^{\textup{(test)}}$ are known constants in $(0,1]$ such that
$a^{\textup{(train)}}r$ and $a^{\textup{(test)}}r$ are integers.

  1. 1.

Compute the singular value decomposition of the observed log-odds matrix,
$\textup{logit}\left(\frac{X^{(\mathrm{train})}_{ij}+0.001}{a^{\textup{(train)}}r+0.002}\right)=\hat{U}\hat{D}\hat{V}^{T}$.
Pseudo-counts prevent taking the logit of 0 or 1.

  2. 2.

Construct the rank-$K$ approximation,
$p^{(K)}=\textup{expit}\left(\hat{U}_{1:K}\hat{D}_{1:K}\hat{V}_{1:K}\right)$.

  3. 3.

Compute the loss function:

NLL $\displaystyle=-\sum_{i,j}\log
f\left(X^{(\mathrm{test})}_{ij}\Big{|}a^{\textup{(test)}}r,p^{(K)}_{ij}\right)$
(1) $\displaystyle=-\sum_{i,j}\left[\log{a^{\textup{(test)}}r\choose
X^{(\mathrm{test})}_{ij}}+X^{(\mathrm{test})}_{ij}\log
p^{(K)}_{ij}+\left(a^{\textup{(test)}}r-X^{(\mathrm{test})}_{ij}\right)\log\left(1-p^{(K)}_{ij}\right)\right],$

where $f(\cdot\mid,r,p)$ is the density function for the
$\mathrm{Binomial}(r,p)$ distribution.

####  5.2.2 Choosing the number of clusters on Gamma data

The goal of the following algorithm is to quantify the NLL that results from
approximating a Gamma-distributed matrix, $X^{(\mathrm{test})}$, using
clusters obtained from a Gamma-distributed matrix, $X^{(\mathrm{train})}$, in
the setting described in Section 5.1.2.

######  Algorithm 4 (Evaluating Gamma clusters with negative log-likelihood
loss).

.
Input: A positive integer $K$, and matrices
$X^{(\mathrm{train})}\in\mathbb{R}^{n\times d}$ and
$X^{(\mathrm{test})}\in\mathbb{R}^{n\times d}$, where
$X^{(\mathrm{train})}_{ij}\overset{\textup{ind}}{\sim}\textup{Gamma}\left(a^{\textup{(train)}}\lambda,\theta_{c_{i},j}\right)$,
$X^{(\mathrm{test})}_{ij}\overset{\textup{ind}}{\sim}\textup{Gamma}\left(a^{\textup{(test)}}\lambda,\theta_{c_{i},j}\right)$.
Here, $\theta\in(0,\infty)^{K^{*}\times d}$ where $\theta_{c_{i},j}$ is the
true but unknown rate parameter for the $c_{i}$th cluster in the $j$th
dimension, $c_{i}\in\\{1,2,\dots,K^{*}\\}$, $\lambda$ is the known shape
parameter, and $a^{\textup{(train)}}$ and $a^{\textup{(test)}}$ are known
constants in $(0,1]$.

  1. 1.

Run $K$-means on $X^{(\mathrm{train})}$ to estimate $K$ clusters. Denote the
cluster assignment of the $i$th observation as $\hat{c}_{i}$.

  2. 2.

Within each cluster, estimate the parameters using the closed-form estimators
of Ye and Chen [2017] and Louzada et al. [2019]. Let $\hat{\lambda}^{(K)}$ and
$\hat{\theta}^{(K)}$ denote the $K\times d$ estimated parameter matrices.

  3. 3.

Compute the loss function:

NLL $\displaystyle=-\sum_{i,j}\log
f\left(X^{(\mathrm{test})}_{ij}\Big{|}\frac{a^{\textup{(test)}}}{a^{\textup{(train)}}}\hat{\lambda}^{(K)}_{\hat{c}_{i},j},\hat{\theta}^{(K)}_{\hat{c}_{i},j}\right)$
(2)
$\displaystyle=-\sum_{i,j}\left[\frac{a^{\textup{(test)}}}{a^{\textup{(train)}}}\hat{\lambda}^{(K)}_{\hat{c}_{i},j}\log\hat{\theta}^{(K)}_{\hat{c}_{i},j}-\log\Gamma\left(\frac{a^{\textup{(test)}}}{a^{\textup{(train)}}}\hat{\lambda}^{(K)}_{\hat{c}_{i},j}\right)\right.$
$\displaystyle\quad\left.+\left(\frac{a^{\textup{(test)}}}{a^{\textup{(train)}}}\hat{\lambda}^{(K)}_{\hat{c}_{i},j}-1\right)\log
X^{(\mathrm{test})}_{ij}-\hat{\theta}^{(K)}_{\hat{c}_{i},j}X^{(\mathrm{test})}_{ij}\right],$

where $f(\cdot\mid\lambda,\theta)$ is the density function for the
$\mathrm{Gamma}(\lambda,\theta)$ distribution.

####  5.2.3 Methods for comparison

For a range of values of $K$, we will apply Algorithms 3 and 4 as follows:

  * Naive method that re-uses data: We apply Algorithms 3 and 4 with $X^{(\mathrm{train})}=X^{(\mathrm{test})}=X$ and $a^{\textup{(train)}}=a^{\textup{(test)}}=1$. That is, in the binomial PCA setting, we compute the first K PCs of the data $X$, and then compute the negative log-likelihood of this rank-$K$ approximation of $X$. In the Gamma clustering setting, we cluster $X$, and then compute the negative log-likelihood of the estimated clusters.

  * Data thinning: We first thin the data as in Algorithm 1 with $\epsilon=0.5$ or $0.8$, to obtain $X^{(1)}$ and $X^{(2)}$. We then apply Algorithms 3 and 4 with $X^{(\mathrm{train})}=X^{(1)}$, $X^{(\mathrm{test})}=X^{(2)}$, $a^{\textup{(train)}}=\epsilon$, and $a^{\textup{(test)}}=1-\epsilon$.

  * Multi-fold thinning: We first thin the data into $M=5$ folds using Algorithm 2 with $\epsilon_{1}=\dots=\epsilon_{M}=0.2$. Then, for $m=1,\dots,M$, we sum the data from all but the $m$th fold into $X^{(-m)}$, and apply Algorithms 3 and 4 with $X^{(\mathrm{train})}=X^{(-m)}$, $X^{(\mathrm{test})}=X^{(m)}$, $a^{\textup{(train)}}=\frac{M-1}{M}$, and $a^{\textup{(test)}}=\frac{1}{M}$. Finally, we average the loss functions in Equations 1 and 2 across the $M$ folds.

Before examining the results of the simulations, we comment on the methods
described above. Our goal is to select the value of $K$ that minimizes the
NLL. Because data thinning produces independent training and test sets, we
expect that the data thinning approaches will produce “U”-shaped NLL curves,
as a function of $K$. By contrast, in the naive approach, the full data $X$ is
used to fit the model and to compute the NLL in Equations 1 and 2, resulting
in monotonically decreasing loss curves, as a function of $K$.

We have chosen to use an NLL loss in both Algorithms 3 and 4, but other loss
functions could also be used and would yield similar results. Appendix D.2
extends Algorithms 3 and 4 to the case of mean squared error loss.

###  5.3 Results

Figure 3 displays the NLL values for all three simulation settings as a
function of $K$; results have been averaged over $2,000$ simulated datasets
and rescaled to the $[0,1]$ interval for ease of comparison. The values of $K$
with the lowest average NLL are circled on the plots. As expected, the data
thinning approaches in Figure 3 exhibit sharp minimum values, as opposed to
the monotonically decreasing curves produced by the naive method. The data
thinning approaches correctly select the true value of $K=K^{*}$ in all three
settings, except for data thinning with $\epsilon=0.5$ in the binomial PCA
setting. In that case, the low value of $\epsilon$ allocates too much
information to the test set, resulting in inadequate signal from the weakest
principal components in the training set. Selecting a larger value of
$\epsilon$ remedies this issue, as seen with $\epsilon=0.8$.

Figure 3: The negative log-likelihood loss averaged over 2,000 simulated data
sets, as a function of $K$, for each simulation setting discussed in Section
5.1. Each curve has been rescaled to take on values between $0$ and $1$, for
ease of comparison. The minimum loss values for each method are circled, and
$K^{*}$ is indicated by the vertical black line.

We further investigate the role of $\epsilon$ by repeating the simulation
study using different values of $\epsilon$ for single-fold data thinning. In
Figure 4, we plot the proportion of simulations that select the correct value
of $K^{*}$ in each of the three settings, as a function of $\epsilon$. We find
that in the Gamma clustering simulations, lower values of $\epsilon$ are
adequate. However, settings with weaker signal, such as the binomial PCA
example, require larger values of $\epsilon$ to identify the true latent
structure. In all settings, we observe that as $\epsilon$ approaches $1$,
performance begins to decay. This is a consequence of inadequate information
remaining in the test set under large values of $\epsilon$, and is consistent
with the discussion of Section 2.3. These findings suggest that in practice,
the optimal value of $\epsilon$ will be context-dependent.

Figure 4: The proportion of simulations for which data thinning selects the
true value of $K^{*}$ with the negative log-likelihood loss, as a function of
$\epsilon$, for the simulation study described in Section 5.1. The optimal
value of $\epsilon$ depends on the problem at hand.

Finally, we examine the benefits of multi-fold data thinning over single-fold
data thinning. Figure 5 displays histograms of the number of simulations that
select each value of $K$. Here we only include data thinning with
$\epsilon=0.8$ and multi-fold thinning with $M=5$, so that both methods use
the same allocation of information between training and test sets. We see that
multi-fold thinning generally selects the correct value of $K$ more often than
single-fold data thinning, mirroring the improvement of $M$-fold cross-
validation using sample splitting over single-fold sample splitting in
supervised settings. However, in the large Gamma setting, the signal is strong
enough that multi-fold thinning does not provide a benefit over single-fold
thinning.

Figure 5: The proportion of simulated data sets in which each candidate value
of $K$ is selected, with the negative log-likelihood loss, under data thinning
with $\epsilon=0.8$ and multifold thinning with $M=5$, for each of the
simulation settings described in Section 5.1. The true value of $K^{*}$ is
indicated by the vertical black line. Multifold thinning tends to select the
true value of $K$ more often than single-fold thinning.

Similar qualitative results are observed under mean squared error loss
(Appendix D.2).

##  6 Selecting the number of principal components in gene expression data

In this section, we revisit an analysis of a dataset from a single-cell RNA
sequencing (scRNA-seq) experiment conducted on a set of peripheral blood
mononuclear cells (PBMC). The dataset is freely available from 10X Genomics,
and was previously analyzed in the “Guided Clustering Tutorial” vignette
[Hoffman et al., 2022] for the popular R package Seurat, which is used for the
analysis of scRNA-seq datasets [Hao et al., 2021, Stuart et al., 2019, Satija
et al., 2015].

The PBMC dataset $X$ is a $2,700\times 32,738$ sparse matrix of non-negative
integers, representing counts from $32,738$ genes in each of $2,700$ cells. We
consider applying principal components analysis to learn a low-dimensional
representation of the data. In the Seurat vignette, filtering, normalization,
log-transformation, feature selection, centering, and scaling are applied to
the data, yielding a transformed matrix $\tilde{Y}\in\mathbb{R}^{2638\times
2000}$. Details are provided in Appendix E. Finally, they compute the singular
value decomposition (SVD) of $\tilde{Y}$, which we denote $\tilde{Y}=UDV^{T}$.
Here we let $U_{k}$ represent the $k$th column of the matrix $U$, and let
$U_{1:K}D_{1:K}V_{1:K}^{T}$ represent the rank-$K$ approximation of
$\tilde{Y}$.

Our goal is to select the number of principal components to use in this low-
rank approximation. In the Seurat vignette, the authors rely on heuristic
solutions such as looking for an “elbow” in the standard deviation of
$U_{K}D_{K}$ as a function of $K$; see Figure 6(a) [James et al., 2013]. Based
on the elbow plot, the authors suggest retaining around $7$ principal
components. Other heuristic approaches suggest as many as $12$ principal
components.

Before introducing the data thinning solution, we introduce a squared-error
based formulation that is mathematically equivalent to the traditional elbow
plot (see Appendix E), but will facilitate a direct comparison with data
thinning. For $K=1,\ldots,20$, we compute the sum of squared errors between
the matrix $\tilde{Y}$ and its rank-$K$ approximation:

$SSE_{K}(\tilde{Y})=\left\|\tilde{Y}-U_{1:K}D_{1:K}V_{1:K}^{T}\right\|_{F}^{2}.$

Because the low-rank approximation $U_{1:K}D_{1:K}V_{1:K}^{T}$ is computed
using $\tilde{Y}$, this loss function monotonically decreases with $K$.
Nonetheless, a heuristic solution for deciding how many principal components
to retain involves looking for the point in which the slope of the curve in
Figure 6(b) begins to flatten. This appears to happen around 5–7 principal
components. While this is consistent with the finding from Figure 6(a), the
exact number of principal components to retain is still unclear. We now show
that data thinning provides a less heuristic approach for estimating the
number of principal components.

Single-cell RNA-sequencing data are often modeled as Poisson random variables
[Townes et al., 2019, Wang et al., 2018, Sarkar and Stephens, 2021]. Thus, we
assume that $X_{ij}\sim\mathrm{Poisson}(\Lambda_{ij})$. Starting with the raw
data matrix $X\in\mathbb{Z}_{\geq 0}^{2700\times 32738}$, we perform Poisson
data thinning with $\epsilon=0.5$ to obtain a training set $X^{(1)}$ and a
test set $X^{(2)}$. Importantly, if our Poisson assumption holds then
$X^{(1)}$ and $X^{(2)}$ are independent. Furthermore, as $\epsilon=0.5$, they
are identically distributed. We then carry out the data processing described
in Appendix E on $X^{(1)}$ to obtain $\tilde{Y}^{(1)}\in\mathbb{R}^{2638\times
2000}$. We obtain $\tilde{Y}^{(2)}\in\mathbb{R}^{2638\times 2000}$ by applying
the same data processing steps to $X^{(2)}$, but retaining only the features
that were selected on $X^{(1)}$, so that the rows and columns of
$\tilde{Y}^{(1)}$ and $\tilde{Y}^{(2)}$ correspond to the same genes and
cells. Details are in Appendix E.

We compute the SVD on the training set,
$\tilde{Y}^{(1)}=U^{(1)}D^{(1)}(V^{(1)})^{T}$. For a range of values of $K$,
we then compute

$SSE_{K}(\tilde{Y}^{(1)},\tilde{Y}^{(2)})=\left\|\tilde{Y}^{(2)}-U_{1:K}^{(1)}D_{1:K}^{(1)}(V_{1:K}^{(1)})^{T}\right\|_{F}^{2}.$
(3)

The results are shown in Figure 6(c). As we not computing and evaluating the
SVD using the same data, the plot of $K$ vs. $SSE_{K}$ is not monotonically
decreasing in $K$. Instead, it reaches a clear minimum at $K=7$, suggesting
that the rank-$7$ approximation provides the best fit to the observed data.
Thus, data thinning provides a simple and non-heuristic way to select the
number of principal components.

Figure 6:  Results for the data analysis in Section 6. _(a)_ An “elbow plot”
of the standard deviation of the principal components, which reproduces the
plot given in the Seurat guided clustering tutorial. _(b)_ Due to the
relationship between $SSE_{K}(\tilde{Y})$ and the standard deviation of the
principal components (see Appendix E), looking for an elbow in (a) is
equivalent to looking for an elbow in (b). _(c)_ The data-thinning version of
(b), which shows a clear minimum in the SSE at 7 principal components.

We conclude with a few remarks.

######  Remark 5 (Choice of $\epsilon$).

Because we used $\epsilon=0.5$, $X^{(1)}$ and $X^{(2)}$ are identically
distributed, as are $\tilde{Y}^{(1)}$ and $\tilde{Y}^{(2)}$. If we had chosen
$\epsilon\neq 0.5$, then while
$\operatorname{\text{E}}[X^{(2)}]=\frac{1-\epsilon}{\epsilon}\operatorname{\text{E}}[X^{(1)}]$
(see Example 2.3), the relationship between
$\operatorname{\text{E}}[\tilde{Y}^{(1)}]$ and
$\operatorname{\text{E}}[\tilde{Y}^{(2)}]$ would depend on the details of the
data processing described in Appendix 1, and the loss function in (3) would
need to be modified accordingly.

######  Remark 6 (Overdispersion).

While we used a Poisson model for scRNA-seq data, there is evidence that a
negative binomial model may be preferable in some settings [Choudhary and
Satija, 2022]. It is possible to modify the analysis in this section using
negative binomial data thinning, as in Neufeld et al. [2023] (in preparation).

##  7 Discussion

In this paper, we proposed _data thinning_ , a general technique for
decomposing a random variable into two or more independent components. We then
applied data thinning to develop a version of cross-validation that is
suitable for unsupervised learning, a problem for which cross-validation with
sample splitting is not an option. Data thinning also has numerous
applications to supervised learning. It might be preferable to sample
splitting in supervised settings, where sample splitting is an option, in
cases where the sample size is small and where we wish to avoid excluding
high-leverage points from the training set [Leiner et al., 2022]. It may also
have power advantages to approaches such as sample splitting or selective
inference [Fithian et al., 2014] in problems such as inference after variable
selection. Finally, it can be used to estimate a model’s test set error for a
range of distributions, along the lines of existing work for the normal
distribution [Oliveira et al., 2021] and the Poisson distribution [Oliveira et
al., 2022].

While we focused on distributions for which the conditional distribution
$G_{\lambda_{1},\lambda_{2},x}$ takes a simple form (see Table 2), the
framework can easily be used to study convolution-closed distributions for
which $G_{\lambda_{1},\lambda_{2},x}$ takes an unfamiliar form; we just need
to sample from it.

In Section 2.4, we considered the impact of using the incorrect value of a
nuisance parameter when performing data thinning. However, we did not consider
what happens when the nuisance parameter is estimated using the data itself.
In future work, we will consider the theoretical and empirical implications of
performing data thinning with an estimated nuisance parameter.

In this paper, we focused on convolution-closed distributions and thus used
additive decompositions where $X=X^{(1)}+X^{(2)}$. For distributions with
bounded support, such as the Beta distribution, non-additive decompositions
are needed. We leave such decompositions to future work.

An R package implementing data thinning and scripts to reproduce the results
in this paper are available at https://anna-neufeld.github.io/datathin/.

##  8 Acknowledgements

Anna Neufeld, Ameer Dharamshi, and Daniela Witten were supported by the Simons
Foundation (Simons Investigator Award in Mathematical Modeling of Living
Systems) and the National Institutes of Health (R01 EB026908, R01 GM123993,
and R01 DA047869). Lucy Gao was supported by the Natural Sciences and
Engineering Research Council of Canada (Discovery Grants).

## References

  * Chen et al. [2021] Fan Chen, Sebastien Roch, Karl Rohe, and Shuqi Yu.  Estimating graph dimension with cross-validated eigenvalues.  _arXiv preprint arXiv:2108.03336_ , 2021.
  * Choudhary and Satija [2022] Saket Choudhary and Rahul Satija.  Comparison and evaluation of statistical error models for scRNA-seq.  _Genome Biology_ , 23(1):1–20, 2022.
  * Fithian et al. [2014] William Fithian, Dennis Sun, and Jonathan Taylor.  Optimal inference after model selection.  _arXiv preprint arXiv:1410.2597_ , 2014.
  * Fu and Perry [2020] Wei Fu and Patrick O Perry.  Estimating the number of clusters using cross-validation.  _Journal of Computational and Graphical Statistics_ , 29(1):162–173, 2020.
  * Gerard [2020] David Gerard.  Data-based RNA-seq simulations by binomial thinning.  _BMC Bioinformatics_ , 21(1):1–14, 2020.
  * Hao et al. [2021] Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M. Mauck III, Shiwei Zheng, Andrew Butler, et al.  Integrated analysis of multimodal single-cell data.  _Cell_ , 2021.
  * Hastie et al. [2009] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman.  _The Elements of Statistical Learning: Data Mining, Inference, and Prediction_ , volume 2.  Springer, 2009.
  * Hoffman et al. [2022] P. Hoffman et al.  Seurat guided clustering tutorial.  _https://satijalab.org/seurat/ articles/pbmc3k_tutorial.html_ , 2022\.  Accessed: 09-12-2022.
  * James et al. [2013] Gareth James, Daniela Witten, Trevor Hastie, and Robert Tibshirani.  _An Introduction to Statistical Learning_ , volume 112.  Springer, 2013.
  * Joe [1996] Harry Joe.  Time series models with univariate margins in the convolution-closed infinitely divisible class.  _Journal of Applied Probability_ , 33(3):664–677, 1996.
  * Jørgensen [1992] Bent Jørgensen.  Exponential dispersion models and extensions: A review.  _International Statistical Review/Revue Internationale de Statistique_ , pages 5–20, 1992.
  * Jørgensen and Song [1998] Bent Jørgensen and Peter Xue-Kun Song.  Stationary time series models with exponential dispersion model margins.  _Journal of Applied Probability_ , 35(1):78–92, 1998.
  * Leiner et al. [2022] James Leiner, Boyan Duan, Larry Wasserman, and Aaditya Ramdas.  Data fission: splitting a single data point.  _arXiv preprint arXiv:2112.11079_ , 2022.
  * Louzada et al. [2019] Francisco Louzada, Pedro L. Ramos, and Eduardo Ramos.  A note on bias of closed-form estimators for the gamma distribution derived from likelihood equations.  _The American Statistician_ , 73(2):195–199, 2019\.
  * Neufeld et al. [2022] Anna Neufeld, Lucy L Gao, Joshua Popp, Alexis Battle, and Daniela Witten.  Inference after latent variable estimation for single-cell RNA sequencing data.  _Biostatistics_ , 2022.
  * Neufeld et al. [2023] Anna Neufeld, Lucy Gao, Joshua Popp, Alexis Battle, and Daniela Witten.  Negative binomial count splitting for single-cell RNA sequencing data.  _To be submitted_ , 2023.
  * Oliveira et al. [2021] Natalia L Oliveira, Jing Lei, and Ryan J Tibshirani.  Unbiased risk estimation in the normal means problem via coupled bootstrap techniques.  _arXiv preprint arXiv:2111.09447_ , 2021.
  * Oliveira et al. [2022] Natalia L Oliveira, Jing Lei, and Ryan J Tibshirani.  Coupled bootstrap test error estimation for Poisson variables.  _arXiv preprint arXiv:2212.01943_ , 2022.
  * Owen and Perry [2009] Art B Owen and Patrick O Perry.  Bi-cross-validation of the SVD and the nonnegative matrix factorization.  _The Annals of Applied Statistics_ , 3(2):564–594, 2009.
  * Rasines and Young [2021] Daniel G Rasines and G Alastair Young.  Splitting strategies for post-selection inference.  _arXiv preprint arXiv:2102.02159_ , 2021.
  * Sarkar and Stephens [2021] Abhishek Sarkar and Matthew Stephens.  Separating measurement and expression models clarifies confusion in single-cell rna sequencing analysis.  _Nature Genetics_ , 53(6):770–777, 2021.
  * Satija et al. [2015] Rahul Satija, Jeffrey A Farrell, David Gennert, Alexander F Schier, and Aviv Regev.  Spatial reconstruction of single-cell gene expression data.  _Nature Biotechnology_ , 33:495–502, 2015.
  * Stuart et al. [2019] Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister, Efthymia Papalexi, William M Mauck III, et al.  Comprehensive integration of single-cell data.  _Cell_ , 177:1888–1902, 2019.
  * Tian [2020] Xiaoying Tian.  Prediction error after model search.  _The Annals of Statistics_ , 48(2):763–784, 2020\.
  * Townes et al. [2019] F. William Townes, Stephanie C Hicks, Martin J Aryee, and Rafael A Irizarry.  Feature selection and dimension reduction for single-cell rna-seq based on a multinomial model.  _Genome Biology_ , 20(1):1–16, 2019.
  * Wang et al. [2018] Jingshu Wang, Mo Huang, Eduardo Torre, Hannah Dueck, Sydney Shaffer, John Murray, et al.  Gene expression distribution deconvolution in single-cell RNA sequencing.  _Proceedings of the National Academy of Sciences_ , 115(28):E6437–E6446, 2018.
  * Ye and Chen [2017] Zhi-Sheng Ye and Nan Chen.  Closed-form estimators for the gamma distribution derived from likelihood equations.  _The American Statistician_ , 71(2):177–181, 2017\.

##  Appendix A Proofs from Section 2

###  A.1 Proof of Theorem 1

This proof is due to Joe [1996] and Jørgensen and Song [1998], but has been
adapted to fit our notation.

Let $x$, our observed data, be a realization of a random variable $X\sim
F_{\lambda}$. Let $\epsilon\in(0,1)$ be chosen such that $\epsilon\lambda$ and
$(1-\epsilon)\lambda$ are in the parameter space $\Lambda$. We draw
$X^{(1)}\mid X=x\sim G_{\epsilon\lambda,(1-\epsilon)\lambda,x}$, where this
notation was defined in Section 2.1, and let $X^{(2)}=X-X^{(1)}$.

Separately, let $X^{\prime}\sim F_{\epsilon\lambda}$ and $X^{\prime\prime}\sim
F_{(1-\epsilon)\lambda}$ be independent, and let
$Y=X^{\prime}+X^{\prime\prime}$. As $F_{\lambda}$ is a convolution-closed
distribution, it follows that $Y\sim F_{\lambda}$ and thus we know that $Y$
has the same marginal distribution as $X$.

We first argue that the joint distribution of $(X^{(1)},X)$ is the same as the
joint distribution of $(X^{\prime},Y)$. The conditional distribution of
$X^{(1)}\mid X$ is the same as the conditional distribution of $X^{\prime}\mid
Y$ by definition of the distribution
$G_{\epsilon\lambda,(1-\epsilon)\lambda,x}$. Furthermore, we have already
noted that the marginal distributions of $X$ and $Y$ are the same. Thus, the
joint distribution of $(X^{(1)},X)$ is the same as the joint distribution of
$(X^{\prime},Y)$.

As $X^{(2)}$ is deterministic given $X^{(1)}$ and $X$, the joint distribution
of $(X^{(1)},X^{(2)})$ is the same as the joint distribution of $(X^{(1)},X)$.
Similarly, the joint distribution of $(X^{\prime},Y)$ is the same as the joint
distribution of $(X^{\prime},X^{\prime\prime})$. Thus, the joint distribution
of $(X^{(1)},X^{(2)})$ is the same as the joint distribution of
$(X^{\prime},X^{\prime\prime})$. As the joint distribution of $X^{\prime}$ and
$X^{\prime\prime}$ is known to be the product of independent distributions
$F_{\epsilon\lambda}$ and $F_{(1-\epsilon)\lambda}$, this completes the proof
of parts 1-3 of Theorem 1.

The final statement of Theorem 1 follows directly from Remark 2.

###  A.2 Proof of Proposition 1

To prove (i), note that since $X^{(1)}\mid X=x$ is normally distributed and
$X$ is normally distributed, a well-known property of the normal distribution
tells us that the marginal distribution of $X^{(1)}$ is normal. We then use
the law of total expectation and the law of total variance to compute its mean
and variance.

$\displaystyle\operatorname{\text{E}}[X^{(1)}]$
$\displaystyle=\operatorname{\text{E}}[\operatorname{\text{E}}[X^{(1)}\mid
X]]=\operatorname{\text{E}}[\epsilon X]=\epsilon\mu$
$\displaystyle\operatorname{\text{Var}}(X^{(1)})$
$\displaystyle=\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[X^{(1)}\mid
X\right]\right)+\operatorname{\text{E}}\left[\operatorname{\text{Var}}\left(X^{(1)}\mid
X\right)\right]$ $\displaystyle=\operatorname{\text{Var}}\left(\epsilon
X\right)+\operatorname{\text{E}}\left(\epsilon(1-\epsilon)\tilde{\sigma}^{2}\right)$
$\displaystyle=\epsilon^{2}\sigma^{2}+\epsilon(1-\epsilon)\tilde{\sigma}^{2}.$

To prove (ii), note that the difference between two normally distributed
variables ($X$ and $X^{(1)}$) is normal. Then note that

$\displaystyle\operatorname{\text{E}}[X^{(2)}]$
$\displaystyle=\operatorname{\text{E}}[X]-\operatorname{\text{E}}[X^{(1)}]=\mu-\epsilon\mu=(1-\epsilon)\mu.$
$\displaystyle\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[X^{(2)}\mid
X\right]\right)+\operatorname{\text{E}}\left[\operatorname{\text{Var}}\left(X^{(2)}\mid
X\right)\right]$
$\displaystyle=\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[X-X^{(1)}\mid
X\right]\right)+\operatorname{\text{E}}\left[\operatorname{\text{Var}}\left(X-X^{(1)}\mid
X\right)\right]$
$\displaystyle=\operatorname{\text{Var}}\left((1-\epsilon)X\right)+\operatorname{\text{E}}\left[\operatorname{\text{Var}}\left(X^{(1)}\mid
X\right)\right]$
$\displaystyle=(1-\epsilon)^{2}\sigma^{2}+\epsilon(1-\epsilon)\tilde{\sigma}^{2},$

which completes the proof of (ii). Finally, to prove (iii), note that

$\displaystyle 2\operatorname{\text{Cov}}(X^{(1)},X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}(X)-\operatorname{\text{Var}}(X^{(1)})-\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=\sigma^{2}-\epsilon^{2}\sigma^{2}-\epsilon(1-\epsilon)\tilde{\sigma}^{2}-(1-\epsilon)^{2}\sigma^{2}-\epsilon(1-\epsilon)\tilde{\sigma}^{2}$
$\displaystyle=2\epsilon(1-\epsilon)\left(\sigma^{2}-\tilde{\sigma}^{2}\right).$

###  A.3 Proof of Proposition 2

Recall that if $A\sim\mathrm{BetaBinomial}\left(r,\alpha,\beta\right)$, then
$\operatorname{\text{E}}[A]=\frac{r\alpha}{\alpha+\beta}$ and
$\operatorname{\text{Var}}(A)=\frac{r\alpha\beta(\alpha+\beta+r)}{(\alpha+\beta)^{2}(\alpha+\beta+1)}$.
Then we can derive the marginal variance of $X^{(1)}$ using the law of total
variance and the fact that $X^{(1)}\mid
X\sim\mathrm{BetaBinomial}(X,\epsilon\tilde{r},(1-\epsilon)\tilde{r})$.

$\displaystyle\operatorname{\text{Var}}(X^{(1)})$
$\displaystyle=\operatorname{\text{E}}[\operatorname{\text{Var}}(X^{(1)}\mid
X)]+\operatorname{\text{Var}}(\operatorname{\text{E}}[X^{(1)}\mid X])$
$\displaystyle=\operatorname{\text{E}}\left[\frac{X\epsilon(1-\epsilon)(\tilde{r}+X)}{(\tilde{r}+1)}\right]+\operatorname{\text{Var}}\left(\epsilon
X\right)$
$\displaystyle=\frac{\epsilon(1-\epsilon)}{\tilde{r}+1}\left(\tilde{r}\operatorname{\text{E}}\left[X\right]+\operatorname{\text{E}}\left[X^{2}\right]\right)+\epsilon^{2}\operatorname{\text{Var}}\left(X\right)$
$\displaystyle=\frac{\epsilon(1-\epsilon)}{\tilde{r}+1}\left(\tilde{r}\operatorname{\text{E}}\left[X\right]+\operatorname{\text{Var}}(X)+E[X]^{2}\right)+\epsilon^{2}\operatorname{\text{Var}}\left(X\right).$

Next note that $\operatorname{\text{Var}}(X^{(2)}\mid
X)=\operatorname{\text{Var}}(X-X^{(1)}\mid
X)=\operatorname{\text{Var}}(X^{(1)}\mid X)$ and that
$\operatorname{\text{E}}[X^{(2)}\mid X]=(1-\epsilon)X$. Thus, we arrive at:

$\displaystyle\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=\frac{\epsilon(1-\epsilon)}{\tilde{r}+1}\left(\tilde{r}\operatorname{\text{E}}\left[X\right]+\operatorname{\text{Var}}(X)+E[X]^{2}\right)+(1-\epsilon)^{2}\operatorname{\text{Var}}\left(X\right).$

To derive the covariance, note that:

$\displaystyle 2\operatorname{\text{Cov}}(X^{(1)},X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}(X)-\operatorname{\text{Var}}(X^{(1)})-\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}(X)-2\frac{\epsilon(1-\epsilon)}{\tilde{r}+1}\left(\tilde{r}\operatorname{\text{E}}\left[X\right]+\operatorname{\text{Var}}(X)+E[X]^{2}\right)-\left(\epsilon^{2}+(1-\epsilon)^{2}\right)\operatorname{\text{Var}}\left(X\right)$
$\displaystyle=-2\frac{\epsilon(1-\epsilon)}{\tilde{r}+1}\left(\tilde{r}r\frac{1-p}{p}+r\frac{1-p}{p^{2}}+r^{2}\frac{(1-p)^{2}}{p^{2}}\right)+2\epsilon(1-\epsilon)r\frac{1-p}{p^{2}}$
$\displaystyle=2\epsilon(1-\epsilon)r\frac{(1-p)^{2}}{p^{2}}\left(1-\frac{r+1}{\tilde{r}+1}\right).$

###  A.4 Proof of Proposition 3

The proof structure is identical to those of Proposition 1 and Proposition 2.
We start by recalling that if $A\sim\mathrm{Gamma}(\alpha,\beta)$, then
$\operatorname{\text{E}}[X]=\frac{\alpha}{\beta}$ and
$\operatorname{\text{Var}}(X)=\frac{\alpha}{\beta^{2}}$. Then:

$\displaystyle\operatorname{\text{Var}}(X^{(1)})$
$\displaystyle=E\left[\operatorname{\text{Var}}\left(X^{(1)}\mid
X\right)\right]+\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[X^{(1)}\mid
X\right]\right)$ $\displaystyle=E\left[\operatorname{\text{Var}}\left(XZ\mid
X\right)\right]+\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[XZ\mid
X\right]\right)$
$\displaystyle=E\left[X^{2}\operatorname{\text{Var}}\left(Z\right)\right]+\operatorname{\text{Var}}\left(X\operatorname{\text{E}}\left[Z\right]\right)$
$\displaystyle=E\left[X^{2}\left(\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\right)\right]+\operatorname{\text{Var}}\left(X\epsilon\right)$
$\displaystyle=\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\left(\operatorname{\text{Var}}(X)+\operatorname{\text{E}}[X]^{2}\right)+\epsilon^{2}\operatorname{\text{Var}}\left(X\right).$

Similarly, we note that $\operatorname{\text{Var}}(X^{(2)}\mid
X)=\operatorname{\text{Var}}(X-X^{(1)}\mid
X)=\operatorname{\text{Var}}(X^{(1)}\mid X)$, while
$\operatorname{\text{E}}(X^{(2)}\mid X)=(1-\epsilon)X$. This allows us to do a
similar derivation and arrive at:

$\displaystyle\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=E\left[\operatorname{\text{Var}}\left(X^{(2)}\mid
X\right)\right]+\operatorname{\text{Var}}\left(\operatorname{\text{E}}\left[X^{(2)}\mid
X\right]\right)$
$\displaystyle=\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\left(\operatorname{\text{Var}}(X)+\operatorname{\text{E}}[X]^{2}\right)+(1-\epsilon)^{2}\operatorname{\text{Var}}\left(X\right).$

Finally,

$\displaystyle 2\operatorname{\text{Cov}}(X^{(1)},X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}(X)-\operatorname{\text{Var}}(X^{(1)})-\operatorname{\text{Var}}(X^{(2)})$
$\displaystyle=\operatorname{\text{Var}}(X)-2\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\left(\operatorname{\text{Var}}(X)+\operatorname{\text{E}}[X]^{2}\right)-\left(\epsilon^{2}+1-2\epsilon+\epsilon^{2}\right)\operatorname{\text{Var}}(X)$
$\displaystyle=-2\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\left(\operatorname{\text{Var}}(X)+\operatorname{\text{E}}[X]^{2}\right)+2\epsilon(1-\epsilon)\operatorname{\text{Var}}(X)$
$\displaystyle=-2\frac{\epsilon(1-\epsilon)}{\tilde{\alpha}+1}\left(\frac{\alpha(1+\alpha)}{\beta^{2}}\right)+2\epsilon(1-\epsilon)\frac{\alpha}{\beta^{2}}$
$\displaystyle=2\epsilon(1-\epsilon)\frac{\alpha}{\beta^{2}}\left(1-\frac{\alpha+1}{\tilde{\alpha}+1}\right).$

##  Appendix B Proof of Theorem 2

The proof is nearly identical to that of Theorem 1. It extends ideas from
Jørgensen and Song [1998] and Joe [1996] to the setting of multiple folds.

Let $x$, our observed data, be a realization of random variable $X\sim
F_{\lambda}$. Let $\epsilon_{1},\ldots,\epsilon_{M}$ be chosen such that
$\sum_{m=1}^{M}\epsilon_{m}=1$, $\epsilon_{m}>0$, and $\epsilon_{m}\lambda$ is
in the parameter space $\Lambda$ for $m=1,\ldots,M$.

Suppose we draw $\left(X^{(1)},\ldots,X^{(M)}\right)\mid X=x\sim
G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\ldots,\epsilon_{M}\lambda,x}$,
where
$G_{\epsilon_{1}\lambda,\epsilon_{2}\lambda,\ldots,\epsilon_{M}\lambda,x}$ was
defined in Section 3.

Separately, let $X_{1},X_{2},\ldots,X_{M}$ be independent random variables,
where $X_{m}\sim F_{\epsilon_{m}\lambda}$, and let $Y=\sum_{m=1}^{M}X_{m}$. As
$F_{\lambda}$ is a convolution-closed distribution, we know that $Y\sim
F_{\lambda}$ and thus $Y$ has the same marginal distribution as $X$.

The conditional distribution of $(X^{(1)},X^{(2)},\ldots,X^{(M)})\mid X$ is
the same as the conditional distribution of $(X_{1},X_{2},\ldots,X_{M})\mid
Y$, by definition of the distribution
$G_{\epsilon_{1}\lambda,\ldots,\epsilon_{M}\lambda,x}$. Furthermore, we have
already seen that the marginal distributions of $X$ and $Y$ are the same.
Thus, the marginal distribution of $(X^{(1)},X^{(2)},\ldots,X^{(M)})$ is the
same as the marginal distribution of $(X_{1},X_{2},\ldots,X_{M})$.
Furthermore, by construction, the marginal distribution of
$(X_{1},X_{2},\ldots,X_{M})$ is that of $M$ independent random variables,
where $X_{m}\sim F_{\epsilon_{m}\lambda}$. This concludes the proof of Theorem
2, parts 1-3. Part 4 of Theorem 2 follows directly from Remark 2.

##  Appendix C Comparing two decompositions of the normal distribution

Consider the task of splitting the $\mathrm{N}(\mu,\sigma^{2})$ distribution
into two independent normally-distributed random variables, with $\sigma$
known. The data thinning proposal is given in Table 2, and leads to
$X^{(1)}\sim\mathrm{N}(\epsilon\mu,\epsilon\sigma^{2})$ and
$X^{(2)}\sim\mathrm{N}((1-\epsilon)\mu,(1-\epsilon)\sigma^{2})$.

The data fission proposal is as follows: given a value of $\tau>0$, we draw
$Z\sim\mathrm{N}(0,\sigma^{2})$, and then let $X^{(1)}=X+\tau Z$ and
$X^{(2)}=X-\frac{1}{\tau}Z$. Then,
$X^{(1)}\sim\mathrm{N}(\mu,(1+\tau^{2})\sigma^{2})$ and
$X^{(2)}\sim\mathrm{N}(\mu,(1+\frac{1}{\tau^{2}})\sigma^{2})$. It follows that
$X^{(1)}\perp\\!\\!\\!\perp X^{(2)}$.

Under the data thinning decomposition, the Fisher Information contained in
$X^{(1)}$ for the parameter $\mu$ is
$I_{X^{(1)}}(\mu)=\frac{\epsilon}{\sigma^{2}}=\epsilon I_{X}(\mu)$. Similarly,
$I_{X^{(2)}}=(1-\epsilon)I_{X}(\mu)$. Under the data fission decomposition,
$I_{X^{(1)}}(\mu)=\frac{1}{(1+\tau^{2})\sigma^{2}}=\frac{1}{1+\tau^{2}}I_{X}(\mu)$,
and $I_{X^{(2)}}(\mu)=\frac{\tau^{2}}{1+\tau^{2}}I_{X}(\mu)$. Therefore, if
$\epsilon=\frac{1}{1+\tau^{2}}$, then the two decompositions result in the
same amount of information about the parameter of interest, $\mu$, in
$X^{(1)}$ and $X^{(2)}$. Hence, we expect the data thinning and data fission
decompositions of a $\mathrm{N}(\mu,\sigma^{2})$ random variable to behave
similarly in applications if $\epsilon=\frac{1}{1+\tau^{2}}$.

##  Appendix D Simulation Study Supporting Details

Here we provide additional details for the simulation study described in
Section 5.

###  D.1 Simulation study parameters

The following three subsections provide additional details about the
simulation studies described in Section 5.1.

####  D.1.1 Section 5.1.1: Selecting the number of principal components for
binomial data

For $K^{*}=10$, we compute $\theta=UDV^{T}$ where $U$ is a $n\times K^{*}$
random orthogonal matrix, $D$ is a $K^{*}\times K^{*}$ diagonal matrix with
diagonal elements equal to $5,6,\dots,14$, and $V$ is a $d\times K^{*}$ random
orthogonal matrix. Then,
$p_{ij}=\frac{\exp{\left(\theta_{ij}\right)}}{1+\exp{\left(\theta_{ij}\right)}}$
for $i=1,\dots,n$ and $j=1,\dots,d$.

####  D.1.2 Section 5.1.2: Selecting the number of clusters in Gamma-
distributed data where $d=2$ and $K^{*}=4$

In the small $d$, small $K^{*}$ clustering setting described in Section 5.1.2,
observations from each cluster are generated as
$X_{ij}\overset{\textup{ind}}{\sim}\text{Gamma}(\lambda,\theta_{c_{i},j})$
where $\lambda=20$,

$\theta=\begin{bmatrix}0.5&5\\\ 5&0.5\\\ 10&10\\\ 0.5&0.5\end{bmatrix},$

and $c_{i}\in\\{1,2,3,4\\}$ is the true cluster membership for the $i$th
observation.

####  D.1.3 Section 5.1.2: Selecting the number of clusters in Gamma-
distributed data where $d=100$ and $K^{*}=10$

In the large $d$, large $K^{*}$ clustering setting described in Section 5.1.2,
observations from each cluster are generated as
$X_{ij}\overset{\textup{ind}}{\sim}\text{Gamma}(\lambda,\theta_{c_{i},j})$
where $\lambda=2$, the $K^{*}\times d$ matrix $\theta$ is constructed such
that for $j=1,\dots,d$ and $k=1,\dots,K^{*}$,

$\theta_{kj}=\begin{cases}0.1&\text{if }k\leq 9\text{ and }10k-9\leq j\leq
10k+10,\\\ 1&\text{otherwise},\end{cases}$

and $c_{i}\in\\{1,2,\dots,10\\}$ is the true cluster membership for the $i$th
observation.

###  D.2 Mean squared error loss function

####  D.2.1 Methods

As an alternative to the negative log-likelihood loss used in Section 5, here
we consider applying a mean squared error (MSE) loss. Algorithms 5 and 6 are
modifications of Algorithms 3 and 4 that make use of MSE instead of NLL.

######  Algorithm 5 (Evaluating binomial PCA with mean squared error loss).

.
Input: A positive integer $K$, and matrices
$X^{(\mathrm{train})}\in\mathbb{R}^{n\times d}$ and
$X^{(\mathrm{test})}\in\mathbb{R}^{n\times d}$, where
$X^{(\mathrm{train})}_{ij}\overset{\textup{ind}}{\sim}\textup{Binomial}\left(a^{\textup{(train)}}r,p_{ij}\right)$,
$X^{(\mathrm{test})}_{ij}\overset{\textup{ind}}{\sim}\textup{Binomial}\left(a^{\textup{(test)}}r,p_{ij}\right)$,
$r$ is the known number of trials, and $a^{\textup{(train)}}$ and
$a^{\textup{(test)}}$ are known constants in $(0,1]$ such that
$a^{\textup{(train)}}r$ and $a^{\textup{(test)}}r$ are integers.

  1. 1.

Compute the singular value decomposition of the observed log-odds matrix,
i.e.
$\textup{logit}\left(\frac{X^{(\mathrm{train})}_{ij}+0.001}{a^{\textup{(train)}}r+0.002}\right)=\hat{U}\hat{D}\hat{V}^{T}$.
Note that the pseudo-counts prevent taking the logit of 0 or 1.

  2. 2.

Construct the rank-$K$ approximation,
$p^{(K)}=\textup{expit}\left(\hat{U}_{1:K}\hat{D}_{1:K}\hat{V}_{1:K}\right)$.

  3. 3.

Compute the loss function:

$\textup{MSE}=\frac{1}{nd}\sum_{i=1}^{n}\sum_{j=1}^{d}\left(X^{(\mathrm{test})}_{ij}-a^{\textup{(test)}}rp_{ij}^{(K)}\right)^{2}.$
(4)

######  Algorithm 6 (Evaluating Gamma clusters with mean squared error loss).

.
Input: A positive integer $K$, and matrices
$X^{(\mathrm{train})}\in\mathbb{R}^{n\times d}$ and
$X^{(\mathrm{test})}\in\mathbb{R}^{n\times d}$, where
$X^{(\mathrm{train})}_{ij}\overset{\textup{ind}}{\sim}\textup{Gamma}\left(a^{\textup{(train)}}\lambda,\theta_{c_{i},j}\right)$,
$X^{(\mathrm{test})}_{ij}\overset{\textup{ind}}{\sim}\textup{Gamma}\left(a^{\textup{(test)}}\lambda,\theta_{c_{i},j}\right)$,
$\theta\in(0,\infty)^{K^{*}\times d}$ where $\theta_{c_{i},j}$ is the true but
unknown rate parameter for the $c_{i}$th cluster in the $j$th dimension where
$c_{i}\in\\{1,2,\dots,K^{*}\\}$, $\lambda$ is the known shape parameter, and
$a^{\textup{(train)}}$ and $a^{\textup{(test)}}$ are known constants in
$(0,1]$.

  1. 1.

Run $K$-means on $X^{(\mathrm{train})}$ to estimate $K$ clusters. Denote the
cluster assignment of the $i$th observation as $\hat{c}_{i}$.

  2. 2.

Compute the matrix of cluster means,
$\hat{\mu}^{\textup{(K)}}\in\mathbb{R}^{K\times d}$.

  3. 3.

Compute the loss function:

$\textup{MSE}=\frac{1}{nd}\sum_{i=1}^{n}\sum_{j=1}^{d}\left(X^{(\mathrm{test})}_{ij}-\frac{a^{\textup{(test)}}}{a^{\textup{(train)}}}\hat{\mu}_{\hat{c}_{i},j}^{\textup{(K)}}\right)^{2}.$
(5)

Algorithms 5 and 6 are then applied in the naive, data thinning, and multi-
fold thinning cases as described in Section 5.2.3.

####  D.2.2 Results

In Figure 7, we plot the average mean squared error curves, as a function of
$K$. As with the negative log-likelihood loss, data thinning approaches
produce curves with sharp minimum values at or near $K=K^{*}$, as opposed to
the naive method’s monotonically-decreasing curves.

Figure 7: The mean squared error loss averaged over 2,000 simulated data sets,
as a function of $K$, for each simulation setting described in Section 5.1.
Each curve has been rescaled to take on values between $0$ and $1$, for ease
of comparison. The minimum loss values for each method are circled, and
$K^{*}$ is indicated by the vertical black line.

In Figure 8, we plot the proportion of simulations that select the correct
value of $K^{*}$ using the mean squared error loss, as a function of
$\epsilon$. Results are largely similar to Figure 4.

Figure 8: The proportion of simulations for which data thinning selects the
true value of $K^{*}$ with the mean squared error loss, as a function of
$\epsilon$, for the simulation study described in Section 5.1. The optimal
value of $\epsilon$ depends on the problem at hand.

Finally, we compare multi-fold to single-fold thinning, under the mean squared
error loss, in Figure 9. As in Figure 5, we find that multi-fold thinning
tends to select the correct value of $K$ more often than single-fold thinning.

Figure 9: The proportion of simulated data sets in which each candidate value
of $K$ is selected, with the mean squared error loss, under data thinning with
$\epsilon=0.8$ and multifold thinning with $M=5$, for each of the simulation
settings described in Section 5.1. The true value of $K^{*}$ is indicated by
the vertical black line. Multifold thinning tends to select the true value of
$K$ more often than single-fold thinning.

##  Appendix E Details for the real data analysis in Section 6

We first explain in detail the preprocessing done to the matrix $X$ in the
Seurat tutorial.

  * (1)

Initial data filtering: We initially filter the data such that only cells with
between 200 and 2500 total Unique Molecular Identifiers (UMIs) remain (with
fewer than 5% of the UMIs coming from from mitochondrial genes) and only genes
that are expressed in at least 200 cells remain. This reduces the size of $X$
from $2,700\times 32,738$ to $2,638\times 13,714$

  * (2)

Log normalization:  Next, the data are normalized and log transformed, such
that

$Y_{ij}=\log\left(\frac{X_{ij}}{\sum_{t=1}^{13714}X_{it}}\times
10,000+1\right).$

  * (3)

Feature selection: Following this transformation, the top 2000 highly variable
genes are selected using the function FindVariableFeatures from the Seurat
package. The goal of the function is to find a subset of features with high
cell-to-cell variation after accounting for the inherent mean-variance
relationship, as these are most likely to be interesting in downstream
analysis, and it implements methodology from Stuart et al. [2019].

  * (4)

Centering and scaling: Finally, the columns of the subsetted matrix
$Y\in\mathbb{R}^{2638\times 2000}$ are centered and scaled to obtain the
matrix $\tilde{Y}$.

After these preprocessing steps, the principal components of $\tilde{Y}$ are
computed.

We now explain the preprocessing for $X^{(1)}$ and $X^{(2)}$ that we use for
our data thinning alternative to the Seurat tutorial. We follow the same four
steps as above, but we are careful to specify what we do on the training set
$X^{(1)}$ as opposed to the test set $X^{(2)}$.

  * (1)

Initial data filtering: We perform the initial data filtering from Step (1)
above on $X^{(1)}$. We then subset $X^{(2)}$ to include the same genes and
cells as those in $X^{(1)}$. After this step, $X^{(1)}$ and $X^{(2)}$ are both
in $\mathbb{Z}_{\geq 0}^{2638\times 13258}$.

  * (2)

Log normalization:  We normalize and log-transform both $X^{(1)}$ and
$X^{(2)}$, such that:

$Y_{ij}^{(1)}=\log\left(\frac{X^{(1)}_{ij}}{\sum_{t=1}^{13258}X^{(1)}_{it}}\times
10,000+1\right),Y_{ij}^{(2)}=\log\left(\frac{X^{(2)}_{ij}}{\sum_{t=1}^{13258}X^{(2)}_{it}}\times
10,000+1\right).$

We note that these random variables are still independent and identically
distributed under our Poisson assumption.

  * (3)

Feature selection: We then apply the Seurat function FindVariableFeatures to
the matrix $Y^{(1)}$ to select the top $2000$ highly variable genes [Stuart et
al., 2019]. We subset both $Y^{(1)}$ and $Y^{(2)}$ to contain only these
genes, such that $Y^{(1)},Y^{(2)}\in\mathbb{R}^{2638\times 2000}$.

  * (4)

Centering and scaling: We center and scale the columns of the subsetted
$Y^{(1)}$ to obtain $\tilde{Y}^{(1)}$. We also center and scale the columns of
the subsetted $Y^{(2)}$ to obtain $\tilde{Y}^{(2)}$.

After these preprocessing steps, the principal components of $\tilde{Y}^{(1)}$
are computed, and the loss function is computed using $\tilde{Y}^{(2)}$.

We now explain the identity that makes Figure 6(a) and Figure 6(b)
mathematically equivalent. In Section 6, we defined

$SSE_{K}(\tilde{Y})=\left\|\tilde{Y}-U_{1:K}D_{1:K}V_{1:K}^{T}\right\|_{F}^{2}.$

We see that:

$\displaystyle\left\|\tilde{Y}-U_{1:K}D_{1:K}V_{1:K}^{T}\right\|_{F}^{2}$
$\displaystyle=\|\tilde{Y}\|_{2}^{2}-2\mathrm{trace}\left(\tilde{Y}^{T}U_{1:K}D_{1:K}V_{1:K}^{T}\right)+\mathrm{trace}\left(V_{1:K}D^{T}_{1:K}U_{1:K}^{T}U_{1:K}D_{1:K}V_{1:K}^{T}\right)$
$\displaystyle=\|\tilde{Y}\|_{F}^{2}-\mathrm{trace}\left(D_{1:K}^{T}D_{1:K}\right)=\left\|\tilde{Y}\right\|_{F}^{2}-\sum_{j=1}^{K}D_{jj}^{2}.$

Thus, if we compute $SSE_{K}(\tilde{Y})$ for $K=1,\ldots,20$, since
$\left\|\tilde{Y}\right\|_{F}^{2}$ is fixed, we can easily obtain the values
of $\sum_{j=1}^{K}D_{jj}^{2}$ for $K=1,\ldots,20$. By taking the differences
between these values for $K$ and $K+1$, we obtain $D_{KK}$ for
$K=1,\ldots,20$. Finally, we note that the standard deviation of the $K$th
principal component $(U_{K}D_{KK})$ can be written as

$\sqrt{(D_{KK}U_{K})^{T}(D_{KK}U_{K})}=\sqrt{D_{KK}^{2}}.$

Since the standard deviation of the $K$th principal component (plotted in
Figure 6(a)) can be obtained directly from the sums of squared error plotted
in Figure 6(b), we say that the two plots are mathematically equivalent.