capbtabboxtable[][]

11footnotetext:  Department of Statistics, Florida State University,
Tallahassee, FL 32306.22footnotetext:  Department of Statistics, University of
Chicago, Chicago, IL, 60637.

# Testing and estimation for clustered signals

Hongyuan Cao1 and Wei Biao Wu2

###### Abstract

We propose a change-point detection method for large scale multiple testing
problems with data having clustered signals. Unlike the classic change-point
setup, the signals can vary in size within a cluster. The clustering structure
on the signals enables us to effectively delineate the boundaries between
signal and non-signal segments. New test statistics are proposed for
observations from one and/or multiple realizations. Their asymptotic
distributions are derived. We also study the associated variance estimation
problem. We allow the variances to be heteroscedastic in the multiple
realization case, which substantially expands the applicability of the
proposed method. Simulation studies demonstrate that the proposed approach has
a favorable performance. Our procedure is applied to an array based
Comparative Genomic Hybridization (aCGH) dataset.

Keywords:  Change-point inference, clustered signal, high dimension, multiple
testing, signal aggregation, variance estimation

##  1 Introduction

Signal detection and multiple testing in a data rich environment have been
important research topics in natural and social sciences. Typical examples
include detecting anomalous traffic in computer networks [33], identifying
voxels that correlate with certain activities [15] in functional Magnetic
Resonance Imaging (fMRI), and associating single nucleotide polymorphisms
(SNPs) with clinical outcomes [29]. The predominant framework in these
research is via individual analysis–testing each hypothesis separately and
declaring statistical significance if the $p$-value is less than certain
threshold [2] or the two-sample $t$-statistic falls into the rejection region
[12, 5]. Various approaches were proposed to improve the power by
incorporating structured or prior information. For example [4] studied group
hypothesis testing; [13, 16] investigated $p$-value weighting; [6] considered
$p$-value aggregating; [10, 19] utilized prior experimental information on
each hypothesis in the inference stage with data from a new experiment; and
[24, 23] harnessed the sparsity of mean vectors with student’s $t$-statistics.

For data with signals having clustered structure, multiple testing approaches
currently in use fall into two general classes. The first approach defines
possible regions of interest in advance, either by field knowledge or an
independent experiment. [32] proposed a spatial testing procedure with pre-
specified regions of interest in a compound theoretical framework; [15]
developed an algorithm specifically tailored for brain imaging data where a
preliminary scan is used to select clusters by grouping highly correlated
voxels; [25] used the supreme statistic in a random field to construct
confidence envelopes for the proportion of false discoveries and [1] used a
two-stage hierarchical testing procedure to test predefined clusters first
followed by a trimming stage to clean locations in which the signal is absent.
The second approach is to adaptively identify a collection of differentially
behaved regions with proven false discovery rate control. For example [30]
mapped the data in the wavelet domain first and removed redundant hypotheses
to reduce the number of hypotheses tested and improve power; [41] studied
multiple testing via false discovery rate control for large scale imaging
data; [34] treated each cluster as a testing unit and defined the false
discovery as the clusters that are falsely declared among all declared
clusters under the assumption that the number of false discoveries is
approximately Poisson. The Poisson approximation requires the sparsity
assumption on the signals. [9] gave a summary of literature in this area and
developed new tools for spatial multiple testing. In this line of research, a
cluster is defined to be a true discovery if it has non-zero overlap with the
support of the signal. Methods that try to incorporate cluster size to improve
power were also explored in [9].

In this paper, we study multiple testing problems for data with clustered
signals. We propose a new test statistic that adaptively recognizes such
clusters. Our test statistic aggregates information along a sliding window to
boost signal noise ratio. At the boundary between signal and non-signal
segments, the test statistic can be much larger than it is within the non-
signal cluster. We investigate the asymptotic distribution of the proposed
test statistic and set up rejection criterion controlling certain type I
errors. A new algorithm is proposed to locate signal clusters for followup
studies. We do not require signals to be sparse, which may be especially
valuable given the current conjecture of polygenic effects on complex disease
[42]. Furthermore, we allow signals to vary within a cluster, which differs
from the popular assumption that the means are identical within the same
cluster [39, 40]. Numerical studies show that when signals have varied sizes,
the proposed method has increased detection accuracy compared to method that
assumes same signal size within a cluster [39, 37, 38]. Computationally, the
speed of our algorithm is linear with number of tests while the algorithm used
in [39] is quadratic. Unlike [6], we present a new approach for variance
estimation under the setup of multiple testing. This is accomplished through
the order statistics of the average squares of the original data across a
sliding window, which is consistent under certain regularity conditions for
the one realization case. In addition, we consider the multiple realization
scenario and allow the variance to be heteroscedastic. New test statistics are
proposed with unknown parameters consistently estimated with available data to
conduct statistical inference. Moreover, the newly proposed algorithms are
more accurate in detecting break-points than algorithms proposed in [6] as an
additional turning parameter $\delta$ is used in the maximization to locate
break-points. Numerical studies show improved detection precision compared to
methods that did not utilize the clustering structure [12].

Recent multiple testing procedures that incorporate covariates require
estimation of the prior probability that the $i$th test corresponds to a null,
$i=1,\ldots,m.$ These weights are then estimated adaptively from available
data. In particular, [35] uses an empirical Bayesian two group mixture model
and proposes to minimize a penalized likelihood function where fused lasso
type of penalty is used to have spatial smoothing [36]. OrderShapeEM proposed
by [7] imposes a monotone increasing constraint on the prior probability of
being null and a monotone decreasing constraint on the density function of
$p$-value under alternative distribution. The implementation is achieved
through combination of EM algorithm and pool-adjacent-violator-algorithm
(PAVA). AdaPT [18] requires an order of the $p$-value to incorporate external
information to boost power. SABHA [20] modifies the BH procedure by
incorporating the probability that the $i$th test corresponds to a null,
$i=1,\ldots,m.$ SABHA further suggests different ways to estimate such
probabilities, including ordering, grouping, and low total variation. AdaPT
and SABHA achieves finite sample control of FDR.

In our work, we impose block signal structure to improve power. Different from
covariate adjusted multiple testing, we do not use covariate for individual
test, instead, we treat clustered signals through aggregation of individual
$p$-values. We do not require external covariate, such as ordering. Our
results are asymptotic in terms of number of test $m$.

An important method for spatial cluster detection is based on scan statistics
[26, 33]. Scan statistic is defined as the maximum number of points in a fixed
window as the window is shifted across the domain. The $p$-value is computed
under the uniform distribution on the domain and the threshold is designed to
control the familywise type I error. This statistic is used for an omnibus
test of the null hypothesis that there is no clustering. If the test rejects
the null hypothesis, then it leaves open the question of where and how much
clustering exists. Our test statistic is devised to compare the observed
information with its expected value under the null hypothesis that there is no
signal and then take the maximum across the domain. If the omnibus test
detects signals, our proposed algorithm can locate such signals which is of
special interest for followup studies.

The rest of the paper is organized as follows. In Section 2, we introduce a
structured hypothesis testing problem with one realization. Section 3 studies
the case that there are multiple realizations. In Section 4, we examine the
performance of the proposed procedure via simulation; we see that our
procedure is better able to detect clustered signals and the variance estimate
has a good performance. Section 5 presents an application of the methodology
to an array based Comparative Genomic Hybridization (aCGH) dataset.

##  2 Test and estimation with one realization

In this section we shall first present a structured hypothesis testing problem
with locally clustered signals. Suppose we are given noisy data of the form

$X_{j}=\mu_{j}+Z_{j},\,1\leq j\leq p,$ (2.1)

where $Z_{j}$ are i.i.d. with mean $0$ and variance $\sigma^{2}$, and
$\mu_{j}$ are means or signals. We say that a signal is present at location
$j$ if $\mu_{j}\not=0$. In the study of aCGH data, we let $X_{j}$ be the
$\log_{2}$ ratio between the test and the reference sample intensities at
locus $j$. Then $X_{j}>0$ (resp. $X_{j}<0$) means copy number duplication
(resp. deletion). In this section we assume that one realization
$(X_{j})_{j=1}^{p}$ is available. In Section 3 we shall deal with the
situation that multiple realizations are available with possibly non-i.i.d.
$Z_{j}$. Based on the observation $(X_{j})_{j=1}^{p}$, we test the null
hypothesis of no signal

$\displaystyle H_{0}:\,\mu_{1}=\cdots=\mu_{p}=0$ (2.2)

versus the alternative hypothesis that signals are clustered: there exist
break-points $1=\tau_{0}\leq\tau_{1}<\cdots<\tau_{l}\leq\tau_{l+1}=p$ such
that

$\displaystyle H_{1}:\mu_{1}=\cdots=\mu_{\tau_{1}-1}=0,$
$\displaystyle\mu_{\tau_{1}},\cdots,\mu_{\tau_{2}-1}\neq 0,$ (2.3)
$\displaystyle\mu_{\tau_{2}}=\cdots=\mu_{\tau_{3}-1}=0,$
$\displaystyle\mu_{\tau_{3}},\cdots,\mu_{\tau_{4}-1}\neq 0,\cdots.$ (2.4)

Let ${\cal S}_{f}=\\{\tau_{f},\ldots,\tau_{f+1}-1\\}$, $f=1,2,\ldots$. We call
sets ${\cal S}_{1},{\cal S}_{3},\ldots$, signal clusters on which $\mu_{j}$s
are non-zero and let ${\cal S}={\cal S}_{1}\cup{\cal S}_{3}\cup\ldots$ be the
signal set. Let ${\cal N}={\cal S}_{0}\cup{\cal S}_{2}\cup\ldots$ be the non-
signal set. Note that our definition of break-points is different from change-
points that are used in change-point analysis, where the alternative
hypothesis is typically formulated as

$\displaystyle
H_{c}:\mu_{1}=\cdots=\mu_{\tau_{1}-1}\not=\mu_{\tau_{1}}=\cdots=\mu_{\tau_{2}-1}\not=\mu_{\tau_{2}}=\cdots=\mu_{\tau_{3}-1}\not=\cdots.$

For example, if there exists a $j$ in the signal cluster ${\cal
S}_{1}=\\{\tau_{1},\ldots,\tau_{2}-1\\}$ of (2.3) such that
$\mu_{\tau_{1}}=\ldots=\mu_{j}\not=\mu_{j+1}=\ldots=\mu_{\tau_{2}-1}$, then
this $j$ is a change-point while it is not a break-point in our sense. While
providing a very general framework, our setting of allowing unequal $\mu_{j}$s
in the signal clusters substantially complicates the related statistical
inference. The primary goal of the paper is to test $H_{0}$ vs $H_{1}$ and to
locate those break-points.

###  2.1 One-sided Test

If in the signal sets ${\cal S}_{1},{\cal S}_{3},\ldots$, all non-zero
$\mu_{i}$ are positive, namely

$\displaystyle H^{\prime}_{1}:\mu_{1}=\cdots=\mu_{\tau_{1}-1}=0,$
$\displaystyle\mu_{\tau_{1}},\cdots,\mu_{\tau_{2}-1}>0,$ (2.5)
$\displaystyle\mu_{\tau_{2}}=\cdots=\mu_{\tau_{3}-1}=0,$
$\displaystyle\mu_{\tau_{3}},\cdots,\mu_{\tau_{4}-1}>0,\cdots,$ (2.6)

then we can use the following test statistic

$\displaystyle R_{i}^{\circ}=\frac{1}{k}\sum_{j=i+1}^{i+k}X_{j},$ (2.7)

where $k$ is the window size parameter. Note that the mean is
$ER_{i}^{\circ}=k^{-1}\sum_{j=i+1}^{i+k}\mu_{j}$. Intuitively, $i$ can be
classified in the signal cluster if $R_{i}^{\circ}$ is big. The cutoff values
can be computed based on Theorem 2.1, which provides a uniform Gaussian
approximation of the distribution of $R_{i}^{\circ}$. Theorem 2.1 follows from
Theorem 3.1 with $n=1$. For completeness we state it here. It asserts that
under $H_{0}$, $R_{i}^{\circ}$ can be uniformly approximated by the Gaussian
process

$\displaystyle\sigma
G^{\circ}_{i}=\frac{1}{k}\sum_{j=i+1}^{i+k}\sigma\eta_{j},\mbox{ where
}\eta_{j}\overset{i.i.d.}{\sim}N(0,1).$ (2.8)

We shall quantify the closeness by the coupling distance

$\displaystyle\Delta^{\circ}=\sqrt{k}\max_{0\leq j\leq
p-k}|R_{j}^{\circ}/\sigma-G^{\circ}_{j}|$ (2.9)

and the distributional distance

