We study the problem of anomaly detection when a large number of data streams is observed. The streams themselves are not structured in any special way. In particular, there is no spatial or other proximity measure between streams, unlike in some areas like syndromic surveillance (Kulldorff et al., 2005) . We work under the following assumptions:
• Entirety. An anomalous stream is affected in its entirety. This assumption focuses the discussion on the unstructured nature of the setting. Although simplistic, this assumption is common in the literature (Aldosari and Moura, 2004; Lexa et al., 2004; Mei, 2008; Patwari and Hero, 2003; Thomopoulos et al., 1989; Yu et al., 2006) .
• Sparsity. An anomaly affects only a small fraction of the streams. This assumption focuses the discussion on the most compelling regimes.
Formally we consider n streams indexed by [n] ≡ {1, . . . , n}, that are observed over a time period [t] ≡ {1, . . . , t}. Specifically we observe X ≡ (X ij : i ∈ [n], j ∈ [t]), where X ij is the value from stream i ∈ [n] collected at time j ∈ [t]. Stream i ∈ [n] corresponds to a vector denoted by X i ≡ (X ij : j ∈ [t]).
Following standard practice, we place ourselves in a statistical decision theory framework and cast detection as a hypothesis testing problem. For simplicity, we assume the variables to be real valued and independent.
In this paper, we consider a stationary setting where the distribution of the observations remain unchanged over time. In detail, under the null hypothesis all X ij are independent and identically distributed (i.i.d.) with common distribution function denoted by F 0 . The alternative hypothesis is similar, but there is a subset of anomalous streams S ⊂ [n] such that for each i ∈ S the distribution of X ij , denoted by F i , stochastically dominates F 0 . In words, larger-than-usual values are observed in the streams indexed by S. The subset S is unknown and has no particular structure. Note that we are still assuming all the random variables X ij are independent. The sparsity assumption means the number of affected streams is believed to be small in comparison with the total number of streams, i.e. |S| n. Succinctly, our hypothesis testing problem is therefore: 
where for all i ∈ S, distribution F i stochastically dominates F 0 .
Despite its apparent simplicity, the above formulation and ensuing hypothesis test is quite suitable in a number of scenarios. For instance, in industrial settings the stream observations X i , i ∈ [n] might correspond to a failure metric measured from n machines. Under nominal conditions we expect these to be i.i.d. samples from some distribution, while under the anomalous conditions some of these machines will occasionally display higher-than-normal values for the failure metric. In more complex settings the observations might correspond to residual values, obtained by fitting a model to the streams in order to ensure these are in a common referential. For instance, in gene expression studies the streams might correspond to the normalized gene expression level of a number of replicates, so that under nominal conditions the expression levels are "equalized", and one is trying to determine if genes are differentially expressed under a specific environmental change. Obviously in this and similar settings some pre-processing of the data is needed to ensure all the streams are comparable. In Section 5.6 we consider yet a different application, where the proposed methodology can be applied directly to the daily counts of diagnosed COVID-19 cases in different municipalities in the Netherlands.
Contributions: Although a thorough characterization of the above testing scenario is far from trivial, it has been addressed quite extensively when the distribution F 0 is assumed known, and in particular when this is a normal distribution (Donoho and Jin, 2004; Ingster, 1997) . Knowledge of that distribution is required to calibrate the proposed tests (e.g., by Monte-Carlo simulation) and ensure the resulting p-value computation is done properly. When F 0 is unknown there are various ways one can proceed. A natural nonparametric approach is calibration by permutation. This approach is, in fact, quite standard in syndromic surveillance (Huang et al., 2007; Kulldorff et al., 2005 Kulldorff et al., , 2009 . It is also common in neuroimaging (Nichols and Holmes, 2002) and other contexts (Flenner and Hewer, 2011; Walther, 2010) .
In this work we propose a variant of the higher criticism statistic that relies on permutations in both its definition and its calibration, leading to an exact test. Furthermore, we show that, in the context of a one-parameter exponential family this test is asymptotically optimal, in the sense that it has essentially the same asymptotic performance as the best possible test with oracle knowledge of the model parameters. We provide practical recommendations so that the test is readily applicable in practice and demonstrate it has good performance in finite sample scenarios. In terms of testing power our procedure compares quite favorably with other proposed variants of higher criticism, particularly when streams are short (i.e., small t), since in its definition in does not rely on asymptotic approximations.
Related work: Early work on distribution-agnostic applications of the higher criticism can be found in Delaigle and Hall (2009) . The setting is different from ours, as the authors consider testing if a vector of observations, with unknown distribution, has the same mean as that of auxiliary data. Different approximation methods are discussed for the null distribution based on proximity measures or normal approximations, as well as the possibility of relaxing independence assumptions. Furthermore, a classification problem is discussed, in which (vector) samples from two populations are used to classify unlabeled vector samples to one of the corresponding classes. In a setting closer to ours, Delaigle et al. (2011) use the higher criticism test with studentized stream means and obtain a p-value by bootstrapping. Asymptotic results for the statistic are obtained, but controlling the type I error in finite samples is not discussed. The guarantees obtained are under a significantly different regime than the one we consider here, namely, the authors only assume that the null distribution has enough finite moments, but in turn require the stream length t to scale polynomially with the number of streams n, while our results hold for much shorter streams whose length scale poly-logarithmically.
The standard higher criticism statistic requires itself knowledge of the null distribution. In the absence of such knowledge several authors proposed the use of asymptotic approximations. In contrast, we use permutation-based estimates to construct the statistic itself, which seems to yield higher power in finite sample scenarios in comparison with approximationbased approaches, in addition to the strong theoretical guarantees we provide. Furthermore, computing the statistic using permutations comes at virtually no extra cost in comparison to the use of asymptotic approximations, when the statistic itself is already calibrated by permutation. Wu et al. (2014) consider a different setting than our, but with some important points-of-contact: in an association setting in genomics the authors aim to use the higher criticism to test if some covariates (gene SNP's) have non-null association to the outcome variables (traits), where a rare-weak assumption on the number and strength of the association is made. Calibration by permutation of the statistic is not the author's main interest, but is mentioned in the numerical section. The higher criticism statistic itself is computed based on asymptotic approximations. This type of methodology has been previously used in applied settings, for example in Sabatti et al. (2009) . In Donoho and Jin (2015) other works in this domain are reviewed, and the possibility of calibration by permutation is also mentioned, but their power properties are not established.
In Arias-Castro and Wang (2017) the authors consider the problem of detecting sparse heterogeneous mixtures from a nonparametric perspective, when the null distribution is assumed to be symmetric and the anomalous observations have a shift in mean. In Zou et al. (2017) the setting is similar to ours, but instead of hypothesis testing, their goal is to identify the set of anomalous streams. The authors then bound the probability of misidentifying the set of anomalous sequences. In Arias-Castro et al. (2018), a permutationbased scan statistic is proposed for the detection of structured anomalies (e.g., an interval) in a single stream, and shown to match to first order the performance of the scan statistic calibrated with (oracle) knowledge of the null distribution. Recently, Kurt et al. (2020) focuses on nonparametric change-point detection in a real-time setting and assumes there is a period of time guaranteed without anomalies. The authors show their decision statistics are asymptotically bounded under the null. Power properties are only shown empirically.
Organization: Section 1 introduces and motivates the problem. Section 2 briefly reviews known results for the normal model and motivates the use of a simple quantized variant of the higher criticism statistic. In Section 3 we introduce the class of distributional models we use as a benchmark for our theoretical power analysis and provide lower bounds on the anomalous signal strength for any test to be asymptotically powerful in the context of any one-parameter exponential family. Section 4 proposes a novel permutation higher criticism test, and we show that it is essentially asymptotically optimal in the context of the one-parameter exponential family. We also propose a permutation max test and establish its power properties. In Section 5 we examine the finite-sample performance of our methodology on simulated data, and also apply that methodology to a real data example. All the proofs are deferred to Section 7. Some supporting results are provided in the Appendix to ensure the manuscript is self-contained.
Notation: Throughout the paper we use standard asymptotic notation. Let n → ∞, then a n = O(b n ) when |a n /b n | is bounded, a n = o(b n ) when a n /b n → 0, and a n = ω(b n ) when b n = o(a n ). We also use probabilistic versions: a n = O P (b n ) when |a n /b n | is stochastically bounded 1 and a n = o P (b n ) when a n /b n converges to 0 in probability. Unless otherwise stated, we consider asymptotic behavior with respect to n → ∞.
This section considers a specific location model, which serves both as a benchmark and provides important insights needed for generalizations. These are well known results, but serve as a stepping stone to propose our methodology in the coming sections. Consider the setting where F 0 is the standard normal distribution with zero mean and unit variance. Under the alternative, there is a set S ⊂ [n] indexing the anomalous streams, and the corresponding observations have an elevated mean µ > 0 corresponding to the signal strength. Specifically, consider the following hypothesis test:
∼ N (0, 1) .
(2)
A natural question to ask is how large does µ need to be to ensure one can distinguish the two hypotheses. In this model it is obvious that the collection of stream means
are jointly sufficient. The setting then reduces to that of the normal means model (Donoho and Jin, 2004; Ingster, 1997) . Note that in those works the alternative hypothesis is a bit different, and often X ij are assumed to be i.i.d. samples from a sparse normal mixture. Nevertheless the asymptotic characterization of the optimal tests can be easily translated to the setting we are considering.
We take an asymptotic point of view. Let ψ n (X) : R nt → {0, 1} be a sequence of tests, where the outcome 1 indicates rejection of the null hypothesis. The risk of a test is simply the sum of type I and worst-case type II error, namely R(ψ n ) = P ∅ (ψ n (X) = 0) + max S:|S|=s P S (ψ n (X) = 1) ,
where we denote P S the probability under the alternative, with anomalous set S having signal strength µ. A sequence of tests is said to be asymptotically powerful if R(ψ n ) → 0 as n → ∞, and it is said to be asymptotically powerless if R(ψ n ) → 1.
To present an asymptotic characterization of the above hypothesis testing problem it is convenient to introduce the following parameterization. Let β ∈ [0, 1] be arbitrary (but fixed) and suppose S has size
where z denotes the smallest integer larger than z. With this parametrization, we can identify two main regimes:
• The dense regime: 0 ≤ β ≤ 1/2. In this scenario a simple test based on i∈[n] Y i (X) is optimal. As stated in the beginning our interest lies in settings where the anomaly affects only a reduced number of streams, and so we do not focus much on this regime. Whenever appropriate we make comments on how our results and methods extend to the dense regime.
• The sparse regime: 1/2 < β < 1. In this case a test based on i∈[n] Y i (X) will be asymptotically powerless. Let r > 0 be fixed and parameterize µ as
Define the detection boundary
It can be shown that, if r < ρ * (β) all tests are asymptotically powerless (see Ingster (1997) and Theorem 1). Furthermore, if r > ρ * (β) there are tests (introduced below) that are asymptotically powerful. Note that there are actually two sub-regimes within the sparse regime, namely:
The moderately sparse regime: 1/2 < β ≤ 3/4. Here Donoho and Jin (2004) show that the higher criticism test, based on the maximum of a normalized empirical process of the (Y i (X), i ∈ [n]) is asymptotically optimal. In fact, this is so in all regimes. Define the stream p-values as
where Φ denotes the cumulative distribution function of a standard normal random variable. Note that, under the null hypothesis these are i.i.d. samples from a uniform distribution in [0, 1]. One variant of the higher criticism rejects the null hypothesis for large values of
where α 0 ∈ (0, 1) and p (1) ≤ · · · ≤ p (n) are the ordered p-values.
The very sparse regime: 3/4 < β < 1. Here Donoho and Jin (2004) show that the test that rejects for large values of max i∈[n] Y i (X) achieves the detection boundary. This is not the case in the moderately sparse regime, as shown in . Expressed in terms of p-values, the max test is entirely equivalent to multiple testing with Bonferroni correction, rejecting the null hypothesis for small values of p (1) = min i∈[n] p i .
Remark 1. Note that the parameterizations in Equations (4) and (5) can be presented in a more general way using asymptotic notation. In particular, we can simply assume |S| = n 1−β+o(1) as n → ∞. For the sparse regime, it suffices to assume µ = 2(r + o(1)) t log(n) .
In this section we show that a quantized higher criticism statistic is also optimal in the normal location model. Quantizing the higher criticism statistic has been done before in different ways -for example, in Arias-Castro et al. (2011), Wu et al. (2014) and others. We discuss a quantization because the adaption and analysis of the higher criticism statistic, when using calibration by permutation, is easier for the quantized form introduced here. In contrast, adaptation to the permutation setting of the original higher criticism statistic in (7) is not straightforward as the computation of the test statistic itself requires knowledge of the null distribution -through the computation of the p-values in Equation (6).
We quantize a different variant of higher criticism. This variant was already introduced for analytical purposes in Donoho and Jin (2004) . Let q ≥ 0 and define
This is simply the number of streams whose average is larger or equal to 2(q/t) log(n). We refer to the latter value as the q-threshold. Define also p q , the probability that a stream mean from the null exceeds that same threshold:
Clearly, under the null hypothesis N q (X) is binomial with parameters n and p q , therefore we can standardize N q (X) and define the statistic
The V q (X) statistic, when maximized over q ∈ [0, ∞), is closely related to the higher criticism statistic (see Remark 4 below).
The following result summarizes several analytical steps in Donoho and Jin (2004) .
Proposition 1 (Donoho and Jin (2004) ). Let q ∈ [0, ∞) be fixed but arbitrary, and h n → ∞.
Furthermore, for S as in (4), when 1/2 < β < 1 (sparse regime), h n = n o(1) , and µ = 2(r/t) log(n) we have:
• If r < 1/4 and r > β − 1/2 then P S (V 4r (X) ≥ h n ) → 1.
The proof of Theorem 3, our main result, hinges crucially in showing that important quantities inside the proof of Proposition 1 have the same asymptotic characterization when considering calibration by permutation. The heart of the proof of Proposition 1 is the realization that, under the alternative, N q (X) is the sum of two independent binomial random variables; N q (X) ∼ Bin(n − s, p q ) + Bin(s, v q ), where
We will show a similar result in our setting, implying the result in this proposition can be used to complete the proof of Theorem 3. For completeness, a proof of Proposition 1 is also provided in Appendix A.1. Remark 2. In the dense regime, it is also possible to show that, when µ = 1 √ t n r−1/2 with r > β, we have P S (V 0 (X) ≥ h n ) → 1. This implies the higher criticism statistic is also optimal in that regime under the normal model.
Proposition 1 shows that, depending on the combination of sparsity |S| and signal strength µ we might want to consider a test based on V q (X) for different values of q. Interestingly, the value q = 1 suffices for the very sparse regime. For the moderately sparse regime the situation is more intricate, and it turns out the choice q = 4r is essentially optimal (unless r ≥ 1/4, as then q = 1 suffices). Obviously we do not know r, so the only way we can capitalize on this knowledge is to scan over a suitable selection of values for q, that hopefully include (or approximate) the optimal value. To keep the analysis simple and avoid the use of advanced tools in empirical processes we scan only over a coarse subset of values for q. The following result is an almost immediate consequence of Proposition 1.
Proposition 2 (Quantized HC). Let Q ⊆ [0, 1] be the set
where k n → ∞ and k n = n o(1) . Define the test statistic
Finally, let h n = ω( √ k n ) and h n = n o(1) , and consider the test that rejects the null hypothesis when T (X) ≥ h n . This test has vanishing type I error. When β > 1/2 and µ = 2(r/t) log(n) this test is asymptotically powerful provided r > ρ * (β).
The proof follows almost immediately by using a union of events bound over the grid elements. Under the alternative, the only intricate case is the moderately sparse regime. However, since the grid is becoming denser as n grows it will contain a value q that is close enough to the optimal value 4r. A proof is provided in Appendix A.2 for completeness. Remark 3. For the dense regime, one can prove that, if µ = 1 √ t n r−1/2 this test is asymptotically powerful provided r > β. This follows immediately using the results from Proposition 1, since the optimal value 0 is guaranteed to be in the grid. Remark 4. The V q (X) statistic is closely related to the more familiar higher criticism statistic as in (7) as
where i + is the largest value of i for which p (i) < 1/2. See Appendix A.3 and Lemma 7.
The results of the previous section were specific for the normal location model. Nevertheless one can envision somewhat natural extensions of the approach and test statistics to the more general context of the exponential family. When studying distribution-free tests, it is customary to compare them with parametric tests (Hettmansperger, 1984) . As in and Arias-Castro et al. (2018) , we consider one-parameter exponential models (in natural form) as parametric models, as these play an important role in the related literature (Donoho and Jin, 2004) .
One-parameter exponential models in natural form: Let F 0 be a probability distribution on the real line, such that all moments are finite. Let µ 0 and σ 2 0 denote respectively the mean and variance of F 0 . The distribution might be continuous, discrete or a combination of the two. In the exponential model there is a parameter θ i associated with each i ∈ [n], and the distribution F i ≡ F θ i is defined through its density f θ i with respect to F 0 : for θ ∈ [0, θ ), define f θ (x) = exp (θx − log ϕ 0 (θ)), where ϕ 0 (θ) = e θx dF 0 (x) and θ = sup{θ > 0 : ϕ 0 (θ) < ∞}, assumed to be strictly positive (and possibly infinite). In other words, f θ i denotes the Radon-Nykodym derivative of F θ i with respect to F 0 . Since a natural exponential family has the monotone likelihood ratio property it follows that F θ is stochastically increasing in θ (Lehmann and Romano, 2005, Lem 3.4.2) . Particular cases of this model are used in a variety of settings. For example, in many signal and image processing applications, a model with F θ corresponding to normal distribution with mean θ and a fixed variance is common. In syndromic surveillance (Kulldorff et al., 2005) , a model with F θ corresponding to a Poisson distribution is popular. Bernoulli models (Walther, 2010) are also a particular case of this class, with F θ corresponding to a Bernoulli distribution.
Let X ij be independent observations, where all observations within stream i are i.i.d. with distribution function F θ i , as defined above. Under the null hypothesis we have θ i = 0 for all i. Under the alternative, there is a subset of streams for which θ i > 0. Our hypothesis testing problem is therefore given by:
Keeping the same minimax stance as before, in the remainder our analysis assumes the anomalous observations have the same distribution so θ i = θ > 0 for i ∈ S. We keep the parametrization of the anomalous set as in (4), namely |S| = n 1−β for β ∈ (0, 1], and use the parameterization θ = 2r σ 2 0 t log n. Like it was the case for the normal model, a detection boundary defines the region in the (r, β) plane where any test will be asymptotically powerless. In fact, when the stream length t is large enough the detection boundary for the exponential family is identical to ρ * (β) defined earlier, as the following result shows.
Theorem 1. Refer to the hypothesis testing problem in (10) and the parameterization |S| = n 1−β with β ∈ (1/2, 1) and θ = 2r σ 2 0 t log n with r > 0. Provided t = ω(log 3 (n)) any test will be asymptotically powerless when r < ρ * (β).
In the next section, we develop a permutation-based test that is able to attain optimal asymptotic performance, giving a converse to the statement of Theorem 1.
The quantized higher criticism statistic proposed in Section 2 is asymptotically optimal in the normal location model, but to compute it and calibrate it (i.e. compute a test p-value) one needs to assume knowledge of the null distribution. With that knowledge, a simple way to proceed would be to calibrate the test by Monte-Carlo simulation. The goal in this section is to construct a test (and test statistic) that can be computed and properly calibrated without the knowledge of the null distribution. A natural and appealing idea is to calibrate these tests by permutation. The idea of using permutations is not new in the context of comparing multiple groups (e.g, consider the classical Kruskal-Wallis test or any test based on runs). However, such approaches cannot deal with sparse alternatives, as in the setting we are considering (particularly when β > 1/2).
As can be expected, using calibration by permutation to calibrate a test based on the value of the maximum stream-mean yields a test that cannot be optimal for all values of β. Nevertheless it is insightful to consider this option for its simplicity.
The max test rejects the null hypothesis for large values of
. In words, X π corresponds to the permuted data over both streams and time. Calibration by permutation amounts to computing a p-value as
In words, this is just the proportion of times the value of the test statistic computed with permuted data is at least as large as the value of the test statistic computed on the original data. Although the computation of the p-value might seem computationally prohibitive, we can approximate it very accurately by considering instead a uniform sample from the set of permutations. We have the following characterization of the statistical properties of this test:
Theorem 2 (Permutation max test). Let α ∈ (0, 1) be arbitrary and refer to the hypothesis testing problem in (10). Let θ = 2r σ 2 0 t log(n) with r > 0. Consider the test that rejects H 0 if P max-perm (X) ≤ α. This test has level at most α. Furthermore let t = ω(log 3 (n)) and β > 0. This test has power converging to one provided r >
This theorem shows the max test calibrated by permutation is only guaranteed to be optimal when β ≥ 3/4, as expected. It is important to note the stream length t starts playing a more prominent role when considering calibration by permutation. Clearly, if t = 1 it is impossible to calibrate such a test by permutation, as the test statistic is actually invariant under permutations in that case. One expects that, for very small values of t there might therefore be a loss of power in comparison to an oracle test calibrated with knowledge of the null distribution. The requirement for t in the above theorem is needed to ensure the stream means have approximately Gaussian tails. For exponential tails this scaling seems necessary. We conjecture that if the tail of F 0 is sub-Gaussian the requirement can be relaxed to t = ω(log 2 (n)) and if it has bounded support t = ω(log n) will suffice. It might seem the above proposition should follow directly from the results in Arias-Castro et al. (2018) as one can interpret the max test as a scan test over intervals of length t within nt observations, but a more careful analysis is needed to get a sharp characterization of the detection boundary.
Although the simplicity of the max test is appealing, the test is far from optimal in the moderately sparse regime (or the dense regime for that matter). Therefore we turn our attention to the possibility of calibrating the higher criticism statistic by permutation. Note that the statistic in Theorem 2 critically depends on the specification of the parameters in our normal setting. Now, since the mean and variance of the null distribution are unknown this needs to be taken into account. Definê
We must also adapt the normalizing term p q accordingly. In principle, the computation of this probability requires knowledge of F 0 . Instead, we obtain a p-value by permutation as followsP
Alternatively one might consider an asymptotic normal approximation to computeP q instead, as used in Sabatti et al. (2009) and Wu et al. (2014) . However, our numerical experiments indicate the permutation-based approach in (13) leads to tests with higher power.
Finally, define the counterpart of V q (X):
where we convention that 0/0 = 0. We are now ready to state the main result.
Theorem 3 (Permutation-based higher criticism test). Let Q ⊆ [0, 1] be the set
with k n → ∞ and |Q| = k n = n o(1) . Define the test statistic:
Finally, define the higher criticism permutation p-value as
Let α ∈ (0, 1) be arbitrary and refer to the hypothesis testing problem in (10). Consider the test that rejects H 0 if P perm-hc (X) ≤ α. This test has level at most α. Furthermore, let β > 1/2 and θ = 2r σ 2
This result shows that, provided the stream length t is large enough this higher criticism test calibrated without knowledge of the null distribution is asymptotically as powerful as any oracle test with that knowledge. This comes at a higher computational cost, but it is not prohibitive; see Section 5.1. Just as with the permutation max test, the requirement on t is driven by the tail behavior of the null distribution.
Remark 5. For the dense case (i.e. β ≤ 1/2) one can show that, if θ = σ 0 √ t n r−1/2 with r > 0 and t/n = Ω(1), then the permutation higher criticism test has power converging to one provided r > β.
This section thoroughly explores the potential of the proposed procedure. We begin by making some important computational remarks, followed by results on simulated data, and concluding with the application of the methodology to a real dataset.
The p-value of the permutation test in Theorem 3 is calculated by considering all possible data permutations. Naturally, this becomes prohibitive even for moderately large datasets, and instead a sample of the permutation distribution can be used to estimate P perm-hc (X).
At first glance computation and calibration of this test statistic might still seem prohibitive: for each permutation π taken from Π, one computes T (X π ), which requires us to compute P q (X π ). Naturally, we can estimate this quantity in a similar fashion using a random set of permutations. However, this would mean that every permutation sample π used to compute T (X π ) requires us to take a new set of permutation samples, which is disastrous computationally. However, note that P q (X) is invariant under permutations, namely P q (X) = P q (X π ). Thus we can estimate P q (X) using a random sample of permutations, and then use these estimates to compute T (X π ). Moreover, numerical experiments indicate we can use the same set of permutations to both estimate P q (X) and compute T (X), without loss of performance or raising the type I error. This approach has therefore been used in the ensuing numerical experiments.
For our experimental results, we consider two different models: the normal location model, discussed in Section 2 and a very standard benchmark, and a model where the underlying distributions are exponential. In the second model, the null distribution F 0 is exponential with mean 1/λ 0 (where we chose λ 0 = 1.5). Following the parametrization set up in Section 3 the observations in the anomalous streams are exponentially distributed with mean 1/(λ 0 −θ), where θ ∈ [0, θ * ) where θ * = λ 0 . This second setting allows us to explore the influence of heavier tails. In fact, the exponential tails are the heaviest tails for which we provide asymptotic guarantees.
In all simulations considered we assume the observations in the anomalous streams have distribution F θ , where
with τ ≥ 0. In words, we consider a signal strength that is a multiple of the asymptotically minimal signal strength ρ * (β) needed for detection. Asymptotically we know that when τ < 1 all tests are powerless, and when τ > 1 the permutation higher criticism test is powerful. In all simulations, test are calibrated at a significance level of 5%.
Although the grid Q proposed in Theorem 3 is valid asymptotically, it can leave out valuable information in finite samples, where the global maximum of V q (X) can be attained for q > 1.
To safeguard the approach against this issue one can instead use a data-dependent grid, such that the grid contains a sufficiently large point q * for whichN q * (X) = 0. A specific choice is
(Recall the definitions given in (12).) The statement in Theorem 3 still holds for this grid, as long as the number of elements in the grid k n is not larger than n o(1) , and the grid Q approximates the range [0, 1] with increasing accuracy as n increases. To satisfy both requirements, a suitable choice for the number of elements k n is given by:
where d n is a slowly diverging function. Note that 1/d n is the spacing between consecutive grid points. We show in this section that d n = log n is a suitable default choice. Nevertheless, provided d n → ∞ and d n = n o(1) , Lemma 4 ensures that k n = n o(1) with high probability. It is then possible to adapt the proof of Theorem 3 to get the same guarantees for this grid choice as explained, but this clutters the proof with undesirable technicalities.
It has been suggested in related works Wu et al., 2014 ) that scanning over a discrete grid for this type of statistic is merely a proof artifact, and the results are also valid when q in (16) is optimized over the continuous interval [0, ∞). We have conducted extensive numerical experiments regarding the choice of grid and observed no deteriorating performance for extremely fine grid choices. However, since we use permutation-based es-timatesP q (X π ) the number of specific values of q we need to consider to optimize V q over [0, ∞) becomes significantly larger than n, in principle as large as nB where B is the number of permutations considered. So, from a practical standpoint it is still rather useful to consider the proposed quantization.
We thus focus on the finite choice as in (18) in all experiments. To recommend the reader a choice of d n , we simulate data from both models outlined in Section 5.2 and compare the performance of the tests for different choices of d n as in (19). We use n = 10 3 streams of length t = (log n) 2 = 48. For the experiments shown we took |S| = 12, (i.e. β ≈ 0.64, moderately sparse regime), but through extensive experimentation we observed our conclusions are not sensitive to the sparsity. The results are given in Figure 1 . We observe that a spacing of 1/d n = 1/ log(n) is a reasonable choice. In all remaining experiments, we therefore extend the grid as in (18) and (19) with d n = log(n).
We now compare our methodology with the oracle test making full use of the null distribution. The oracle test statistic is defined in (9) for the normal model. For the exponential model the stream averages are Gamma-distributed, and so one uses the quantiles of that distribution to compute p q instead. The oracle-test is then calibrated by Monte-Carlo simulation. For this calibration, 10 4 simulations were used. In the presentation we use HC-permutation and HC-oracle to refer to the two tests.
We compare the oracle test with our methodology with respect to two aspects: (19). The power of the test is depicted as a function of a τ in the parameterization in (17). We took n = 10 3 streams with length t = 48 and set |S| = 12. For each test 10 3 permutations were used. Each level was repeated 10 3 times, resulting in the 95% confidence bars depicted.
• Stream lengths. Given that calibration by permutation is meaningless for t = 1 we expect the performance gap between an oracle and the proposed test to be large for small streams. However, we observe the finite-sample performance is already quite comparable to the oracle test even for very small streams.
• Performance loss. Our test makes no use of knowledge of the null distribution. While asymptotically equivalent to first order, we expect a performance loss in finite samples compared to the oracle.
We use n = 10 3 streams in our simulations. To study the influence of stream length we considered |S| = 12 (i.e. β ≈ 0.65, moderately sparse regime) but the results are qualitatively identical for other choices of sparsity. The signal strength was chosen slightly above the asymptotic detectability threshold, but ensuring the difference in performance was highlighted. For the exponential model very short streams would lead to a value θ ≥ λ 0 in the parameterization in (17), so the shortest stream considered is t = 4. Figure 2 show that the performance is comparable with respect to the oracle tests even at very short lengths, suggesting that finite-sample performance is much more forgiving than what our theoretical characterization might suggest, as there is only a very small loss of power in comparison to an oracle. Naturally, the permutation test has no power for streams of length t = 1. Note that the strength of the anomalous signal is decreasing in t, such that in the normal model the oracle test performs identically for any value of t. In contrast, for the exponential model, the difference in null and alternative means is θ/(λ 0 (λ 0 − θ)) and is larger for smaller streams, so all tests perform worse at larger values of t.
We now turn our attention to the power characteristics of the proposed tests. We fix n = 10 3 (17) with a fixed level for τ . Note that the strength of the anomalous signal is decreasing in t as discussed in the text. For the exponential model, the parametrization of the alternative is undefined for t ≤ 3, so the t-axis starts at t = 4. We took n = 10 3 streams and set |S| = 12. For the permutation tests, 10 3 permutations were used. For the oracle test, the Monte-Carlo calibration is based on 10 4 samples. Each level was repeated 10 3 times, resulting in the 95% confidence bars depicted. and t = (log n) 2 = 48 and consider |S| = 12 (i.e. β ≈ 0.64, moderately sparse regime), and |S| = 3 (i.e. β ≈ 0.84, very sparse regime). The results are plotted in Figure 3 , and we see that the proposed permutation higher criticism test has only minute loss of power when compared to the oracle test.
As mentioned earlier, several authors have considered practical variants of higher criticism where the permutation stream means are assumed to be approximately normal, and therefore the test statistic is computed based on a normal approximation (e.g., see Sabatti et al. (2009); Wu et al. (2014) and references therein). In essence, instead of computing (or, more appropriately, estimating)P q (X) as we do in Equation (13), one could instead compute:
, which arises from the assumption that the standardized permutation stream means Y i (X π ) are approximately standard normal. The statistic can then be computed like before, and calibrated by permutation. We refer to this methodology as HC-approximation in what follows. We also present results of the max test.
To compare the two approaches, we consider the exponential setting outlined in Section 5.2 with n = 100 streams and t = 4 and t = 6 as stream length. The reason to take a smaller (17). We took n = 10 3 streams with length t = 48. For the permutation tests, 10 3 permutations were used. For the oracle test, the Monte-Carlo calibration is based on 10 4 samples. Each level was repeated 10 3 times, resulting in the 95% confidence bars depicted. (17). We took n = 10 2 streams and set |S| = 12. For each test 10 3 permutations were used. Each level was repeated 10 3 times, giving rise to the 95% confidence bars depicted.
number of streams than before is to clearly emphasize the differences between the normal approximation and permutation approaches. We set |S| = 12 (i.e. β ≈ 0.64), but similar results can be obtained in the other regimes. Figure 4 depicts the results.
We see that there is a clear improvement in finite-sample performance when computingP q (X) by permutation, as opposed to the normal approximation. As expected, the difference in power decreases as the streams get longer, since the normal approximation becomes more accurate. In any case, the approximation p Φ q underestimates the true value of p q when q is large, as the tail of a Gamma distribution is much heavier than that of a normal distribution. SinceP q (X) does not rely on asymptotic approximations it is able to more accurately estimate p q , resulting on a more powerful test. Extensive experimentation showed that the permutation-based approach appears to be always superior when the underlying model is not normal.
After an initial explosive spread of COVID-19 several countries are entering a more stable regime. When in such a regime, it is of high importance to policy makers to detect signs of a new impending outbreak: for instance, when a few municipalities of the country reporting higher-than-usual number of cases.
When the virus is circulating in a reasonably controlled fashion, the dynamics of the disease transmission are homogeneous across different municipalities and a relatively short time-frame is considered, it is not unreasonable to model the daily diagnoses per capita of each municipality as originating i.i.d. from some unknown distribution. However, when the virus is spreading fast in a few municipalities those will display a slightly higher-than-usual number of daily-cases. As such, the framework we have introduced in this paper can be useful here to detect signs of local outbreaks, as these would lead to rejection of our null hypothesis (1).
We consider data from the Netherlands. In this country, it is not unreasonable to assume municipalities are sufficiently comparable: the country is very small and population density is not too different among municipalities. We consider a timeframe of 5 days. Note that the incubation period of COVID-19 is believed to be between two and fourteen days -see, for example Lauer et al. (2020) .
We use data on new COVID-19 cases per 100.000 inhabitants from 13th of March up until 10th of August 2020, for each of the 355 municipalities of the Netherlands. This data was processed from two sources; the data on the number of diagnoses (uncorrected for municipality population) was retrieved from the Dutch national institute for public health and the environment (RIVM) at https://data.rivm.nl 2 , and the data on municipality population was retrieved from the Dutch central agency for statistics (CBS) at https://opendata.cbs.nl 3 . These were combined to obtain the number of COVID-19 cases per 100.000 inhabitants.
When tests are done on the raw data we often reject the null hypothesis due to one or a few very large observations. While these cases are important in the context of the application, one needs no powerful test to mark them as anomalous as their abnormality is so clear. Especially in the context of COVID-19, these "clear" outliers will be investigated regardless. A more interesting question is then if, after removing "clear" outliers, the other municipalities still contain anomalous signals or not.
To remove "clear" outliers in a rigorous way, we use the max test as described in Theorem 2. We obtain the 95% quantile of the permutation maximum stream mean distribution, and mark all observations exceeding this threshold as "clear" anomalies. Then, we apply our permutation higher criticism test as in Theorem 3, as well as the higher criticism test using a normal approximation as in Section 5.5 on the remaining data to detect possible signals.
We apply our test sequentially over a sliding window of five days across the time domain, and obtain results for each window. The sequential p-values obtained in this fashion are plotted in Figure 5 , along with the virus' nationwide progression. Obviously, the obtained p-values are dependent and care must be taken when interpreting them -the figure is given merely to aid the presentation.
As can be seen in Figure 5 , our test results in much smaller p-values than the approximation method at each window. There are cases where our test indicates anomalies in seemingly stable periods. In these periods, while hard to see in aggregated data, we thus have some evidence that some municipalities have larger-than-usual values. This does not necessarily lead to a nationwide outbreak, since local measures, such as a restricting access to specific nursery homes, might have been taken to prevent further spread. As our data is aggregated per municipality and local measures are very hard to identify, we are not able to take this into account in our analysis.
Nonetheless, it would be interesting to identify the anomalous municipalities contributing to the low p-value of our test. A natural, but perhaps naive approach here would be to compute p-values of each stream by permutation, and apply suitable multiple-testing methods to control the family-wise error rate or the false discovery rate. Note, however, that such permutation p-values will be dependent, so standard methods as Benjamini and Hochberg (1995) might lead to misleading results. Furthermore, note that the signal picked up by our test may very well be too faint to be reliably estimated, but large enough to be detectable. Donoho and Jin (2004) include a discussion on this in the normal model we consider. This remains an interesting area for future research. Work related to this problem can be found, for example, in Bartroff and Song (2016) ; Romano et al. (2008) ; Romano and Wolf (2007) .
This paper considered the problem of anomaly detection when one observes multiple streams of data. This problem is closely related to the detection of sparse mixtures, where it is known that a test based on a higher criticism statistics is asymptotically optimal. Computation of this test statistics and calibration of the ensuing test requires the knowledge of the null (nominal) distribution. In this work we propose a distribution-free version of the statistic that is based on permutation ideas, and we show that it is also asymptotically optimal in the context of one-parameter exponential families. The proposed test is exact, and appears to be superior to alternative proposals that make use of asymptotic approximations. The numerical results in Section 5.4 show that the calibration by permutation results in nearly no loss of power in finite samples compared to an oracle test constructed and calibrated with full distribution knowledge. Its usefulness on real statistical problem is demonstrated in Section 5.6.
Our asymptotic results are proven for the class of exponential models. In Delaigle et al. (2011) a less restrictive class of models is considered, for which the proposed testing methodology attains a similar detection boundary when streams are of length t = n. We conjecture that, when we consider a similar class of models (i.e., a location family where the base distribution satisfies some tail-decaying condition) the detection boundary ρ * (β) is attained as well, provided the streams are long enough. Specifically, when the tails of the base distribution are decaying exponentially, a stream length t = ω(log 3 n) suffices, as shown in our work. For Pareto tails, we conjecture that at least t = ω(n) is needed to attain the ρ * (β) detection boundary, using large and moderate deviations theory results from Mikosch and Nagaev (1998) . One should keep in mind, however, that in an arbitrary location family the stream means are in general not a sufficient statistic for the location parameter, therefore the optimal detection boundary might not be characterized in the same way as for the normal model. For example, if F θ is the uniform distribution in [θ, θ + 1] a test based on the stream means will surely be suboptimal.
Further research could discuss the possibility of adapting the Berk-Jones statistic (Berk and Jones, 1979) in a similar fashion. Alternatively, one could investigate the possibility of replacing the observations by their ranks. This has numerical advantages and in some cases leads only to a minor loss of power (Arias-Castro et al., 2018) . A thorough analysis of this proposal is challenging as one must carefully account for the dependencies induced by the use of ranks, and remains an interesting avenue for future work.
In this section we provide the proofs of our main results. We begin with a simple technical lemma that greatly facilitates the presentation. The proof is a trivial consequence of standard tail-bounds for the normal approximation (Feller, 1968, Section 7 .1, Lemma 2).
Lemma 1. Let Φ be the cumulative distribution function of the standard normal distribution and x ≥ 0. Then (1)) log(n) = n −x+o(1) as n → ∞ .
Next, we argue that for the proofs of our main results, we can assume F 0 has zero mean and unit variance.
To prove the results in this paper it suffices to consider the case where the nominal distribution F 0 has zero mean and unit variance. To see this suppose F 0 has arbitrary mean and variance µ 0 and σ 2 0 . DefineF 0 (x) = F 0 (µ 0 + σ 0 x). It is easy to see the distributionF 0 has zero mean and unit variance. Using this we can easily re-parameterize the hypothesis test in (10).
Let X be a random variable with distribution F θ for some θ ∈ [0, θ * ) and defineX = X−µ 0 σ 0 . Define alsoφ 0 (θ) = eθ x dF 0 (x), the moment generating function ofF 0 . It is easy to check thatX has density with respect toF 0 given by exp θ x − log(φ(θ)) whereθ = σ 0 θ (equivalently θ = 1 σ 0θ ). Therefore, statements inθ pertaining a zero mean and unit variance distribution can be translated to a general distribution by simple multiplication by a factor 1/σ 0 .
Proof. Without loss of generality and as explained in Section 7.1 we assume that F 0 has mean zero and variance one, as this makes the arguments easier and less cluttered.
Let ψ(X) : R nt → {0, 1} denote an arbitrary test function. We begin by bounding the worst case risk of this test by the average risk, namely R(ψ) = P ∅ (ψ(X) = 0) + max S:|S|=s P S (ψ(X) = 1) ≥ P ∅ (ψ(X) = 0) + 1 n s S:|S|=s P S (ψ(X) = 1) .
The average risk can naturally be interpreted as the risk of testing the simple null hypothesis against a simple alternative, where S is chosen uniformly at random over the class of all subsets of [n] with cardinality s. Since we are doing a test between two simple hypotheses the optimal test (i.e., the test minimizing the average risk) is given by the Neyman-Pearson lemma, namely ψ(X) = 1 {L ≥ 1} where L is the likelihood ratio given by
The risk of this test can be easily expressed as 1 − 1 2 E (|L − 1|), where the expectation is with respect to the null hypothesis (so all X ij are i.i.d. with distribution F 0 ). To proceed we need to get an upper bound on E (|L − 1|). A simple, but often useful way to proceed is to use Jensen's inequality to get
where the equality above follows since E (L) = 1. This approach is generally referred to as the second moment method. To show any test is asymptotically powerless it suffices therefore to show that E (L 2 ) converges to one as n → ∞.
To simplify the presentation let S and S denote two independent random variables, and both independent from X. Both S and S are sampled uniformly from the set {S ⊂ [n] : |S| = s}. Then clearly L = E (exp (θX S − ts log(ϕ 0 (θ))) |X) and therefore E L 2 = E (exp (θX S − ts log(ϕ 0 (θ))) exp (θX S − ts log(ϕ 0 (θ)))) = E (exp (θ(X S + X S ) − 2ts log(ϕ 0 (θ)))) = E (exp (t|S ∩ S |(log ϕ 0 (2θ) − 2 log ϕ 0 (θ)))) .
For the last equality we used the fact that for all S and S we have X S +X S = 2X S∩S +X S S (in the previous expression denotes the symmetric set difference).
The beauty of the above result is that is reduces quantification of the risk to a statement about the moment generating function of the random variable K ≡ |S ∩ S |. Given the distribution S and S we conclude that K has an hypergeometric distribution with parameters (n, s, s), and therefore K is stochastically bounded from above by the binomial distribution with parameters (s, s n−s ). Using the well-know expression for the moment generating function of a binomial distribution we conclude that
Consider now the specific parameterizations of s and θ in the theorem statement. Note that when t = ω(log 3 n) then necessarily θ → 0, so we can conveniently use a Taylor expansion of the moment generating function ϕ 0 (θ) around θ = 0:
as θ → 0. Using the fact that F 0 has mean zero and unit variance we get ϕ 0 (θ) = 1 + 1 2 θ 2 + O(θ 3 ) as θ → 0. Simple asymptotic algebra yields that κ(θ) = 1 + θ 2 + O(θ 3 ). Since 1 + x ≤ e x we conclude that κ(θ) ≤ exp (θ 2 + O(θ 3 )).
When β > 1/2 we conclude that
The last expression converges to 0 provided 1 − 2β + 2r < 0 meaning that when r < β − 1/2 any test is asymptotically powerless. This lower bound is tight when β ∈ (1/2, 3/4] (the moderately sparse regime) but it is a bit loose for the very sparse regime. However, a modification of the above argument allows us to get a tight lower bound when β > 3/4.
The very sparse regime: the main limitation of the second moment method as presented above has to do with the fact that the likelihood ratio statistic L might take rather large values. Although this might be a rare occurrence, it can be enough to ensure the second moment is much larger than the first moment. A way to mitigate this issue is to consider a so-called truncated second moment method. Let Ω denote an arbitrary event and define the truncated likelihood ratioL ≡ L1 {Ω}. ClearlyL ≤ L and therefore
where we used the triangle inequality and the fact that E (L) = 1, followed by Jensen's inequality. Therefore, to show a test is powerless it suffices to show that both E L and E L 2 converge to one as n → ∞. The choice of event Ω is therefore quite crucial. In the present context we are going to consider the event
where η > 0 must be carefully chosen.
Truncated first moment: Note first that E L = E (L1 {Ω}) is the probability of Ω under the alternative hypothesis (where there is a set S of anomalous streams and S is chosen uniformly at random over the subsets of [n] with cardinality s). Given the symmetry of the definition of Ω we see that E L = P S (Ω) where S is an arbitrary set with cardinality s.
using the union bound in the last line. Using Lemma 6 we conclude that (1) ,
Therefore, when r ≤ 1 + η
This means the first truncated moment converges to one for any η > 0, provided
Truncated second moment: bounding this term requires significantly more work. Begin by noting that
The steps above mimic the derivation for the regular second moment, and the main difference is that we now need to consider the moment generating function of a truncated distribution, instead of the original distribution. Clearlỹ ϕ 0 (θ) ≤ ϕ 0 (θ) and so we conclude that E L 2 ≤ E (exp (t|S ∩ S |(logφ 0 (2θ) − 2 log ϕ 0 (θ)))) .
Defineκ(θ) ≡φ 0 (2θ)/ϕ 0 (θ) 2 . As before, to show the truncated moment converges to zero it suffices to show that s 2 n − s (κ(θ) t − 1) → 0 .
Note that, the argument based on the untruncated second moment method indicates all tests are powerless if r < β − 1/4. Since we are considering the case β ≥ 3/4 this means that it suffices to treat only the case where r ≥ 3/4 − 1/2 = 1/4. To get an upper bound onφ 0 (2θ) t we make use of the following technical result.
Lemma 2. Let X be a real-valued random variable and let f : R → [0, ∞) be one-to-one increasing and differentiable. Then, for any τ ∈ R,
To use this lemma we must get a good upper bound on P ∅ (Y 1 > x) for x ≤ τ (η). When x ≤ 0 we trivially bound this probability by one, and for 0 ≤ x < θ * we make use of a simple Chernoff bound. In a similar fashion to the proof of Lemma 6 we have
Let n be large enough so that τ (η) < θ * . Applying the above result and Lemma 2 we get
where in the second to last step we used the fact that t = ω(log 3 n). At this point note that
log n < 0 when r ≥ 1/4, provided we choose η > 0 sufficiently small. Therefore using Lemma 1 we conclude thatφ
With an analogous argument to the one used forφ(2θ) t one can show that (ϕ 0 (θ)) 2t = n 2r+o(1) . Therefore
The above converges to zero when 1 − 2β + 2r − ( √ 4r − √ 1 + η) 2 < 0. This is the case η is small enough and r < (1 − √ 1 − β) 2 , since in that case 1 − 2β + 2r − ( √ 4r − 1) 2 < 0, completing the proof.
Proof. By the arguments in Section 7.1, the proof continues under the assumption that F 0 has zero mean and unit variance.
Under the null and for a given (but arbitrary) permutation π ∈ Π it is clear that X π and X have exactly the same distribution. Therefore max i {Y i (X)} is uniformly distributed on the set {max i {Y i (X π }, π ∈ Π} (with multiplicities) conditionally on the order statistics of X. So, for a given α > 0
If there are no ties, the first inequality above is an equality, but with ties present the test becomes slightly more conservative. This argument is completely standard and for more details on permutation tests the reader is referred to (Lehmann and Romano, 2005) .
What remains to be proven is the behavior of the test under the alternative. Namely we must show that, provided r is large enough (as stated in the proposition) then for any α > 0
For convenience, let π be a uniformly distributed permutation of Π and let this be independent from X. We can rewrite our permutation p-value as a conditional probability:
To get a good upper-bound on the permutation p-value we use the following concentration inequality (see Shorack and Wellner (1986) and Arias-Castro et al. (2018) , for instance).
Lemma 3 (Bernstein bound for sampling without replacement). Let (Z 1 , . . . , Z m ) be sampled without replacement from the set {z 1 , . . . , z n }. Define z max = max j {z j }, z = 1 n n j=1 z j , Z = 1 m n j=1 Z j and σ 2 z = 1 n n j=1 (z j − z) 2 . Then, for all τ ≥ 0:
Using this lemma, we find that
The first inequality is a consequence of a simple union bound, and we used Lemma 3 in the second inequality.
At this point it is clear that, to control the p-value of our test we need to characterize the behavior of X, σ 2 X , max i,j X ij and max i Y i (X) under the alternative hypothesis. Since |S| = o(n) most of the elements of X are samples from the null distribution. Therefore we intuitively expect that X and σ 2 X should be good estimators for the mean and variance of F 0 . The behavior of the term max i,j X ij is a bit more delicate, but one can see that for the given parameterization of the null the dominant contribution is still given by the null distribution. In contrast, the term max i Y i (X) − X really depends on the alternative -the largest stream mean is surely driven by the anomalous observations. Formally, we can show the following result.
Lemma 4. Let β ∈ ( 1 2 , 1), θ = 2r(log n)/t with r > 0 and consider the alternative hypothesis in (10). Assume F 0 has zero mean and variance one and t = ω(log(n)). Then
(iii) Assume further that t = ω(log 3 (n)) and let ε > 0. Provided r > (
as n → ∞.
The first result of the lemma provides a rate at which the bound on our sample variance can decrease. For the analysis of the max test a much simpler result (already proved in Arias-Castro et al. (2018)) suffices: for any ε > 0 (i) implies that
Note also that the bound in Lemma 3 is monotonically decreasing in τ . This ensures that for ε > 0 and with probability tending to one under the alternative, provided r > (
Written differently, with probability tending to 1 under the alternative:
To ensure P max-perm (X) → 0, or equivalently, log P max-perm (X) → −∞ it suffices to ensure 1 c (log(n) + log(t)) 2(1 + ε) log(n) t < ε/2 .
However, since we assume t = ω(log 3 (n)) this is immediately satisfied, since the l.h.s. converges to zero.
We have just proved that, for ε > 0 and r > ( √ 1 + ε − √ 1 − β) 2 , P max-perm (X) → 0 as n → ∞. Since ε > 0 is arbitrary this implies the result in the theorem, concluding the proof.
Proof. By the arguments in Section 7.1, the proof continues under the assumption that F 0 has zero mean and unit variance. Furthermore the conservativeness of this test follows by the standard argument already presented in the proof of Theorem 2.
For the rest of the proof consider alternative hypothesis. We must show that P S (P perm-hc (X) < α) → 1 as n → ∞. Like before, we can write our permutation p-value as a conditional probability as follows:
where π is independent from X and uniformly distributed over Π. The first step is to understand and simplify the role of π in the above expression. This mirrors the analysis under the null hypothesis in the proof of Theorem 2, as we need to show that the values of the test statistic computed with permuted data are a good surrogate for the values of the test statistic under the null. However, the argument becomes more complex due to the dependencies introduced by the permutation. Like before, we will use the union bound for the max-operator in the permuted statistic, inducing a multiplicity by the grid-size:
To proceed recall that quantifyingV q (X π ) requires the quantification of two terms:N q (X π ) andP q (X π ). Note, however, thatP q (X) is invariant under permutations of X, and thereforê P q (X π ) =P q (X), as explained before. As such the only random quantity inside the probability symbol above (conditionally on X) isN q (X π ). Noting that E N q (X π ) X = nP q (X) we have
To apply Chebyshev's inequality we need the right-hand-side of the inequality inside the probability to be positive, andT (X) might be negative. In the latter case we simply bound the probability by one. Therefore, using Chebyshev's inequality we get
where we convention that 0/0 = 0. To continue, we must quantify the conditional variance ofN q (X π ). The permutation on X causes dependencies, but these are benign when realizing the conditional permuted stream means are negatively associated conditional on the data. Using Theorem 2.11 and the properties P 6 and P 4 in Joag-Dev and Proschan (1983), we find that Y i (X π )|X and Y j (X π )|X are negatively associated if i = j. For ease of notation, define
Then,
where we used the definition of negative association (Definition 2.1 from Joag-Dev and Proschan (1983) ). In conclusion we get the following simple bound for the p-value:
To continue the proof we must show thatT (X) is of order larger then k n = n o(1) . This mimics the approach in Proposition 1 and Theorem 2 under the alternative. Recall thatT (X) = max q∈QVq (X), so it suffices to get a good lower bound forV q (X). This requires controlling bothP q (X) (which is random), and the counting termN q (X) that is more complicated than before, as it involves estimates of the mean and variance of F 0 . Furthermore, we are no longer assuming normality. At this point we consider separately the cases β > 1/2 and β ≤ 1/2:
Let θ = 2r(log n)/t with r > 0. Provided t is large enough, the following lemma gives us a high-probability upper bound forP q (X).
Lemma 5. Consider the setting of Lemma 4 and let q > 0 and t = ω(log 3 (n)). There exists sequence g n → 0 such that under both the null and alternative hypothesis P P q (X) ≤ n −q+gn → 1 .
Letp q = n −q+gn where g n is the sequence in the lemma above. SinceV q (X) is decreasing in P q (X) (see (14)) we find that, under the alternative, with probability converging to 1
The next step is to get a high-probability lower-bound for the counting term. Our first aim is to bound the effect of the sample statistics. For q > 0 definẽ
By Lemma 4, we have X = O P 1/ √ nt and σ 2 X = 1 + O P ( 1/nt). Note that we could replace the terms log n/(nt) with g(n)/(nt), where g(n)/(nt) = o(1), but we chose this for concreteness. Therefore, we conclude thatN q (X) ≥Ñ q (X) with probability tending to 1. Thus, under the alternative, we have the following high-probability lower bound forV q (X):
where the inequality holds with probability tending to 1.
All the remains to be show is thatṼ q (X) is sufficiently large for an adequate value of q. It is clear thatÑ q (X) is a sum of two independent binomial random variables. Note that
Note that this probability is easily characterized by Lemma 6. Letp q ≡ w i,q if i / ∈ S and v q ≡ w i,q if i ∈ S, so that under the alternativẽ N q (X) ∼ Bin(n − s,p q ) + Bin(s,ṽ q ), Note that our probabilitiesp q andṽ q resemble the probabilities p q and v q considered in the proof of Proposition 1. Using Lemma 6 we immediately conclude thatp q = n −q+o (1) and
Thereforep q ,p q andṽ q all have the same asymptotic characterization as in the normal case, so we can continue in an analogous fashion as in the proof of Proposition 1 and Theorem 2 and conclude that when r > ρ * (β)
We have now all the pieces needed to bound the p-value in (26). The result just shown trivially implies thatT (X) → 1 and (k n + 1)/T 2 (X) → 0 with probability tending to one, and therefore
completing the proof.
Proof.
We have
where in the last line we changed the integration variable to x := f −1 (t).
Proof. Let us start by characterizing the overall sample mean X. By Chebyshev's inequality we have
as n → ∞. Now note that, under the alternative hypothesis
where X ∼ F θ . A Taylor expansion of the function inside the integral around θ = 0 yields (n)). Putting all this together yield the first result stated in (i).
For the second result in (i) note that
For the first term we see that 1 nt
Var (X ij ) + 1 nt i∈S,j∈ [t] Var
as n → ∞. In the above the variance of the anomalous streams was characterized with a Taylor expansion of θ around 0, similarly to what was done for the average term.
For the second term note first that all the moments of F 0 are finite, in particular the fourth moment. Therefore by Chebyshev's inequality we have
Finally, we know that X 2 = O P (1/nt) when β > 1/2. Putting everything together yields the second result stated in (i).
The argument needed to prove (ii) is the same already used in Arias-Castro et al. (2018), and presented here for completeness. Letting x > 0, a union bound gives
Now, let c ∈ (0, θ * − θ). We have that:
Therefore, with considerable slack, we can take x = 2 c log(nt) guarantee that both terms in (27) converge to zero. Together with the characterization of X in (i) we conclude that
To show part (iii) very different argument is needed as this is a lower-bound on the tail probability, rather than an upper bound. The following lemma gives a precise characterization of the tail probability.
Lemma 6. Consider the setting of Lemma 4. Let i ∈ S and let q n → q > r as n → ∞ and t = ω(log 3 n). Then
To show (iii) begin by noting that
Consider the case where ( √ 1 + ε − √ 1 − β) 2 < r < 1 + ε. Note that, in that case, we can use Lemma 6 with (28) which gives: (1) , where in last inequality we simply used the fact that log(1 + x) ≥ x. Finally, provided
The statement for r > 1 + ε follows immediately since this probability is monotonically increasing in r. To get the statement in (iii) we just need to use this result together with the characterization of X, concluding the proof.
Proof. We begin by using Lemma 3 to obtain an upper-bound onP q (X)
Now, part (i) of Lemma 4 implies that σ 2 X ≤ 1 + ε for ε > 0, with probability tending to 1. Using part (ii) of the same lemma we get
with probability converging to 1. Therefore, provided t = ω(log 3 (n)) there is a function g n → 0 such thatP q (X) ≤ exp (−(q + g n ) log(n)) .
To streamline the presentation let W 1 , . . . , W t be i.i.d. with distribution F θ and denote by ϕ θ (x) the moment generating function of F θ . Define also τ = (2q n /t) log n. We start by getting an upper bound for the said probability when r < q. A simple Chernoff bounding argument yields
We must now characterize ϕ θ (λ). First note that ϕ θ (λ) = ϕ 0 (λ + θ)/ϕ 0 (θ). Now, we develop a Taylor expansion of ϕ 0 (λ) around λ = 0, as we did in Equation (20). Note that F 0 has zero mean and unit variance. We obtain:
as λ → 0. A similar expansion can be developed for ϕ 0 (λ + θ) around λ + θ = 0. Combining all this yields
as both λ, θ → 0, where we used a Taylor expansion for the fraction around θ 2 /2+O(θ 3 ) = 0. This suggests the choice λ * = τ − θ, which is positive provided n is large enough since r < q. This choice yields the bound
since τ > θ and we used the basic inequality log(1 + x) ≤ x. When t = ω(log 3 n) the first term dominates, and therefore we conclude that
To lower-bound the probability in (29) we use a tilting argument. Let θ τ be such that F θτ has mean τ . Such a choice exists for n large enough and necessarily θ τ > θ for large n, since r < q. DefineW 1 , . . . ,W t to be i.i.d. with distribution F θτ . Then
With this change of measure we can conveniently use the central limit theorem to get a meaningful bound. Begin by noting that
where σ 2 θτ denotes the variance of F θτ . By the central limit theorem we know that
converges in distribution to a standard normal distribution and therefore the probability in the expression above converges to Φ(1) − Φ(0) ≈ 0.34 > 1/4. Therefore we conclude that, for n large enough
To control the remaining terms recall that ϕ 0 (λ) = 1 + λ 2 /2 + O(λ 3 ) as λ → 0 (see Equation (30)). Note also that τ = θ τ + O(θ 2 τ ), which implies (after some manipulation) that θ τ = τ + O(τ 2 ). Finally, note that both τ and θ have the same order of magnitude. Putting all this together we conclude that
For the other term note that σ θτ = 1 + o(1), and therefore σ θτ / √ t = o(θ). This implies that
In conclusion
When t = ω(log 3 n) the term (τ − θ) 2 term dominates, and we see we get the asymptotic behavior as in the upper bound, concluding the proof.
For the case r ≥ q we see that necessarily P 1 t j∈[t] W j ≥ τ ≥ n o(1) , so the tail probability cannot be extremely small. In fact, when r > q this probability will be lower bounded by a constant.
Proof. For convenience, define |S| = s. The first statement in the proposition is a simple consequence of Chebyshev's inequality. Under the null hypothesis N q (X) ∼ Binomial(n, p q ), therefore
Since h n → ∞ the statement of the proposition follows.
Consider first the sparse regime 1/2 < β < 1. Under the alternative hypothesis N q (X) is a sum of two independent binomial random variables, that is, N q (X) ∼ Bin(n − s, p q ) + Bin(s, v q ). Each of the variables corresponds respectively to the counts of the null and anomalous stream means exceeding the threshold 2q log(n), where v q = 1 − Φ 2 log(n)( √ q − √ r) .
In words, v q is the probability an anomalous stream mean exceeds the threshold 2q log(n). For convenience of presentation define also a n ≡ s(v q − p q ) − h n np q (1 − p q ) .
Using this definition we see that P S (V q (X) < h n ) = P S N q (X) ≤ h n np q (1 − p q ) + np q = P S N q (X) − E (N q (X)) ≤ h n np q (1 − p q ) + s(p q − v q ) = P S (N q (X) − E (N q (X)) ≤ −a n ) = P S (− (N q (X) − E (N q (X))) ≥ a n ) ≤ P S (|N q (X) − E (N q (X))| ≥ a n ) .
When a n > 0 we can easily bound the above expression using Chebyshev's inequality. To make the presentation and derivations clear and simple we make use of asymptotic notation.
All the statements below are taken when n → ∞. Using Lemma 1 we see that p q = n −q+o (1) and
Now, since h n = n o(1) , we have that a n is given by: a n = n 1−β−( √ q− √ r) 2 +o(1) − n max{ 1−q 2 ,1−β−q}+o (1) if r < q n 1−β+o(1) − n max{ 1−q 2 ,1−β−q,1−β−( √ q− √ r) 2 }+o (1) if r ≥ q .
We can now easily retrieve the conditions under which a n > 0 for large enough n, namely the exponent of the positive terms must be larger than the exponent of the negative terms in the right-hand-side of Equation (33). As a n differs depending on whether r < q or r ≥ q we obtain a different set of conditions of each of these cases. However, many of such conditions are satisfied trivially, and for each case we obtain only one single nontrivial condition:
If r < q then a n → ∞ provided 1
If r ≥ q then a n → ∞ provided 1 − β > 1 − q 2 .
Under these conditions we can continue from (31) via Chebyshev's inequality and obtain P S (V q (X) < h n ) ≤ Var (N q (X)) a 2 n .
To analyze the bound, first observe that Var (N q (X)) = (n − s)p q (1 − p q ) + sv q (1 − v q ) ≤ np q + sv q (1 − v q ).
Note that v q and 1 − v q are equal for r < q and r > q respectively, and thus v q (1 − v q ) is equal for both cases asymptotically. Moreover, the case r = q can be merged with r = q by observing that 1 4 n 1−β = 1 4 n 1−β−( √ q− √ r) 2 = n 1−β−( √ q− √ r) 2 +o(1) . We therefore obtain for all r > 0:
Var (N q (X)) = n max{1−q,1−β−( √ q− √ r) 2 }+o(1) . (1) if r < q, n max{1−q,1−β−( √ q− √ r) 2 }−2(1−β)+o(1) if r ≥ q.
All that needs to be done at this point is to identify under which conditions the above bound converges to zero. Clearly, when r ≥ q this is the case. If r ≤ q we require that:
Therefore, for the results in the proposition to hold we require that conditions (34), (35), (36) and (37) are simultaneously satisfied. We can easily see that (34) and (36) are equivalent. After rewriting condition (36) slightly, we are then left with the following three conditions:
if r ≥ q, then 1 − β > 1 − q 2 .
if r < q,
if r < q, then 1 − β − ( √ q − √ r) 2 > 0.
At this point it is a matter of algebra to check these conditions are satisfied for the statements in the proposition. Note that we do not provide results for β = 1, as this would require taking q > 1.
Proof. Consider first the analysis under the null hypothesis. A simple application of the result in Proposition 1 and a union of events bound yields
where the last statement follows from the assumptions on h n and k n . Now we turn our attention to the alternative hypothesis, and must show that P S (T (X) ≥ h n ) → 1 , as n → ∞, under the conditions in the theorem for r and β. For the regimes where the choices q = 0 or q = 1 are optimal it is straightforward to ensure that P S (V q ≥ h n ) → 1, so the theorem result is obviously true. For instance, P S max q∈Q V q (X) ≥ h n ≥ P S (V 1 ≥ h n ) → 1 provided µ = 2r log(n) and r > (1 − √ 1 − β) 2 . Likewise a similar argument holds in the dense case, where the choice q = 0 is the best one.
The only situation that is slightly more intricate is that when we must consider the choice q = 4r, as in general that value is not in the grid Q. However, since k n → ∞ we will be able to find a value in Q that is close enough to 4r. Consider the case µ = 2r log(n), r < 1/4 and r > β −1/2. Let q * = min q∈Q |4r −q|. Clearly |q * −4r| ≤ 1 kn so q * = 4r +o(1). Therefore p q * = n −4r+o(1) and v q * = n −r+o(1) and we can follow exactly the same steps as in the proof of Proposition 1. Therefore
which is what we wanted to show.
The following lemma formalized the equivalence between the traditional higher-criticism statistic as in (7) and the V q (X) statistic introduced in (9).
where i + is the largest value of i for which p (i) < 1/2.
The proof follows by noting that the supremum of V q (X) is attained at a finite number of points. This result is essentially stated in Donoho and Jin (2004) but without a complete proof. Note that the right-hand-side is not entirely equal to the higher criticism statistic, but it is very similar (particularly when taking α 0 = 1/2, which is a common recommendation).
Taking α 0 means that essentially almost all streams that have p-values larger than 1/2 will be ignored.
Proof. We show our assertion (41) in two steps. First, we will show that the supremum of V q (X) is attained at a finite number of points. We will then show that the maximum of these points equals the HC statistic, if α 0 is chosen appropriately. For ease of notation, define:
such that we can succinctly write our classical higher-criticism statistic as the maximum over h i (p (i) ): max i∈[n]
√ n i n − p (i) p (i) (1 − p (i) ) = max i∈ [n] h i (p (i) )
To start, define:
for all Y i > 0. Denote X (i) the ordered stream means, with X (1) being the largest. Note that q i is monotone in Y i > 0, and thus the ordering is preserved.
Note that N q (X) is a stepwise, nonincreasing function in q, constant on intervals (q (i) , q (i+1) ]. Next, note that V q (X) is decreasing in p q , and p q is decreasing in q. Therefore, as N q (X) is constant on (q (i) , q (i+1) ], V q (X) is increasing on (q (i) , q (i+1) ].
Finally, as N q (X) is nonincreasing, V q (i) (X) ≥ V q (i) +δ (X). Therefore, the supremum of V q (X) must be attained at the points q i . Now, note that N q (i) (X) = i and p q (i) = p (i) , and thus V q (i) = h i (p (i) ). However, we only have these q (i) defined for Y i > 0, so: 