$\rho^{\circ}=\sup_{u}|P(\sqrt{k}\max_{0\leq j\leq
p-k}R_{j}^{\circ}/\sigma\geq u)-P(\sqrt{k}\max_{0\leq j\leq
p-k}G^{\circ}_{j}\geq u)|.$ (2.10)

We first introduce a moment condition.

######  Condition 2.1.

$Z_{1},Z_{2},\ldots$, are i.i.d. with mean $0$ and variance $\sigma^{2}$, and
the $\theta$th norm
$\|Z_{i}\|_{\theta}:=(E|Z_{i}|^{\theta})^{1/\theta}<\infty$, where $\theta>2$.
Write $K_{\theta}:=\|Z_{i}\|_{\theta}$.

######  Theorem 2.1.

Assume Condition 2.1 and $\mu_{i}=0,1\leq i\leq p.$ (i) Let $\theta>2$. Then
there exists a possibly larger probability space on which one can define
$(Z_{j})_{j}$ and $(\eta_{j})_{j}$ such that, for all $u>0$ and any positive
integer $k$,

$\displaystyle P\left[k^{1/2}\Delta^{\circ}\geq
c_{0}u\right]\leq{\frac{pK_{\theta}^{\theta}}{u^{\theta}\sigma^{\theta}}},$
(2.11)

where $c_{0}$ is a constant only depending on $\theta$. (ii) Let $\theta>3$.
The distributional distance

$\rho^{\circ}\lesssim{k^{-1/6}(\log
p)^{7/6}+(pk^{-\theta/2})^{1/(\theta+1)}(\log p)^{(3\theta-2)/(2+2\theta)}}$
(2.12)

where $a\lesssim b$ means a = O(b) and the multiplicative constant in
$\lesssim$ only depends on $\theta$, $\sigma^{2}$ and $K_{\theta}$. Namely
there exists a constant $C>0$ depending on $\theta$, $\sigma^{2}$ and
$K_{\theta}$ such that $\rho^{\circ}\leq C(k^{-1/6}(\log
p)^{7/6}+(pk^{-\theta/2})^{1/(\theta+1)}(\log p)^{(3\theta-2)/(2+2\theta)}).$

Theorem 2.1 implies that, if the window size $k$ satisfies
$p^{2/\theta}=o(k)$, then $\Delta^{\circ}=o_{P}(1)$ by letting
$u=(pk)^{1/(\theta+2)}$. Under the slightly stronger condition

$\displaystyle p^{2/\theta}(\log p)^{3-2/\theta}=o(k),$ (2.13)

we have $\rho^{\circ}=o(1)$, suggesting that $R_{j}^{\circ}$ and $\sigma
G^{\circ}_{j}$ are uniformly close.

Let $\hat{\sigma}^{2}$ be an estimate of $\sigma^{2}$ and $g_{1-\alpha}$ be
the $(1-\alpha)$th quantile of $\max_{0\leq j\leq p-k}G_{j}^{\circ}$,
$\alpha\in(0,1)$. The latter can be computed by Monte Carlo simulations. In
Section 2.4.1 we shall propose a consistent estimate of $\sigma^{2}$ when
$(\mu_{j})_{j}$ has form (2.3). Theorem 2.1 suggests rejecting $H_{0}$ and
accepting the alternative hypothesis $H^{\prime}_{1}$ of (2.5) at level
$\alpha$ if $\max_{0\leq j\leq p-k}R^{\circ}_{j}>\hat{\sigma}g_{1-\alpha}$.
Alternatively, let $T=p/k$, by Corollary A1 in [3] we can also have the Gumbel
convergence

$P[{{\max_{0\leq j\leq p-k}\sqrt{k}G^{\circ}_{j}}\over\sqrt{2\log
T}}-1-{{\log\log T-{1\over 2}\log(4\pi)}\over{4\log T}}\leq v]\to
e^{-e^{-v}},$ (2.14)

which gives an approximate solution for $g_{1-\alpha}$ by letting
$v=-\log\log(1-\alpha)^{-1}$. We do not recommend the latter since the
convergence of (2.14) is very slow. A bootstrap calibration procedure is
proposed in Section 3.2 which has better finite sample properties.

###  2.2 Two-sided Test

Under the general alternative $H_{1}$ of (2.3), the test statistic (2.7) is no
longer applicable since the $\mu_{j}$s in the signal clusters can potentially
cancel each other out. As a simple remedy, assuming at the outset that
$\sigma^{2}$ is known, we define the modified version

$\displaystyle
R_{i}^{\dagger}=\frac{1}{k}\sum_{j=i+1}^{i+k}(X_{j}^{2}-\sigma^{2}),$ (2.15)

which, since $\epsilon_{j}=Z_{j}^{2}-\sigma^{2}+2\mu_{j}Z_{j}$ has mean $0$
under $H_{0}$, mimics $R_{i}^{\circ}$ in (2.7) in view of

$\displaystyle
X_{j}^{2}-\sigma^{2}=\mu_{j}^{2}+(Z_{j}^{2}-\sigma^{2}+2\mu_{j}Z_{j})=\mu_{j}^{2}+\epsilon_{j}.$
(2.16)

Hence a location $i$ with a big value of $R_{i}^{\dagger}$ will likely be in
signal clusters, regardless of signs of $\mu_{j}$. Then we can apply the one-
sided test procedure in Section 2.1. Note that the other modified version
$R_{i}^{\star}:=k^{-1}\sum_{j=i+1}^{i+k}(|X_{j}|-m_{1})$, where
$m_{1}=E|Z_{j}|$, does not have the property that
$E(|X_{j}|-m_{1})=E(|\mu_{j}+Z_{j}|-|Z_{j}|)>0$ for non-zero $\mu_{j}$. So in
general $R_{i}^{\star}$ cannot be used in the two-sided test. Assume that
$E(Z^{4}_{j})<\infty$. Similar to (2.9) and (2.10), we define

$\displaystyle\Delta^{\dagger}=\sqrt{k}\max_{0\leq j\leq
p-k}|R_{j}^{\dagger}/\kappa-G^{\circ}_{j}|,\mbox{ where
}\kappa=\|Z_{j}^{2}-\sigma^{2}\|_{2}=[E(Z_{j}^{2}-\sigma^{2})^{2}]^{1/2},$
(2.17)

and the distributional distance

$\displaystyle\rho^{\dagger}=\sup_{u}|P(\sqrt{k}\max_{0\leq j\leq
p-k}R_{j}^{\dagger}/\kappa\geq u)-P(\sqrt{k}\max_{0\leq j\leq
p-k}G_{j}^{\dagger}\geq u)|.$

Note that under $\mu_{j}=0,$ we have
$\mbox{Var}(\epsilon_{j})=E(Z_{j}^{2}-\sigma^{2}+2\mu_{j}Z_{j})^{2}=\kappa.$

######  Corollary 2.1.

Assume Condition 2.1 hold with $\theta>4$ and $\mu_{i}=0,1\leq i\leq p.$ Then
there exists a larger probability space on which one can define $(Z_{j})_{j}$
and $(\eta_{j})_{j}$ such that for all $u>0$,

$\displaystyle P\left[k^{1/2}\Delta^{\dagger}\geq
c_{0}u\right]\leq{\frac{pK_{\theta}^{\theta}}{(\kappa u)^{\theta/2}}},$ (2.18)

where $c_{0}$ is a constant only depending on $\theta$, and the distributional
distance

$\rho^{\dagger}\lesssim{k^{-1/6}(\log
p)^{7/6}+(pk^{-\theta/4})^{2/(\theta+2)}(\log p)^{(3\theta-4)/(4+2\theta)},}$
(2.19)

where the constant in $\lesssim$ only depends on $\theta$, $\kappa$ and
$K_{\theta}$.

Corollary 2.1 follows from Theorem 2.1 by replacing $\theta$ in the latter by
$\theta/2$ in view of (2.16). In comparison with (2.13), Corollary 2.1
requires the stronger condition $p^{4/\theta}(\log p)^{3-4/\theta}=o(k)$ to
ensure that $\rho^{\dagger}=o(1)$.

Estimation of $\sigma^{2}$ and $\kappa$ is discussed in Section 2.4.1. Recall
that $g_{1-\alpha}$ is the $(1-\alpha)$th quantile of $\max_{0\leq j\leq
p-k}G_{j}^{\circ}$, $\alpha\in(0,1)$. Corollary 2.1 suggests rejecting $H_{0}$
and accepting the alternative hypothesis $H_{1}$ of (2.5) at level $\alpha$ if
$\max_{0\leq j\leq p-k}R_{j}^{\dagger}>\hat{\kappa}g_{1-\alpha}.$

######  Remark 1.

Denote by $\Delta^{\circ}_{k}$ the quantity $\Delta^{\circ}$ in (2.9). A
careful analysis of the proof of Theorem 3.1 (which implies Theorem 2.1 with
$n=1$) indicates that Theorem 2.1 is still valid with $\Delta^{\circ}_{k}$
(resp. $R_{k}^{\bullet}:=\sqrt{k}\max_{0\leq j\leq p-k}R_{j}^{\circ}/\sigma$
and $G_{k}^{\bullet}:=\sqrt{k}\max_{0\leq j\leq p-k}G^{\circ}_{j}$) therein
replaced by the uniform version $\max_{k\leq m\leq p}\Delta^{\circ}_{m}$
(resp. $\max_{k\leq m\leq p}R_{m}^{\bullet}$ and $\max_{k\leq m\leq
p}G_{m}^{\bullet}$). Similarly, for the two-sided test, Corollary 2.1 also
holds with the uniform version $\max_{k\leq m\leq p}m^{-1/2}\max_{0\leq j\leq
p-m}$
$\sum_{i=j+1}^{j+m}(X_{i}^{2}-\sigma^{2})$. The latter quantity has an
interesting connection with the adaptive Neyman’s high dimensional
multivariate normal mean test which has the form $\max_{1\leq m\leq
p}(2m)^{-1/2}\sum_{i=1}^{m}$
$(X_{i}^{2}-\sigma^{2})$, which was considered in Section 2.1 in [11] in the
setting that large values of $\mu$ concentrate on the first $m$ dimensions and
$X_{i}\overset{i.i.d.}{\sim}N(\mu_{i},1)$. Here $m$ is estimated by the
maximizer $\hat{m}={\rm argmax}_{1\leq m\leq
p}(2m)^{-1/2}\sum_{i=1}^{m}(X_{i}^{2}-\sigma^{2})$.

###  2.3 An Algorithm for Locating Break-points

Once the null hypothesis is rejected, we need to locate break-points. We
propose Algorithms 2.1 and 2.2 for locating break-points based on the one- and
the two-sided tests, respectively. Theoretical properties of Algorithm 2.1
(resp. 2.2) are given in Theorem 2.2 (resp. 2.3).

####  2.3.1 Locating break-points based on one-sided test

We first present an algorithm based on the one-sided test.

######  Algorithm 2.1.

Step 1. Let $L_{j}^{\circ}=R_{j-k}^{\circ},j=k,\ldots,p$. Compute
$Q_{j}^{\circ}=1(R_{j}^{\circ}>\gamma)+1(L_{j}^{\circ}>\gamma)$ for a pre-
specified cutoff value $\gamma$, $j=k,\ldots,p-k.$ We use a majority vote
approach to smooth $Q_{j}^{\circ}.$ Specifically, denote
$j_{0}=\sum_{j=i-k}^{i+k}I\\{Q_{j}^{\circ}=0\\},j_{1}=\sum_{j=i-k}^{i+k}I\\{Q_{j}^{\circ}=1\\},$
and $j_{2}=\sum_{j=i-k}^{i+k}I\\{Q_{j}^{\circ}=2\\}.$ Let
$\tilde{Q}_{j}^{\circ}=\\{k,\mbox{such that}\
j_{k}=\max_{l\in\\{0,1,2\\}}j_{l}\\}.$

Step 2. Decompose $\\{1,\ldots,p\\}=W_{0}\cup W_{1}\cup W_{2},$ where $j\in
W_{0}$ if $\tilde{Q}_{j}^{\circ}=0,j\in W_{1}$ if $\tilde{Q}_{j}^{\circ}=1$
and $j\in W_{2}$ if $\tilde{Q}_{j}^{\circ}=2.$ Let ${\cal M}_{1},\ldots,{\cal
M}_{\hat{l}}$ be connected components of $W_{1}.$

Step 3. Given $\delta<\gamma$, the break-points are defined as
$\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{R_{j}^{\circ}:\,L_{j}^{\circ}\leq\delta\\}$ if ${\cal M}_{i}$ is the
transition region from $W_{0}$ to $W_{2}$. If ${\cal M}_{i}$ is the transition
region from $W_{2}$ to $W_{0},$ $\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{L_{j}^{\circ}:\,R_{j}^{\circ}\leq\delta\\}$.

The estimated signal sets are $\hat{\cal
S}_{1}=\\{\hat{\tau}_{1},\ldots,\hat{\tau}_{2}-1\\},\hat{\cal
S}_{3}=\\{\hat{\tau}_{3},\ldots,\hat{\tau}_{4}-1\\},\ldots$. The rationale
behind Algorithm 2.1 is that if $\mu_{j}=0$, then $R_{j}^{\circ}$ is close to
$0;$ on the other hand, in the signal clusters, $R_{j}^{\circ}$ tends to be
large. By locally averaging the data, we can reduce the variability, which has
the effect of boosting the signal noise ratio. If there are many weak signals,
we are able to detect them by the aggregation. On the other hand, if sporadic
large values of $X_{j}$ arise, they can be smoothed out through
$R_{j}^{\circ}$ to avoid false discoveries. Therefore, we can effectively de-
noise the data to achieve better inference. In the signal cluster,
$Q_{j}^{\circ}$ is most likely to be $2$ and in the non-signal cluster,
$Q_{j}^{\circ}$ is most likely to be $0.$ In the boundary between signal and
non-signal cluster, $Q_{j}^{\circ}$ is most likely to be $1.$ After Step 1, we
get smoothed $\tilde{Q}_{j}^{\circ}$ that are in clusters of $0,1$ and $2.$
Step 2 focuses on the signal and non-signal cluster boundary regions, where
$\tilde{Q}_{j}^{\circ}=1.$ Step 3 locates break-points. The basic idea is that
without noise at the true break-points, $R_{j}^{\circ}$ reaches the maximum as
there is no noise to dilute the summation if we are transiting from non-signal
cluster to signal cluster. The constraint $L_{j}^{\circ}\leq\delta$ prevents
the detected break-points to be too far from the true break-points when
$\mu_{j}$ increases in the signal cluster. Similarly, without noise,
$L_{j}^{\circ}$ obtains the maximum if we are transiting from signal to non-
signal cluster at the true break-points. The constraint
$R_{j}^{\circ}\leq\delta$ prevents the detected break-points to be too far
from the true break-points when $\mu_{j}$ decreases in the signal cluster.
With two thresholds $\delta<\gamma$, Algorithm 2.1 has more flexibility and
produces more accurate estimates of the break-points than the procedure in [6]
which only uses one threshold $\gamma$.

Our method depends on the choice of window size $k$ and thresholds $\gamma$
and $\delta$. Theoretically speaking, the allowable range of $k$ is specified
in (2.13). Our simulation studies show that the proposed method is relatively
robust to different choices of $k.$ In practice, following the idea of the
adaptive Neyman’s high dimensional multivariate normal mean test mentioned in
Remark 1, as a simple rule of thumb choice we can let $\hat{m}=\mathop{\rm
arg\max}_{m\geq\sqrt{p}}R_{m}^{\bullet}$ and $k=\lfloor\hat{m}/2\rfloor$. For
a data-driven selection of $\gamma$ and $\delta$, we can choose
$\gamma=\hat{\sigma}g_{1-\alpha}$ and $\delta=\hat{\sigma}g_{1,1-\alpha}$,
where $g_{1-\alpha}$ and $g_{1,1-\alpha}$ are the $(1-\alpha)$th quantiles of
$\mbox{max}_{0\leq j\leq p-k}G_{j}^{\circ}$ and $\mbox{max}_{j\in
W_{1}}G_{j}^{\circ}$, respectively, with $\alpha$ close to $0$. They can be
obtained by simulations. Section 2.4.1 gives an estimate $\hat{\sigma}$ of
$\sigma$.

######  Condition 2.2.

Recall ${\cal S}_{f}=\\{\tau_{f},\ldots,\tau_{f+1}-1\\}$ and ${\cal S}={\cal
S}_{1}\cup{\cal S}_{3}\cup\ldots$ is the signal set. Assume
$d:=\min_{i\in{\cal S}}\mu_{i}>0$ and $2k<\min_{f}(\tau_{1+f}-\tau_{f})$.

######  Condition 2.3.

We say that a random variable $Z$ is $\sigma^{2}$-sub-Gaussian if
$E(\exp(uZ/\sigma))\leq\exp(u^{2}/2)$ for all $u\in{\mathbb{R}}$. Note that
$N(0,\sigma^{2})$ is $\sigma^{2}$-sub-Gaussian.

To state Theorem 2.2, we need to introduce truncated moment functions. For a
random variable $X$ with $E(X^{2})<\infty$, define the truncated moment

${\cal M}_{\upsilon}(X)=E\min(|X|^{\upsilon},\,X^{2})<\infty,\,\,\upsilon>2.$
(2.20)

If $X$ has finite $\theta$th moment with $2<\theta<\upsilon$, then ${\cal
M}_{\upsilon}(X)\leq E(|X|^{\theta})$. Theorem 2.2(i) (resp. (ii)) concerns
sub-Gaussian (resp. polynomial-tailed) noises.

######  Theorem 2.2.

Assume Condition 2.2 and $d/2>\gamma>\delta$. (i) Assume Condition 2.3 holds
for $Z_{j}$. Denote by $\hat{l}$ the estimated number of break points. Then

$\displaystyle 1-P\left[\hat{l}=l,\,\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq{2k\delta\over d}\right]\leq c_{3}({p\over
k}e^{-c_{1}k\gamma^{2}/\sigma^{2}}+le^{-c_{2}k\delta^{2}/\sigma^{2}}),$ (2.21)

where $c_{1},c_{2},c_{3}$ are absolute constants. (ii) Assume Conditions 2.1
and let $\upsilon\geq\theta$. Then

$\displaystyle 1-P\left[\hat{l}=l,\,\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq{2k\delta\over d}\right]$
$\displaystyle\lesssim$ $\displaystyle p{\cal
M}_{\upsilon}(Z_{1}/(k\gamma))+{p\over k}e^{-c_{1}k\gamma^{2}/\sigma^{2}}$
(2.23) $\displaystyle+lk{\cal
M}_{\upsilon}(Z_{1}/(k\delta))+le^{-c_{2}k\delta^{2}/\sigma^{2}}$
$\displaystyle\leq$ $\displaystyle
K_{\theta}^{\theta}(p\gamma^{-\theta}+lk\delta^{-\theta})k^{-\theta}$ (2.25)
$\displaystyle+pk^{-1}e^{-c_{1}k\gamma^{2}/\sigma^{2}}+le^{-c_{2}k\delta^{2}/\sigma^{2}},$

where $c_{1}$ and $c_{2}$ are absolute constants and the constant in
$\lesssim$ only depends on $\theta$ and $\upsilon$.

Theorem 2.2 is proved in the Supplementary Material. In comparison with
(2.21), the extra term
$K_{\theta}^{\theta}(p\gamma^{-\theta}+lk\delta^{-\theta})k^{-\theta}$ in
(2.23) is due to polynomial tails, which are heavier than the sub-Gaussian
ones. In the sub-Gaussian case (i) with unbounded $l$ (namely $l\to\infty$),
choose $\gamma=C_{1}(k^{-1}\log p)^{1/2}$, and $\delta=C_{2}(k^{-1}\log
l)^{1/2}$, where $C_{1}$ and $C_{2}$ are sufficiently large constants, we have
the uniform bound $\max_{j\leq l}|\hat{\tau}_{j}-\tau_{j}|={O_{P}(d^{-1}(k\log
l)^{1/2})}$. The condition $d/2>\gamma$ requires that $k\geq C_{3}d^{-2}\log
p$ for a sufficiently large constant $C_{3}$. When $l$ is bounded, by letting
$k=\lfloor Cd^{-2}\log p\rfloor$ for a sufficiently large $C$, we can
similarly obtain the uniform bound $\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|={O_{P}(d^{-2}(\log p)^{1/2})}$. In the context of
detecting a deterministic signal with unknown spatial extent in the univariate
sampled data model with standard white Gaussian noises, [8] dealt with the
special case $\mu_{j}=d{\bf 1}_{\tau_{1}\leq j<\tau_{2}}$ and considered the
consistency of detection based on scan statistics under the condition
$\tau_{2}-\tau_{1}\geq c_{p}d^{-2}\log p$, where $c_{p}=2+\iota_{p}$ and
$\iota_{p}^{2}\log p\to\infty$. The latter observation has a similar flavor as
our condition $k\geq C_{3}d^{-2}\log p.$

The polynomial-tailed case (2.23) is more involved. To ensure that the right
hand side of (2.23) goes to $0$, we can choose
$\gamma=C_{1}(k^{-1}p^{1/\theta}+(k^{-1}\log p)^{1/2})$ and
$\delta=C_{2}(k^{-1}(lk)^{1/\theta}+(k^{-1}\log l)^{1/2})$, where
$C_{1},C_{2}>0$ are sufficiently large constants. Assume $k\geq
C_{3}(d^{-2}\log p+d^{-1}p^{1/\theta})$ for a sufficiently large constant
$C_{3}$, we have the uniform consistency $\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq O_{P}(k\delta/d)$. Thus the numbers of false
discoveries and missed discoveries are bounded by $lO_{P}(k\delta/d)$. If $l$
is bounded, then the latter bound becomes $O_{P}[(k\log p)^{1/2})/d]$.

####  2.3.2 Locating break-points based on two-sided test

We next present an algorithm based on the two-sided test. It is similar to
Algorithm 2.1. With the square form (2.15), it can pick up signals with
alternating positive and negative signs. Same simulation assisted choice of
$\gamma$ and $\delta$ as in the one-sided test can be used.

######  Algorithm 2.2.

Step 1: Calculate $R_{i}^{\dagger}$ and let
$L_{i}^{\dagger}=R_{i-k}^{\dagger},i=k,\ldots,p.$ For a pre-specified
$\gamma,$ let
$Q_{i}^{\dagger}=1(R_{i}^{\dagger}>\gamma)+1(L_{i}^{\dagger}>\gamma),i=k,\ldots,p-k.$
The same majority vote approach as in Algorithm 2.1 is used to smooth
$Q_{i}^{\dagger},$ denoted as $\tilde{Q}_{i}^{\dagger}.$

Step 2: Decompose $\\{1,\ldots,p\\}=W_{0}\cup W_{1}\cup W_{2},$ where $i\in
W_{0}$ if $\tilde{Q}_{i}^{\dagger}=0,i\in W_{1}$ if
$\tilde{Q}_{i}^{\dagger}=1$ and $i\in W_{2}$ if $\tilde{Q}_{i}^{\dagger}=2.$
Let ${\cal M}_{1},\ldots,{\cal M}_{\hat{l}}$ be connected components of
$W_{1}.$

Step 3. Given $\delta<\gamma$, the break-points are estimated as
$\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{R_{j}^{\dagger}:\,L_{j}^{\dagger}\leq\delta\\}$ if ${\cal M}_{i}$ is
the transition region from $W_{0}$ to $W_{2}$. If ${\cal M}_{i}$ is the
transition region from $W_{2}$ to $W_{0},$
$\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{L_{j}^{\dagger}:\,R_{j}^{\dagger}\leq\delta\\}$.

######  Condition 2.4.

Recall Condition 2.2 for ${\cal S}$. Let $d=\min_{i\in{\cal S}}|\mu_{i}|>0$
and assume $2k<\min_{f}(\tau_{1+f}-\tau_{f})$.

######  Theorem 2.3.

Assume Conditions 2.3, 2.4, and $(k^{-1}\log p)^{1/2}=o(d^{2})$. Let
$\gamma=c_{1}(k^{-1}\log p)^{1/2}$ and $\delta=c_{2}(k^{-1}\log l)^{1/2}$,
where $c_{1}$ and $c_{2}$ are sufficiently large constants. Then there exists
a constant $c>0$ independent of $k$ and $p$ such that

$\displaystyle P\left[\hat{l}=l,\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq{ck\delta\over d^{2}}\right]\to 1.$ (2.26)

Theorem 2.3 provides a bound for uniform deviations of the estimated break-
points. It is proved in the Supplementary Material, where the polynomial-
tailed case is also studied. Same choice of $\gamma$ and $\delta$ can be used
as in the one-sided test scenario.

###  2.4 Variance Estimation

####  2.4.1 Estimation of $\sigma^{2}$

To apply Theorem 2.1 and Corollary 2.1 for computing the cutoff values based
on $R_{j}^{\circ}$ and $R_{j}^{\dagger}$, we need to deal with the key issue
of estimating the variance $\sigma^{2}$. Furthermore, to use
$R_{j}^{\dagger}$, we need to estimate $\kappa^{2}$. For the nonparametric
regression model $X_{i}=\mu_{i}+Z_{i}=f(i/p)+Z_{i},\,1\leq i\leq p$, where
$\mu_{i}=f(i/p)$, $f$ is a smooth function and $Z_{i}$ are i.i.d. with mean
$0$ and variance $\sigma^{2}$, the problem of estimating $\sigma^{2}$ has a
long history; see [14] and references therein. However the difference-based
method in the latter paper does not work here. Due to the presence of the
nonzero $\mu_{j}$s, the problem of estimating $\sigma^{2}$ is highly
nontrivial. The latter problem is further complicated by the fact that the
nonzero $\mu_{j}$s in the signal segments can change wildly. Here we shall use
order statistics and obtain a consistent estimator. Let

$\hat{\sigma}^{2}_{i}=\frac{1}{m}\sum_{j=i}^{i+m-1}X_{j}^{2},\quad 1\leq i\leq
p^{\prime},\mbox{ where }p^{\prime}=p-m+1.$ (2.27)

Let
$\hat{\sigma}^{2}_{(1)}\leq\hat{\sigma}^{2}_{(2)}\leq\ldots\leq\hat{\sigma}^{2}_{(p^{\prime})}$
be the order statistics of
$\hat{\sigma}^{2}_{1},\ldots,\hat{\sigma}^{2}_{p^{\prime}}$. Theorem 2.4 shows
that, for any $k\leq p^{\prime}/2$, $\hat{\sigma}^{2}_{(k)}$ is a consistent
estimator of $\sigma^{2}$ under suitable conditions of $m$. The intuition is
as follows: for large $m$, we expect that $\hat{\sigma}^{2}_{i}\approx
E\hat{\sigma}^{2}_{i}=\sigma^{2}+m^{-1}\sum_{j=i}^{i+m-1}\mu_{j}^{2}$. The
latter uniform closeness relation will be made rigorous in the proof of
Theorem 2.4, which is proved in the Supplementary Material. Under Condition
2.5 below, we expect that majority of $\sum_{j=i}^{i+m-1}\mu_{j}^{2}$ will be
$0.$ Thus the median or any lower quantile of $E\hat{\sigma}^{2}_{i}$ is
$\sigma^{2}$.

In practice we choose the sample median estimate with $k=p^{\prime}/2$.

######  Condition 2.5.

There exists a constant $c>0$ such that the length of non-signal clusters
$\tau_{i+1}-\tau_{i}\geq cp$ for all even $i$, and the total length ${\sum_{i\
\mbox{is even}}}\ (\tau_{1+i}-\tau_{i})\geq\lambda p$ with constant
$\lambda>1/2$.

Condition 2.5 implies the natural requirement that the proportion of non-
signals (namely $j$ with $\mu_{j}\neq 0$) is larger than $1/2$.

######  Theorem 2.4.

Assume $(\mu_{j})$ satisfies (2.3), Condition 2.5, $Z_{i}\in{\cal
L}^{\theta},\theta>4$, $p^{2/\theta}=o(m)$ and $m=o(p)$. Then we have for any
$k\leq p^{\prime}/2$ that

$\displaystyle\hat{\sigma}^{2}_{(k)}=\sigma^{2}+O_{P}(\gamma_{p}),\mbox{ where
}\gamma_{p}=({\log p\over m})^{1/2}+{p^{2/\theta}\over m}.$ (2.28)

If Condition 2.3 holds, $\log p=o(m)$ and $m=o(p)$, then
$\hat{\sigma}^{2}_{(k)}=\sigma^{2}+O_{P}(({m^{-1}\log p})^{1/2})$.

####  2.4.2 Estimation of $\kappa$

The estimation of $\kappa$ in (2.17) is much more involved. The key issue is
to estimate the fourth order moment $E(Z_{i}^{4})$. Unlike (2.27), we cannot
simply use order statistics of the moving window sample averages
$\Xi_{i}:=m^{-1}\sum_{j=i}^{i+m-1}X_{j}^{4},1\leq i\leq p-m+1,$ to estimate
$E(Z_{i}^{4})$, since the median or lower quantile of
$E(\Xi_{i})=E(Z_{i}^{4})+m^{-1}\sum_{j=i}^{i+m-1}(\mu_{j}^{4}+6\mu_{j}^{2}\sigma^{2}+4\mu_{j}^{3}E(Z_{j}^{3}))$
is generally not $E(Z_{i}^{4})$ if $E(Z_{j}^{3})\not=0$. $E(\Xi_{i})$ can be
greater or less than $E(Z_{i}^{4})$ depending on what
$\mu_{j},j=i,\ldots,i+m-1$ and $E(Z_{j}^{3})$ are. The reason is that the
function $E(\mu+Z_{j})^{4}$ may not be minimized at $\mu=0$. For example, if
$Z_{j}=E_{j}-1$ with $E_{j}\sim\exp(1)$, then $E(\mu+Z_{j})^{4}$ is minimized
at $\mu\approx-0.596072$. To circumvent the latter problem, we introduce

$\displaystyle\hat{\nu}_{i}={1\over
m}\sum_{j=i}^{i+m-1}(X_{j}-X_{j-1})^{4},\quad 2\leq i\leq p-m+1,$ (2.29)

and $\nu=E(Z_{1}-Z_{0})^{4}=2\kappa^{2}+8\sigma^{4}$. Note that
$E(\mu+Z_{1}-Z_{0})^{4}$ is indeed minimized at $\mu=0$. The above estimate
resembles the first order difference based estimate; see [14]. However the
setting and the motivation are quite different. Let
$\hat{\nu}_{(2)}\leq\ldots\leq\hat{\nu}_{(p-m+1)}$ be the order statistics.
Corollary 2.2 below concerns asymptotics for $\hat{\nu}_{(k)}$. It is proved
in the Supplementary Material. Then we can estimate $\kappa^{2}$ by
$\hat{\kappa}^{2}=\hat{\nu}_{(k)}/2-4\hat{\sigma}^{4}_{(k)}$. In practice we
can choose $k=p^{\prime}/2$, which corresponds to the sample median. By
Theorems 2.4 and Corollary 2.2, we have
$\hat{\kappa}^{2}=\kappa^{2}+O_{P}({\phi_{p,m}})$, where $\phi_{p,m}$ is a
function of $p$ and $m,$ given in the following corollary.

######  Corollary 2.2.

Assume (2.3), Condition 2.5 and that $Z_{i}\in{\cal L}^{q},q>4$,
$p^{4/q}=o(m)$ and $m=o(p)$. Then we have for any $k\leq p^{\prime}/2$ that

$\displaystyle\hat{\nu}^{2}_{(k)}=\nu+O_{P}({\phi_{p,m}}),\mbox{ where
}{\phi_{p,m}}=({\log p\over m})^{1/2}+{p^{4/q}\over m}.$ (2.30)

If Condition 2.3 holds, $(\log p)^{2}=o(m)$ and $m=o(p)$, then
$\hat{\nu}^{2}_{(k)}=\nu+O_{P}(({m^{-1}\log p})^{1/2})$.

##  3 Test and estimation with multiple realizations

In Section 2, only one realization $(X_{j})_{j=1}^{p}$ is available, under the
assumption that the errors $Z_{j}$ are i.i.d. When we have more than one
realization, we will be able to detect clustered signals even if the variances
change along the sequence. The allowance of heteroscedasticity substantially
expands the application of our methods. Let $n({\geq 2})$-realizations
$Y_{i}=(Y_{i1},\ldots,Y_{ip})^{T}$ be observed, $i=1,\ldots,n$, with

$\displaystyle Y_{ij}=\mu_{j}+Z_{ij},\,\,1\leq j\leq p,$ (3.1)

where $Z_{ij}$ has mean $0$, variance $\sigma_{j}^{2}$ and independent across
both $i$ and $j$. We are interested in testing (2.5) and (2.3). To this end,
we propose a new test statistic and derive an omnibus test under the global
null hypothesis $H_{0}$ in (2.2). Let
$\hat{\mu}_{j}=n^{-1}\sum_{i=1}^{n}Y_{ij}$.

###  3.1 One-sided Test

Given a window size $k$, define

$R_{j}^{\star}={{\sum_{l=j+1}^{j+k}\sqrt{n}\hat{\mu}_{l}}\over{v_{j}^{1/2}}},\,\mbox{
where }v_{j}=\sum_{l=j+1}^{j+k}\sigma^{2}_{l},\,0\leq j\leq p-k.$ (3.2)

Let $(G^{\star}_{j})_{0\leq j\leq p-k}$ be a mean zero Gaussian vector which
has the same covariance structure as $(R^{\star}_{j})_{0\leq j\leq p-k}$. As a
stochastic realization, we can let

$G^{\star}_{j}={W_{j}\over{v_{j}^{1/2}}},\mbox{ where
}W_{j}=\sum_{l=j+1}^{j+k}\sigma_{l}\eta_{l},\,{v_{j}=E(W_{j}^{2})}\ \mbox{ and
}\eta_{l}\overset{i.i.d.}{\sim}N(0,1).$ (3.3)

Let $\sigma=(\sigma_{1},\ldots,\sigma_{p})$. Then $G^{\star}_{j}$ has marginal
variance $1$ and covariance matrix
$\Gamma(\sigma)=(\gamma_{j,j^{\prime}}(\sigma))_{0\leq j,j^{\prime}\leq p-k}$
with
$\gamma_{j,j^{\prime}}(\sigma)=v^{-1/2}_{j}v^{-1/2}_{j^{\prime}}E(W_{j}W_{j^{\prime}}).$
Note that $\gamma_{j,j^{\prime}}(\sigma)=0$ if $|j-j^{\prime}|\geq k$ and
$(W_{j})$ are $(k-1)$-dependent. Let the coupled distance

$\displaystyle\Delta^{\star}=\max_{0\leq j\leq
p-k}|R^{\star}_{j}-G^{\star}_{j}|.$ (3.4)

Theorem 3.1 below concerns the Gaussian approximation in terms of the
closeness of $R_{j}^{\star}$ and $G_{j}^{\star}$ with various metrics. It is
proved in the Supplementary Material. Relation (3.6) is a coupling statement
which provides a tail probability inequality for the maximum distance
$\Delta^{\star}$ on some common probability space, while (3.7) is for the
distributional distance

$\displaystyle\rho^{\star}:=\sup_{u}|P(\max_{0\leq j\leq p-k}R^{\star}_{j}\geq
u)-P(\max_{0\leq j\leq p-k}G^{\star}_{j}\geq u)|.$ (3.5)

We shall impose the following regularity condition.

######  Condition 3.1.

Let $\theta>2$. Assume that there exist positive constants $\sigma_{*}$,
$\sigma^{*}$ and $K_{\theta}$ such that, for all $1\leq j\leq p$,
$\sigma_{*}\leq\sigma_{j}\leq\sigma^{*}$, and $E|Z_{ij}|^{\theta}\leq
K_{\theta}^{\theta}$.

######  Theorem 3.1.

Assume Condition 3.1 and $\mu_{i}=0,1\leq i\leq p$. (i) Let $\theta>2$. Then
there exists a Gaussian process $(G^{\star}_{j})_{0\leq j\leq p-k}$ such that
on a possibly larger probability space, for all $u>0$,

$\displaystyle P\left[(nk)^{1/2}\Delta^{\star}\geq
c_{0}u\right]\leq{{npK_{\theta}^{\theta}}\over
u^{\theta}\sigma_{*}^{\theta}},$ (3.6)

where $c_{0}$ is a constant only depending on $\theta$. (ii) Let $\theta>3$.
Then the distributional distance

${\rho^{*}}\lesssim(nk)^{-1/6}(\log
p)^{7/6}+(np/(nk)^{\theta/2})^{1/(\theta+1)}(\log
p)^{(3\theta-2)/(2+2\theta)},$ (3.7)

where the constant in $\lesssim$ only depends on $\theta$, $\sigma_{*}$,
$\sigma^{*}$ and $K_{\theta}$.

We emphasize that our theorem does not require $n\to\infty$ and it is also
applicable when $n$ is finite. For example, when $n=2$ observations are
available, if we choose the window size $k$ be sufficiently large such that
$p(\log p)^{3\theta/2-1}=o(k^{\theta/2})$, then by (3.7) and elementary
manipulations we still have ${\rho^{*}}\to 0$. Under the slightly weaker
condition $p^{2/\theta}=o(k)$, $R_{j}^{\star}$ and $G_{j}^{\star}$ are
uniformly close to each other in the sense of $\max_{k\leq j\leq
p-k}|R_{j}^{\star}-G_{j}^{\star}|=o_{P}(1)$ in view of (3.6).

###  3.2 Calculating cutoff values

If the variances $\sigma_{j}^{2}$ are known, given the level $0<\alpha<1$, we
can choose the cutoff value $u=u_{1-\alpha}$ such that

$\displaystyle P(\max_{0\leq j\leq p-k}G_{j}^{\star}\geq
u_{1-\alpha})=\alpha.$ (3.8)

The above can be done by Monte Carlo simulations. Assuming that $p,n,k$
satisfy the relation ${(np)^{2/\theta}(\log p)^{3-2/\theta}=o(nk).}$ Then the
right hand side of (3.7) goes to $0$. By Theorem 3.1, the test $\max_{0\leq
j\leq p-k}R_{j}>u_{1-\alpha}$ has the asymptotically correct size $\alpha$.

####  3.2.1 Estimation of block sum variances

In general, however, the variances $\sigma_{j}^{2}$ are not known. Since we
have multiple realizations, we can estimate them by the classical unbiased
variance estimate

$\displaystyle\hat{\sigma}_{j}^{2}={1\over{n-1}}\sum_{i=1}^{n}(Y_{ij}-\hat{\mu}_{j})^{2}.$
(3.9)

Correspondingly, our test statistic $R_{j}^{\star}$ in (3.2) now becomes

$\displaystyle\hat{R}_{j}={{\sum_{l=j+1}^{j+k}\sqrt{n}\hat{\mu}_{l}}\over{(\sum_{l=j+1}^{j+k}\hat{\sigma}^{2}_{l})^{1/2}}},\,\,0\leq
j\leq p-k.$ (3.10)

At first glance, if $n$ is small, $\hat{\sigma}_{j}^{2}$ may deviate
substantially from $\sigma_{j}^{2}$. For example, if $n=2$, then
$\hat{\sigma}_{j}^{2}=(Y_{1j}-Y_{2j})^{2}/2,$ which may be quite different
from $\sigma_{j}^{2}$. This difference might suggest that replacing
$\sigma_{l}^{2}$ in $R_{j}$ by $\hat{\sigma}_{l}^{2}$ can be problematic.
However, interestingly, under suitable conditions on $p,k,n$, $R_{j}^{\star}$
and $\hat{R}_{j}$ can still be uniformly close. This can be intuitively
explained by the fact that, in $R_{j}$, it is the block sum variance
$v_{j}=\sum_{l=j+1}^{j+k}\sigma^{2}_{l}$ that is directly involved, not just a
single $\sigma_{l}^{2}$. The sum
$\hat{v}_{j}=\sum_{l=j+1}^{j+k}\hat{\sigma}^{2}_{l}$ can still be a good
estimate of $v_{j}$, despite that individually the difference
$\hat{\sigma}_{j}^{2}-\sigma_{j}^{2}$ can be big due to a small $n$. The
convergence rate is given in the following Proposition 3.1. It implies that,
under Condition 3.1, if $p=o(n^{\theta/2-1}k^{\theta/2})$, then
$\hat{v}_{j}/v_{j}$ is uniformly close to $1$. It is proved in the
Supplementary Material.

######  Proposition 3.1.

Let Condition 3.1 be satisfied. If $\theta>4$, we have

$P(n\max_{0\leq j\leq
p-k}|\hat{v}_{j}-v_{j}|>u)\lesssim{{npK_{\theta}^{\theta}}\over{u^{\theta/2}}}+{p\over
k}\exp(-c_{3}{u^{2}\over{nkK_{4}^{4}}}),$ (3.11)

where the constant in $\lesssim$ and $c_{3}>0$ only depend on $\theta$. If
$2<\theta\leq 4$, then

$\displaystyle P(n\max_{0\leq j\leq
p-k}|\hat{v}_{j}-v_{j}|>u)\lesssim{{npK_{\theta}^{\theta}}\over{u^{\theta/2}}}.$
(3.12)

Note that (3.11) of Proposition 3.1 implies that we have the uniform
convergence rate

$n\max_{0\leq j\leq
p-k}|\hat{v}_{j}-v_{j}|=O_{P}((np)^{2/\theta}+(nk)^{1/2}\log p).$

Under Condition 3.1, $k\sigma_{*}^{2}\leq v_{j}\leq k(\sigma^{*})^{2}$. Thus
the term $n\max_{0\leq j\leq p-k}|\hat{v}_{j}-v_{j}|$ in (3.11) can be
replaced by the ratio normalized version $nk\max_{0\leq j\leq
p-k}|\hat{v}_{j}/v_{j}-1|$ so that (3.11) is still valid with the constants in
$\lesssim$ and $c_{3}$ therein depending on $\theta$, $\sigma_{*}$ and
$\sigma^{*}$. By elementary calculations, if
$p=o(n^{\theta/2-1}k^{\theta/2})$, the ratios $\hat{v}_{j}/v_{j}$ are
uniformly close to $1$ in the sense that $\max_{0\leq j\leq
p-k}|\hat{v}_{j}/v_{j}-1|=o_{P}(1)$.

####  3.2.2 A bootstrap calibration procedure

To perform the test for $H_{0}:\mu_{j}\equiv 0$ based on $\hat{R}_{j}$ with
$\sigma_{j}^{2}$ replaced by their estimates $\hat{\sigma}_{j}^{2}$, we need
to estimate the corresponding cutoff value $u_{1-\alpha}$ based on (3.8).
Recall that $\Gamma(\sigma)=(\gamma_{j,j^{\prime}}(\sigma))_{k\leq
j,j^{\prime}\leq p-k}$ is the covariance matrix for the vector $(Z_{j})_{k\leq
j\leq p-k}$. Write $u_{1-\alpha}=q_{\alpha}(\sigma)$ as a function of
$\sigma=(\sigma_{1},\ldots,\sigma_{p})$. Write
$\hat{\sigma}=(\hat{\sigma}_{1},\ldots,\hat{\sigma}_{p})$ and
$\hat{u}_{1-\alpha}=q_{\alpha}(\hat{\sigma})$ which satisfies

$\displaystyle P^{*}(\max_{0\leq j\leq
p-k}G^{*}_{j}\geq\hat{u}_{1-\alpha})=\alpha,$ (3.13)

where $P^{*}$ is the probability measure given $Y=(Y_{1},\ldots,Y_{n})$ and,
given $\hat{\sigma}$, $(G^{*}_{j})_{k\leq j\leq p-k}$ is mean $0$ Gaussian
vector with covariance matrix $\Gamma(\hat{\sigma})$. In particular, as (3.3),
we can define

$\displaystyle G^{*}_{j}={W_{j}^{*}\over{\hat{v}_{j}^{1/2}}},\mbox{ where
}W_{j}^{*}={\sum_{l=j+1}^{j+k}\hat{\sigma}_{l}\eta_{l}}$ (3.14)

and $\eta_{l},l\in\mathbb{Z}$, are i.i.d. $N(0,1)$ random variables that are
independent of $Y=(Y_{1},\ldots,Y_{n})$. Given $\hat{\sigma}$, the cutoff
value $\hat{u}_{1-\alpha}$ in (3.13) can also be computed by extensive
simulations.

The following theorem shows the validity of the above plug-in method in the
sense that the size of our test is close to $\alpha$. It is proved in the
Supplementary Material.

######  Theorem 3.2.

Let $t_{1}=(\log p)^{1/2}(nk)^{-1/2}$, $t_{2}=(pn(nk)^{-\theta/2}(\log
p)^{-2/3})^{1/(1/3+\theta/2)}$ and $t_{*}=\max(t_{1},t_{2})$. Recall (3.7) for
$\rho^{*}$. Let $0<\alpha<1$. Then under Condition 3.1, we have

$\displaystyle\quad|P(\max_{0\leq j\leq
p-k}\hat{R}_{j}\geq\hat{u}_{1-\alpha})-\alpha|\lesssim\rho^{*}+t_{*}^{1/3}(\log
p)^{2/3},$ (3.15)

where the constant in $\lesssim$ only depends on
$\sigma_{*},\sigma^{*},\theta$ and $K_{\theta}$. In particular, the right hand
side of (3.15) is $o(1)$ if $pn(\log p)^{3\theta/2-1}=o((nk)^{\theta/2}).$

###  3.3 Estimating break-points based on one-sided test

Algorithm 3.1 shows estimating break-points based on the one-sided test. It
uses $R_{j}^{\star}$ assuming that $\sigma_{j}$, $1\leq j\leq p$, are known.
If not known, we shall use the estimates $\sigma_{j}^{2}$ in (3.9). Same
simulation assisted $\gamma$ and $\delta$ can be used as in the one
realization one-sided test case. Theorem 3.3 provides theoretical properties
of the break-point estimates.

######  Algorithm 3.1.

Step 1. Let $L_{i}^{\star}=R_{i-k}^{\star},i=k,\ldots,p-k$ and denote
$Q_{i}^{\star}=1(R_{i}^{\star}>\gamma)+1(L_{i}^{\star}>\gamma)$ for a pre-
specified cutoff value $\gamma$. We use a majority vote approach to smooth
$Q_{i}^{\star}.$ Specifically, denote
$j_{0}^{\star}=\sum_{i=j-k}^{j+k}I\\{Q_{i}^{\star}=0\\},j_{1}^{\star}=\sum_{i=j-k}^{j+k}I\\{Q_{i}^{\star}=1\\},$
and $j_{2}^{\star}=\sum_{i=j-k}^{j+k}I\\{Q_{i}^{\star}=2\\}.$ Let
$\tilde{Q}_{j}^{\star}=\\{k,\mbox{such that}\
j_{k}^{\star}=\mbox{max}_{l\in\\{0,1,2\\}}j^{\star}_{l}\\}.$

Step 2. Decompose $\\{1,\ldots,p\\}=W_{0}\cup W_{1}\cup W_{2},$ where $i\in
W_{0}$ if $\tilde{Q}_{i}^{\star}=0,i\in W_{1}$ if $\tilde{Q}_{i}^{\star}=1$
and $i\in W_{2}$ if $\tilde{Q}_{i}^{\star}=2.$ Let ${\cal M}_{1},\ldots,{\cal
M}_{\hat{l}}$ be connected components of $W_{1}.$

Step 3. Let
$R_{j}^{\flat}={{\sum_{f=j+1}^{j+k}\sqrt{n}\hat{\mu}_{f}}/\sqrt{k}}$ and
$L_{j}^{\flat}=R_{j-k}^{\flat}$. Given $\delta<\gamma$, the break-points are
defined as $\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{R_{j}^{\flat}:\,L_{j}^{\star}\leq\delta\\}$ if ${\cal M}_{i}$ is the
transition region from $W_{0}$ to $W_{2}$. If ${\cal M}_{i}$ is the transition
region from $W_{2}$ to $W_{0},$ $\hat{\tau}_{i}=\mbox{argmax}_{j\in{\cal
M}_{i}}\\{L_{j}^{\flat}:\,R_{j}^{\star}\leq\delta\\}$.

Differently from Algorithm 2.1, in Step 3 of Algorithm 3.1 we use
$R_{j}^{\flat}$ instead of $R_{j}^{\star}$ in the argmax function. The reason
is for technical convenience: one has monotonicity
$E(R_{j}^{\flat})<E(R_{i}^{\flat})$ for $\tau_{1}-k<j<i\leq\tau_{1}$, which
tends to make the estimated break-point closer to $\tau_{1}$. In comparison
$E(R_{j}^{\star})$ is generally not monotone, since the variances
$\sigma^{2}_{j}$ can be unequal.

######  Theorem 3.3.

Assume Conditions 2.2, 2.5, $Z_{ij}$ are $\sigma^{2}$-sub-Gaussian, and
$2\sigma\gamma\leq d\sqrt{nk}$. Let
$m=\lfloor 2k^{1/2}\delta\sigma n^{-1/2}d^{-1}\rfloor$. Then

$\displaystyle 1-P\left[\hat{l}=l,\,\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq m\right]\lesssim{p\over
k}\exp(-c_{1}\gamma^{2})+l\exp(-c_{2}\delta^{2}),$ (3.16)

where the constant in $\lesssim$ and $c_{1},c_{2}>0$ are independent of
$k,d,n$ and $p$.

Theorem 3.3 is proved in the Supplementary Material. Assume that $(\log
p)(\log l)=o(n^{2}d^{4})$ and $k$ satisfies $(d^{2}n)^{-1}\log p=o(k)$ and
$k=o(nd^{2}/\log l)$. Let $\gamma=C_{1}(\log p)^{1/2}$, $\delta=C_{2}(\log
l)^{1/2}$, where $C_{1},C_{2}>0$ are constants. Then the right hand side of
(3.16) can be arbitrarily small by letting $C_{1},C_{2}$ sufficiently large.
Theorem 3.3 implies that we can have exact recovery with probability
$P[\hat{l}=l,\,\max_{j\leq l}|\hat{\tau}_{j}-\tau_{j}|=0]\to 1$ since
$k^{1/2}\delta n^{-1/2}d^{-1}\to 0$.

###  3.4 Two-sided test: A U-statistic approach

In the one-realization case, we use (2.15) to test the two-sided alternative
$H_{1}$ of (2.3). If multiple realizations $Y_{i}=(Y_{i1},\ldots,Y_{ip})^{T}$,
$1\leq i\leq n$, are available, we can use the $U$-statistic

$\displaystyle W_{j}={2\over{n(n-1)}}\sum_{1\leq i<i^{\prime}\leq
n}Y_{ij}Y_{i^{\prime}j},$ (3.17)

which is an unbiased estimate of $\mu_{j}^{2}$. This is different from (2.15)
in that $X_{j}^{2}$ in the latter is not an unbiased estimate of
$\mu_{j}^{2}$. As an important consequence, we remark that unlike the two-
sided test in Section 2.2, here we do not need to use $\kappa$ of form (2.17).
Under $H_{0}$, the variance of $W_{j}$ is $2\sigma_{j}^{4}/(n(n-1))$. Define

$\displaystyle
R_{j,4}=((n^{2}-n)/2)^{1/2}{{W_{j+1}+\ldots+W_{j+k}}\over{(\sigma^{4}_{j+1}+\ldots+\sigma^{4}_{j+k})^{1/2}}}.$
(3.18)

Let $\eta_{i},i\in\mathbb{Z},$ be i.i.d. $N(0,1)$. Define the Gaussian process

$\displaystyle
G_{j,4}={{\sigma_{j+1}^{2}\eta_{j+1}+\ldots+\sigma_{j+k}^{2}\eta_{j+k}}\over{(\sigma^{4}_{j+1}+\ldots+\sigma^{4}_{j+k})^{1/2}}}.$
(3.19)

######  Theorem 3.4.

Assume Condition 3.1 and $\mu_{i}=0,1\leq i\leq p$. Then the distributional
distance

$\displaystyle\sup_{u}|P(\max_{0\leq j\leq p-k}R_{j,4}\geq u)-P(\max_{0\leq
j\leq p-k}G_{j,4}\geq u)|$ (3.20) $\displaystyle\lesssim k^{-1/6}(\log
p)^{7/6}+(pk^{-\theta/2})^{1/(\theta+1)}(\log p)^{(3\theta-2)/(2+2\theta)}.$
(3.21)

In $R_{j,4}$, the quantity $\sigma_{j}^{4}$ is typically unknown. Here we
shall propose an unbiased estimate. Note that the natural estimate
$(\hat{\sigma}_{j}^{2})^{2}$ with $\hat{\sigma}_{j}^{2}$ given in (3.9) is not
unbiased. Let

$\displaystyle\quad\hat{\omega}_{j}={{1/4}\over{n(n-1)(n-2)(n-3)}}\sum(Y_{ij}-Y_{i^{\prime}j})^{2}(Y_{hj}-Y_{h^{\prime}j})^{2},$
(3.22)

where the sum is over mutually different indexes
$i,i^{\prime},h,h^{\prime}\in\\{1,\ldots,n\\}$. Clearly
$E(\hat{\omega}_{j})=\sigma_{j}^{4}$. Similar to (3.10), consider the realized
version

$\displaystyle
R^{*}_{j,4}=((n^{2}-n)/2)^{1/2}{{W_{j+1}+\ldots+W_{j+k}}\over{(\hat{\omega}_{j+1}+\ldots+\hat{\omega}_{j+k})^{1/2}}}.$
(3.23)

To test $H_{0}$ vs $H_{1}$ in (2.3), we reject $H_{0}$ at level
$\alpha\in(0,1)$ if $\max_{0\leq j\leq p-k}R^{*}_{j,4}\geq q_{1-\alpha}$ for
some cutoff value $q_{1-\alpha}$. As in (3.13), $q_{1-\alpha}$ can be
approximated by $\hat{q}_{1-\alpha}$, which satisfies $P^{*}(\max_{0\leq j\leq
p-k}G^{*}_{j,4}\geq\hat{q}_{1-\alpha})=\alpha$, where

$\displaystyle
G^{*}_{j,4}={{\hat{\omega}_{j+1}^{1/2}\eta_{j+1}+\ldots+\hat{\omega}_{j+k}^{1/2}\eta_{j+k}}\over{(\hat{\omega}_{j+1}+\ldots+\hat{\omega}_{j+k})^{1/2}}}$
(3.24)

a Gaussian process conditioning on $(Y_{1},\ldots,Y_{n})$. As a slightly
modified version, noting that for i.i.d. $N(0,1)$ random variables
$Z_{1},\ldots Z_{n}$, the $U$-statistic $2\sum_{1\leq i<i^{\prime}\leq
n}Z_{i}Z_{i^{\prime}}=(n^{2}-n)\bar{Z}_{n}^{2}-\sum_{i=1}^{n}(Z_{i}-\bar{Z}_{n})^{2}$
is identically distributed as $\zeta=(n-1)\chi^{2}_{1}-\chi^{2}_{n-1}$, where
the $\chi^{2}$ random variables $\chi^{2}_{1}$ and $\chi^{2}_{n-1}$ are
independent, we can use

$\displaystyle
G^{\diamond}_{j,4}={{\hat{\omega}_{j+1}^{1/2}\zeta_{j+1}+\ldots+\hat{\omega}_{j+k}^{1/2}\zeta_{j+k}}\over{(\hat{\omega}_{j+1}+\ldots+\hat{\omega}_{j+k})^{1/2}}},$
(3.25)

where $\zeta_{i}$ are independent and identically distributed as
$\zeta/\sqrt{2n(n-1)}$. If $n$ is big, $G^{\diamond}_{j,4}$ gives a better
approximation.

In the definition of $\hat{\omega}_{j}$ in (3.22), it involves a 4-fold
summation with $O(n^{4})$ computation complexity. Interestingly, we can have
the following expression which allows computing $\hat{\omega}_{j}$ within only
$O(n)$ steps: elementary but tedious calculations show

$\hat{\omega}_{j}={{(n-1)(4S_{j,3}S_{j,1}-nS_{j,4}-3S_{j,2}^{2})+(nS_{j,2}-S_{j,1}^{2})^{2}}\over{n(n-1)(n-2)(n-3)}},\mbox{
where }S_{j,l}=\sum_{i=1}^{n}Y_{ij}^{l}.$ (3.26)

To compute $W_{j}$ in (3.17), we use the well-known formula
$W_{j}={{S_{j,1}^{2}-S_{j,2}}/(n(n-1))}.$

####  3.4.1 Estimating break-points based on two-sided test

Similar to Algorithm 2.2, we can adjust Algorithm 3.1 for locating break-
points based on two-sided test. Theorem 3.5 shows theoretical properties of
Algorithm 3.2 and it is proved in the Supplementary Material.

######  Algorithm 3.2.

Step 1. Let $Q_{j,4}=1(R_{j,4}>\gamma)+1(L_{j,4}>\gamma)$ for a pre-specified
cutoff value $\gamma$, $j=k,\ldots,p-k.$ We use a majority vote approach to
smooth $Q_{j,4}.$ Specifically, denote
$l_{0}=\sum_{i=j-k}^{j+k}I\\{Q_{j,4}=0\\},l_{1}=\sum_{i=j-k}^{j+k}I\\{Q_{j,4}=1\\},$
and $l_{2}=\sum_{i=j-k}^{j+k}I\\{Q_{j,4}=2\\}.$ Let
$\tilde{Q}_{j,4}=\\{k,\mbox{such that}\
l_{k}=\mbox{max}_{j\in\\{0,1,2\\}}l_{j}\\}.$

Step 2. Decompose $\\{1,\ldots,p\\}=W_{0}\cup W_{1}\cup W_{2},$ where $j\in
W_{0}$ if $\tilde{Q}_{j,4}=0,j\in W_{1}$ if $\tilde{Q}_{j,4}=1$ and $j\in
W_{2}$ if $\tilde{Q}_{j,4}=2.$ Let ${\cal M}_{1},\ldots,{\cal M}_{\hat{l}}$ be
connected components of $W_{1}.$

Step 3. Let
$R^{\ddagger}_{j}=((n^{2}-n)/2)^{1/2}\sum_{h=1}^{k}W_{j+h}/{k^{1/2}}$ and
$L^{\ddagger}_{j}=R^{\ddagger}_{j-k}.$ Given $\delta$, the break-points are
defined as $\hat{\tau}_{f}=\mbox{argmax}_{j\in{\cal
M}_{f}}\\{R_{j}^{\ddagger}:\,L_{j,4}\leq\delta\\}$ if ${\cal M}_{f}$ is the
transition region from $W_{0}$ to $W_{2}$. If ${\cal M}_{f}$ is the transition
region from $W_{2}$ to $W_{0},$ $\hat{\tau}_{f}=\mbox{argmax}_{j\in{\cal
M}_{f}}\\{L_{j}^{\ddagger}:\,R_{j,4}\leq\delta\\}$.

######  Theorem 3.5.

Assume Condition 2.4 and $Z_{ij}$ are $\sigma^{2}$-sub-Gaussian. Let
$\gamma=c_{1}(\log p)^{1/2}$, $\delta=c_{2}(\log l)^{1/2}$, where
$c_{1},c_{2}>0$ are sufficiently large constants. Assume that $(\log
p)^{1/2}=o(d^{2}n\sqrt{k})$. Then there exists a constant $c>0$ independent of
$n$, $k$ and $p$ such that

$\displaystyle P\left[\hat{l}=l,\,\max_{j\leq
l}|\hat{\tau}_{j}-\tau_{j}|\leq{ck^{1/2}(\log l)^{1/2}\over nd^{2}}\right]\to
1.$ (3.27)

##  4 Numerical studies

In this section, we present simulation studies, assess the finite sample
performance of the proposed methods and compare them with competing methods
[12, 39, 2]. We look at one- and two-sided tests with one realization in
Section 4.1 and Section 4.2, respectively. One- and two-sided tests with more
realizations are presented in the Supplementary Material.

###  4.1 Simulation study 1

Consider the model $X_{i}=\mu_{i}+Z_{i},1\leq i\leq p.$ The number of tests
are $p=600,2000$ and $6000.$ There are $2$ break-points
$\tau_{1}=1+0.4p,\tau_{2}=0.6p,$ $1$ signal cluster $[\tau_{1},\tau_{2}]$ and
the configuration is displayed in Table 1 and Figure 3. We compare it with the
change point detection for epidemic alternative proposed in [39]. We simulate
data with three different error terms: standard normal distribution, rescaled
student $t$ distribution with $6$ degree of freedom ($t(6)/1.5^{0.5}$) and
rescaled Laplace distribution ($\mbox{LP}(0,1)/2^{0.5}$) so that their
variances are all $1$.

The sliding window length $k=\lfloor p^{1/2}\rfloor$ is used in the
calculation of $R_{i}^{\circ}=k^{-1}\sum_{j=i+1}^{i+k}X_{j}.$ We also show
results for other choices of $k.$ In order to estimate the variance
$\sigma^{2}$, we choose the tuning parameter $m=k.$ Let $p^{\prime}=p-m+1$,
$\hat{\sigma}_{i}^{2}=m^{-1}\sum_{j=i}^{i+m-1}X_{i}^{2},1\leq i\leq
p^{\prime}.$ Theoretically speaking, any statistics $\hat{\sigma}^{2}_{(j)}$
with $j\leq p^{\prime}/2$ are consistent and we use
$\hat{\sigma}^{2}_{(\lfloor{p^{\prime}/2\rfloor})}$ as the estimate.

Table 1: Signal configuration for the one-sided test. seq 1: the linear sequence from 0.4 to 1.6; seq 2: the linear sequence from 1.6 to 0.4. Segment means percentage of the sequence. Segment | 40 | 10 | 10 | 40
---|---|---|---|---
Signal | 0 | seq 1 | seq 2 | 0
Figure 1: Signal configuration when $p=600.$

We implemented algorithm 2.1. Thresholding values $\gamma$ and $\delta$ are
chosen as $0.95$th quantile of $\hat{\sigma}\mbox{max}_{0\leq j\leq
p-k}G_{j}^{\circ}$ and and $\hat{\sigma}\mbox{max}_{j\in W_{1}}G_{j}^{\circ},$
respectively, where
$G_{j}^{\circ}=\sum_{i=j+1}^{j+k}\eta_{i}/k,\eta_{i},i\in\mathbb{Z},$ are
i.i.d. $N(0,1)$ and $W_{1}$ are the major connected components which include
indices $j$ such that
$Q_{j}^{\circ}=1(R_{j}^{\circ}>\gamma)+1(L_{j}^{\circ}>\gamma)=1.$

In implementing [39], we use $L_{1}$, the likelihood ratio statistic as an
example for illustration. Similar results can be obtained for other test
statistics. Specifically,

${L_{1}}=\mbox{max}_{1\leq i<j\leq
p}\\{\sum_{k=i+1}^{j}X_{k}-\frac{j-i}{p}\sum_{k=1}^{p}X_{k}-\frac{1}{2}\delta_{0}(j-i)\\},$
(4.1)

where $\delta_{0}$ is the signal magnitude, which is assumed to be the same
within a cluster in [39]. In our setup, we take $\delta_{0}=1,$ which is the
average of signal magnitude within the cluster $[\tau_{1},\tau_{2}].$ We
identify the region $[\hat{I},\hat{J}]$ as the epidemic alternative, where
$\sum_{k=\hat{I}+1}^{\hat{J}}X_{k}-p^{-1}(\hat{J}-\hat{I})\sum_{k=1}^{p}X_{k}-\frac{1}{2}\delta_{0}(\hat{J}-\hat{I})$
is the obtained maximum value in (4.1). Note that the computational speed is
quadratic with number of tests $p.$ Our evaluation criterion is the combined
error rate (CER), which is the expected value of the ratio of the number of
falsely rejected hypotheses and falsely accepted hypotheses over total number
of tests, the estimated number of break points $\hat{l}$ and the average
difference between the estimated break points and true break points. For the
proposed method, we also look at false discovery rate (FDR), which is the
expected value of the ratio of false rejections over total rejections and the
power, which is the expected value of the the ratio of true rejections over
total number of non-nulls.

Table 2 summarizes results based on $10^{3}$ replications. We can see that
across different error distributions, the variance estimate $\hat{\sigma}^{2}$
has a decent performance and, as expected from our asymptotic theory, it is
close to the true ones. The proposed method has smaller CER compared to method
based on [39], especially with large number of tests. Both procedures
correctly identified 2 break points. The difference between estimated break
points and true ones are smaller based on the proposed method especially with
large samples. Our results are robust to different error terms and the sliding
window length $k.$ For different error distributions the respective values of
CER are quite close, as expected from our theoretical result.

Table 2: Summary statistics for one-sided test with $1,000$ simulations. $N(0,1)$: standard normal; $t(6)/1.5^{0.5}$: rescaled student $t$ distribution with df $6$; $\mbox{LP}(0,1)/2^{0.5}$: rescaled Laplace distribution; $k$ is the window size; CER is computed based on the proposed method; $\mbox{CER}_{Y}$ is based on [39]; $\hat{l}$ is estimated number of break points based on the proposed method; $\hat{l}_{Y}$ is estimated number of change points based on [39]; Diff is the average distance between estimated break points and true break points based on the proposed method; $\mbox{Diff}_{Y}$ is the average distance between estimated change points and true change points based on [39]; FDR is the expected value of the ratio of false rejections over total rejections and Power is the expected value of the the ratio of true rejections over total number of non-nulls.  $k$ | $\hat{\sigma}^{2}$ | CER | $\mbox{CER}_{Y}$ | $\hat{l}$ | $\hat{l}_{Y}$ | Diff | $\mbox{Diff}_{Y}$ | FDR | Power
---|---|---|---|---|---|---|---|---|---
$p=600$ |  |  |  |  |  |  |  |  |
|  |  |  | $N(0,1)$ |  |  |  |
24 | 1.0533 | 0.0503 | 0.0538 | 2 | 2 | 15.35 | 16.63 | 0.0016 | 0.75
30 | 1.0665 | 0.0475 | 0.0508 | 2 | 2 | 14.30 | 15.68 | 0.0029 | 0.77
36 | 1.0605 | 0.0528 | 0.0492 | 2 | 2 | 15.85 | 15.24 | 0.0021 | 0.74
|  |  |  | $t(6)/1.5^{0.5}$ |  |  |  |
24 | 1.0363 | 0.0489 | 0.0513 | 2 | 2 | 14.66 | 15.84 | 0.0015 | 0.76
30 | 1.0312 | 0.0511 | 0.0533 | 2 | 2 | 15.33 | 16.44 | 0.0019 | 0.75
36 | 1.0425 | 0.0554 | 0.0543 | 2 | 2 | 16.61 | 16.75 | 0.0020 | 0.73
|  |  |  | $\mbox{LP}(0,1)/2^{0.5}$ |  |  |  |
24 | 1.0128 | 0.0517 | 0.0528 | 2 | 2 | 18.50 | 16.33 | 0.0033 | 0.74
30 | 1.0377 | 0.0532 | 0.0548 | 2 | 2 | 17.08 | 16.85 | 0.0051 | 0.74
36 | 1.0630 | 0.0528 | 0.0497 | 2 | 2 | 15.84 | 15.36 | 0.0010 | 0.74
$p=2,000$ |  |  |  |  |  |  |  |  |
|  |  |  | $N(0,1)$ |  |  |  |
44 | 1.0469 | 0.0262 | 0.0495 | 2 | 2 | 29.54 | 50.05 | 0.0021 | 0.87
55 | 1.0535 | 0.0244 | 0.0514 | 2 | 2 | 25.96 | 51.90 | 0.0016 | 0.88
66 | 1.0420 | 0.0251 | 0.0499 | 2 | 2 | 25.15 | 50.36 | 0.0025 | 0.88
|  |  |  | $t(6)/1.5^{0.5}$ |  |  |  |
44 | 1.0342 | 0.0279 | 0.0505 | 2 | 2 | 31.46 | 50.97 | 0.0015 | 0.86
55 | 1.0394 | 0.0248 | 0.0503 | 2 | 2 | 26.64 | 50.82 | 0.0019 | 0.88
66 | 1.0355 | 0.0265 | 0.0518 | 2 | 2 | 26.52 | 52.28 | 0.0032 | 0.87
|  |  |  | $\mbox{LP}(0,1)/2^{0.5}$ |  |  |  |
44 | 1.0382 | 0.0278 | 0.0524 | 2 | 2 | 32.84 | 52.91 | 0.0006 | 0.86
55 | 1.0569 | 0.0228 | 0.0500 | 2 | 2 | 22.77 | 50.46 | 0.0018 | 0.89
66 | 1.0475 | 0.0111 | 0.0498 | 2 | 2 | 26.15 | 50.33 | 0.0017 | 0.87
$p=6,000$ |  |  |  |  |  |  |  |  |
|  |  |  | $N(0,1)$ |  |  |  |
60 | 1.0433 | 0.0170 | 0.0495 | 2 | 2 | 51.51 | 148.85 | 0.0007 | 0.92
77 | 1.0396 | 0.0116 | 0.0489 | 2 | 2 | 40.67 | 147.34 | 0.0009 | 0.94
100 | 1.0395 | 0.0108 | 0.0509 | 2 | 2 | 32.34 | 153.22 | 0.0015 | 0.95
|  |  |  | $t(6)/1.5^{0.5}$ |  |  |  |
60 | 1.0308 | 0.0168 | 0.0505 | 2 | 2 | 66.74 | 152.02 | 0.0008 | 0.92
77 | 1.0330 | 0.0127 | 0.0505 | 2 | 2 | 43.87 | 151.99 | 0.0014 | 0.94
100 | 1.0345 | 0.0103 | 0.0507 | 2 | 2 | 30.94 | 152.71 | 0.0011 | 0.95
|  |  |  | $\mbox{LP}(0,1)/2^{0.5}$ |  |  |  |
60 | 1.0429 | 0.0181 | 0.0493 | 2 | 2 | 65.60 | 148.27 | 0.0004 | 0.91
77 | 1.0465 | 0.0136 | 0.0499 | 2 | 2 | 46.42 | 150.07 | 0.0009 | 0.93
100 | 1.0452 | 0.0111 | 0.0498 | 2 | 2 | 33.16 | 150.05 | 0.0014 | 0.95

Per the request of a referee, we implement the BH procedure [2] with Gaussian
error term and summarize the results in Table 3. The simulation set up is the
same as that in Table 2. At FDR level $5\%,$ we can see that the BH procedure
always controls FDR but with low power for clustered signals.

Table 3: BH procedure with Gaussian error term. The definition of FDR, Power and CER are the same as that in Table 2. $p$ | $\mbox{FDR}_{\mbox{BH}}$ | $\mbox{Power}_{\mbox{BH}}$ | $\mbox{CER}_{\mbox{BH}}$
---|---|---|---
600 | 0.0435 | 0.0011 | 0.0724
6000 | 0.0433 | 0.0003 | 0.1101

We also conduct simulation studies to check the empirical type-I error rates
under the global null with Gaussian error term. The results are summarized in
Table 4. At significance level $5\%,$ the proposed method has a similar type-I
error rate to BH procedure under the global null as evidenced from Table 4.

Table 4: Type-I error rates under the global null with Gaussian error term. FDR represents FDR based on the proposed method and $\mbox{FDR}_{\mbox{BH}}$ represents FDR based on BH procedure.  $p$ | $k$ | FDR | $\mbox{FDR}_{\mbox{BH}}$
---|---|---|---
600 | 36 | 0.0594 | 0.0495
6000 | 60 | 0.0396 | 0.0495

###  4.2 Simulation study 2

In this section, we examine the two-sided test procedure. Data is generated
through model (2.1). Let $p=600,2000$ and $6000$. The signal configuration is
summarized in Table 5. We look at the robustness of our procedure with
different error terms ($N(0,1)/2^{0.25}$, $t(10)/(75/16)^{0.25}$ and
$LP(0,1)/20^{0.25}$), which are standardized to have $\kappa=1$. Window size
$k=\lfloor p^{1/2}\rfloor$ and $m=\lfloor p^{1/2}\rfloor$ are used for
illustration. The calculation of $\hat{\sigma}^{2},\gamma$ and $\delta$ are
the same as that in simulation study 1 except that $\hat{\kappa}$ is used
instead of $\hat{\sigma}$ and the calculation of $\hat{\kappa}$ is through
$\hat{\kappa}^{2}=\hat{\nu}_{(k)}/2-4\hat{\sigma}^{4}_{(k)}$. We follow
Algorithm 2.2 to implement our method. As a comparison, results based on true
values of $\sigma^{2}$ and $\kappa$ are presented as well.

Table 5: Signal configuration for the two-sided test. "-1 and 1 alternating": $\mu_{i}$ is $-1$ if $i$ is odd and $1$ if $i$ is even, $\mbox{seq}(0.5,1.5)$: a linear sequence from $0.5$ to $1.5$ and $\mbox{seq}(1.5,0.5)$: a linear sequence from $1.5$ to $0.5.$ Segment ($\%$ ) | 30 | 10 | 20 | 5 | 5 | 30
---|---|---|---|---|---|---
Signal strength | 0 |  $-1$ and $1$ alternating | 0 | $\mbox{seq}(0.5,1.5)$ | $\mbox{seq}(1.5,0.5)$ | 0

From Table 6 we can see that procedures using the estimated parameters and the
true ones have a comparable performance in terms of CER, FDR, power, estimated
number of break points and the difference between estimated break points and
true break points. This is consistent with our large sample theory. The
results are relatively robust across different error terms. As numbers of
tests increase, CER and FDR decrease and power and the difference between
estimated break points and true break points increase.

Table 6: Summary statistics for two-sided test with $1,000$ simulations. Underscore $e$ is based on estimated $\sigma^{2}$ and $\kappa$, and underscore $t$ is based on true $\sigma^{2}$ and $\kappa$. p | $\mbox{CER}_{e}$ | $\mbox{CER}_{t}$ | $\mbox{FDR}_{e}$ | $\mbox{FDR}_{t}$ | $\mbox{Power}_{e}$ | $\mbox{Power}_{t}$ | $\hat{l}_{e}$ | $\hat{l}_{t}$ | $\mbox{Diff}_{e}$ | $\mbox{Diff}_{t}$
---|---|---|---|---|---|---|---|---|---|---
$N(0,1)/2^{0.25}$ |  |  |  |  |  |  |  |  |  |
600 | 0.0822 | 0.0743 | 0.0207 | 0.0323 | 0.61 | 0.65 | 4 | 4 | 19.67 | 21.53
2000 | 0.0390 | 0.0354 | 0.0091 | 0.0142 | 0.81 | 0.84 | 4 | 4 | 19.47 | 19.81
6000 | 0.0223 | 0.0202 | 0.0048 | 0.0101 | 0.89 | 0.91 | 4 | 4 | 35.64 | 38.68
$t(10)/(75/16)^{0.25}$ |  |  |  |  |  |  |  |  |  |
600 | 0.0817 | 0.0753 | 0.0266 | 0.0367 | 0.61 | 0.65 | 4 | 4 | 17.62 | 18.84
2000 | 0.0378 | 0.0339 | 0.0129 | 0.0203 | 0.82 | 0.85 | 4 | 4 | 36.37 | 30.40
6000 | 0.0208 | 0.0192 | 0.0073 | 0.0109 | 0.90 | 0.91 | 4 | 4 | 51.26 | 64.54
$LP(0,1)/20^{0.25}$ |  |  |  |  |  |  |  |  |  |
600 | 0.0758 | 0.0642 | 0.0217 | 0.0295 | 0.64 | 0.71 | 4 | 4 | 18.53 | 13.03
2000 | 0.0350 | 0.0327 | 0.0140 | 0.0173 | 0.84 | 0.85 | 4 | 4 | 35.70 | 41.91
6000 | 0.0199 | 0.0178 | 0.0078 | 0.0103 | 0.91 | 0.92 | 4 | 4 | 79.76 | 79.63

##  5 Applications to real data

We now apply our procedure to an array-based Comparative Genomic Hybridization
(array CGH) data. Array CGH is a powerful technology for measuring copy
numbers at thousands of loci simultaneously. The output of array CGH
experiment is usually a long vector, spanning each chromosome, recording the
$\mbox{log}_{2}$ ratios of the normalized probe intensities from the test
samples vs. the reference samples. These ratios of intensities are used to
approximate the ratios of DNA copy numbers in the test samples vs. the
reference samples. A $\mbox{log}_{2}$ ratio far from $0$ (either positive or
negative) indicates a possible DNA copy number amplification or deletion for
the probe. Identification of chromosomal alteration regions will provide
valuable information to elucidate disease etiology and to discover novel
disease related genes.

In the study conducted by [28], cDNA microarray CGH was profiled across $6691$
mapped human genes in $44$ breast tumor samples and $10$ breast cancer cell
lines. The raw data can be downloaded from the PNAS website
(https://www.pnas.org/content/suppl/2002/09/23/162471999.DC1/4719CopyNoDatasetLegend.html).
We picked the breast cancer cell line BT474 as an example, and applied our
method to detect DNA copy number amplification. Details of one realization are
in the Supplementary Material and the results are presented in Table 7

Table 7: Results based on one sequence and multiple sequence with one-sided
test

One realization | Multiple realizations
---|---
Chromosome | beginning | ending | Chromosome | beginning | ending
number | loci | loci | number | loci | loci
11 | 68434309 | 81603744 | 11 | 46512342 | 81603744
|  |  | 14 | 16522721 | 106822024
|  |  | 15 | 17156123 | 18891425
17 | 28552955 | 82172608 | 17 | 28552955 | 42040770
20 | 43585793 | 66314778 | 20 | 44457372 | 66314778
21 | 12430025 | 15830914 | 21 | 12430025 | 15889676

Note: Chromosome 14 and 15 are connected as one cluster with very short
segments in chromosome 15 with multiple realizations analysis, chromosome 20
and 21 are connected as one cluster with both one realization and multiple
realizations analysis.

For multiple realization analysis, we consider the one-sided test using cell
line 1 in addition to BT474 for analysis. We compute $\hat{\mu}_{i}$ and
$\hat{\sigma}_{i}^{2}$ for $i=1,\ldots,p$, and test statistics
$\hat{R}_{j}=\sum_{l=j}^{j+k-1}\sqrt{2}\hat{\mu}_{l}/(\sum_{l=k}^{j+k-1}\hat{\sigma}_{l}^{2})^{1/2},j=1,\ldots,p-k+1.$
We use the same window length as in the one realization case $k=\lfloor
p^{1/2}\rfloor=78$, and compute
${Q}_{j}^{\star}=1({R}_{j}^{\star}>\gamma)+1(L_{j}^{\star}>\gamma)$ following
algorithm 3.1. Critical values $\gamma=3.8907$ and $\delta=1.0992$ are
obtained through the $0.95$th quantile of the empirical distribution of
$\mbox{max}_{1\leq j\leq p-k+1}{G}^{\star}_{j}$ and $\mbox{max}_{j\in
W_{1}}{G}^{\star}_{j},$ respectively, where
$G_{j}^{\star}=\sum_{l=j}^{j+k-1}\hat{\sigma}_{l}\eta_{l}/(\sum_{l=j}^{j+k-1}\hat{\sigma}^{2}_{l})^{1/2},j\in\mathbb{Z},\eta_{j}$
are i.i.d. $N(0,1)$ random variables and $W_{1}$ is the transition region
which includes indices $j$ such that the smoothed $\tilde{Q}_{j}^{\star}=1$.

The results are summarized in Table 7. We can see that four clustered regions
are detected by the multiple realizations analysis, three of which overlap
with those detected by one realization analysis, which shows that
amplifications in these genome regions are shared among the two breast cancer
patients. The identified chromosomal amplification regions are implicated in
the literature to harbor genes associated with breast cancer [27]. In cancer
studies, “passenger" mutations tend to occur more or less randomly throughout
the genome, and “driver" mutations tend to cluster and favor certain genome
positions containing functionally relevant genes. An important goal in the
analysis of tumor cell lines is to find the “driver" mutations, which play a
functional role in driving tumor progression [31]. Thus our analysis can
suggest followup studies and intervention strategies. We choose to conduct our
data analysis at the genome level, rather than at the chromosome level because
genome scale analysis allows the detection of copy number aberrations
involving entire chromosome arms, which might be missed in chromosome-level
analyses for which no actual changepoints exist.

Acknowledgments. We are grateful to two referees for their many helpful
comments. The research is partially supported by an NSF grant.

## References

  * [1] Benjamini, Y. and Heller, R. (2007). False discovery rates for spatial signals. J. Amer. Statist. Assoc., 102 1272–1281.
  * [2] Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. B, 57 289–300
  * [3] Bickel, P. J. and Rosenblatt, M. (1973) On some global measures of the deviations of density function estimates. Ann. Statist., 1 1071–1095.
  * [4] Cai, T.T. and Sun, W. 2009\. Simultaneous testing of groups hypotheses: finding needles in multiple haystacks. J. Amer. Statist. Assoc., 104 1467–1481.
  * [5] Cao, H. and Kosorok, M.R. (2011). Simultaneous critical values for $t$-tests in very high dimensions. Bernoulli,17 347–394.
  * [6] Cao, H. and Wu, W. B. (2015). Changepoint estimation: another look at multiple testing problems. Biometrika, 102 974–980.
  * [7] Cao, H., Chen, J. and Zhang, X.(2021). Optimal false discovery rate control for large scale multiple testing with auxiliary information. arXiv:2103.15311.
  * [8] Chan, Hock Peng and Walther, Guenther (2013). Detection with the scan and the average likelihood ratio Statistica Sinica, 23, 409–428,
  * [9] Chouldechova, A. (2014). False discovery rate control for spatial data. Ph.D. Dissertation.
  * [10] Du, L. and Zhang, C. (2014). Single-index modulated multiple testing. Ann. Statist., 42 1262–1311.
  * [11] Fan, J. (1996). Test of significance based on wavelet thresholding and Neyman’s truncation. J. Amer. Statist. Assoc., 91, 674–688.
  * [12] Fan, J., Hall, P. and Yao, Q. (2007). To how many simultaneous hypothesis tests can normal, student’s $t$ or bootstrap calibration be applied? J. Amer. Statist. Assoc., 102 1282–1288.
  * [13] Genovese, C. R., Roeder, K., and Wasserman, L. (2006). False discovery control with $p$-value weighting. Biometrika, 93 509–524.
  * [14] Hall, P., Kay, J. and Titterington, D. (1990). Asymptotically optimal difference-based estimation of variance in nonparametric regression. Biometrika, 77 521–528.
  * [15] Heller, R., Stanley, D., Yekutieli, D., Rubin, N. and Benjamini, Y.(2006) Cluster-based analysis of FMRI data. NeuroImage, 33 599–608.
  * [16] Hu, J. X., Zhao, H. and Zhou, H. (2010). False discovery rate control with groups. J. Amer. Statist. Assoc., 105 1215–1227.
  * [17] Lai, T.-L., Xing, H. and Zhang, N. (2008). Stochastic segmentation models for array-based comparative genomic hybridization data analysis. Biostatistics, 9 290–307.
  * [18] Lei, L. and Fithian, W. (2018). AdaPT: An interactive procedure for multiple testing with side information Journal of the Royal Statistical Society, Series B, 80, 649–679.
  * [19] Li, A. and Barber, R. F. (2015). Accumulation tests for FDR control in ordered hypothesis testing. J. Amer. Statist. Assoc., 112 837–849.
  * [20] Li, A. and Barber, R. F. (2019). Multiple testing with the structure adaptive Benjamini-Hochberg algorithm Journal of the Royal Statistical Society, Series B, 85, 45–74.
  * [21] Tansey, W., Koyejo, O., Poldrack, R. A. and Scott, J. G. (2018). False discovery rate smoothing. Journal of the American Statistical Association, 13, 1156–1171.
  * [22] Tibshirani, R. and Wang, P. (2008). Spatial smoothing and hot spot detection for CGH data using the fused lasso. Biostatistics, 9, 18–29.
  * [23] Liu, W.-D. (2015) Incorporation of sparsity information in large-scale multiple two-sample t tests. ArXiv e-prints, 2015. arXiv: 1410.4282.
  * [24] Liu, W.-D. and Shao, Q.-M. (2014). Phase transition and regularized bootstrap in large-scale $t$-tests with false discovery rate control. Ann. Statist., 42, 2003–2025.
  * [25] Pacifico, M. P., Genovese, C., Verdinelli, I. and Wasserman, L. (2004). False discovery control for random fields. J. Amer. Statist. Assoc., 99 1002–1014.
  * [26] Patil, G.P. and Taillie, C. (2003). Geographic and network surveillance via scan statistics for critical area detection. Statist. Sci., 18 457–465.
  * [27] Pollack, J. R., Perou, C. M., Alizadeh, A.A., Eisen, M. B., Pergamenschikov, A., Williams, C. F., Jeffrey, S. S., Botstein, D. and Brown, P. O. (1999). Genome-wide analysis of DNA copy-number changes using cDNA microarrays. Nature Genetics, 23 41–46.
  * [28] Pollack, J. R., Sorlie, T., Perou, C. M., Rees, C. A., Jeffrey, S. S., Lonning, P. E., Tibshirani, R., Botstein, D., Borresen-Dale, A.-L. and Brown, P. O. (2002). Microarray analysis reveals a major direct role of DNA copy number alteration in the transcriptional program of human breast tumors. Proc. Natl. Acad. Sci. U.S.A., 99 12963–12968.
  * [29] Pritchard, J. (2001). Are rare variants responsible for susceptibility to complex diseases? Am. J. Hum. Genet., 1 123–137.
  * [30] Shen, X., Huang, H. and Cressie, N. (2002). Nonparametric hypothesis testing for a spatial signal. J. Amer. Statist. Assoc., 97 1122–1140.
  * [31] Stratton, M. R., Campbell, P. J. and Futreal, P. A. (2009). The cancer genome. Nature, 458 719–724.
  * [32] Sun, W., Reich, B. J., Cai, T.T., Guindani, M. and Schwartzman, A. (2015). False discovery control in large-scale spatial multiple testing. J. R. Stat. Soc. B, 77 59–83.
  * [33] Szor, P. (2005). The art of computer virus research and defense. Addison-Wesley Professional.
  * [34] Siegmund, D., Zhang, N.R. and Yakir, B. (2011). False discovery rate for scanning statistics. Biometrika, 98 979–985.
  * [35] Tansey, W., Koyejo, O., Poldrack, R. A. and Scott, J. G. (2018). False discovery rate smoothing. Journal of the American Statistical Association, 13, 1156–1171.
  * [36] Tibshirani, R. and Wang, P. (2008). Spatial smoothing and hot spot detection for CGH data using the fused lasso. Biostatistics, 9, 18–29.
  * [37] Wu, W.-C. and Zhou, Z.(2018) Gradient-based structural change detection for nonstationary time series M-estimation. Ann. Statist., 46, 1197–1224.
  * [38] Wu, W.-C. and Zhou, Z.(2020) Multiscale jump testing and estimation under complex temporal dynamics. ArXiv e-prints, 2020, arXiv: 1909.06307.
  * [39] Yao, Q. (1993) Tests for change-points with epidemic alternatives. Biometrika, 80 179–191.
  * [40] Yao, Y.-C. (1988) Estimating the number of change-points via Schwarz’ criterion. Statistics and Probability Letters, 6 181–189.
  * [41] Zhang, C. M., Fan, J. and Yu, T. (2010) Multiple testing via FDRL for large scale imaging data. Ann. Statist. 39 613–642.
  * [42] Zhou, X., Carbonetto, P. and Stephens, M. (2013). Polygenic modeling with Bayesian sparse linear mixed models. PLoS Genet., 9 e1003264.