A super scalable algorithm for short segment detection
Ning Hao Yue Selena Niu Feifei Xiao Heping Zhang
stat.ME stat.AP stat.CO


equationsection







1]Ning Hao
1]Yue Selena Niu
2]Feifei Xiao
3]Heping Zhang
[1]University of Arizona
[2]University of South Carolina
[3]Yale University
                     



A super scalable algorithm for short segment detection
    [
    
======================================================





In many applications such as copy number variant (CNV) detection, the goal is to identify short segments on which the observations have different means or medians from the background. Those segments are usually short and hidden in a long sequence, and hence are very challenging to find. We study a super scalable short segment (4S) detection algorithm in this paper. This nonparametric method clusters the locations where the observations exceed a threshold for segment detection. It is computationally efficient and does not rely on Gaussian noise assumption. Moreover, we develop a framework to assign significance levels for detected segments. We demonstrate the advantages of our proposed method by theoretical, simulation, and real data studies.



Keywords: copy number variation, inference, nonparametric method, signal detection.



§ INTRODUCTION


Chromosome copy number variant (CNV) is a type of structural variation with abnormal copy number changes involving DNA fragments <cit.>. CNVs result in gains or losses of the genome, therefore interfering downstream functions of the DNA contents. Accounting for a substantial amount of genetic variation, CNVs are considered to be a risk factor for human diseases. Over the past decade, advances in genomic technologies have revealed that CNVs underlie many human diseases, including autism <cit.>, cancer <cit.>, schizophrenia <cit.>, and major depressive disorder <cit.>. It is fundamental to develop fast and accurate CNV detection tools.

A variety of statistical tools have been developed to discover structural changes in CNV data during last 20 years. Popular algorithms include circular binary segmentation <cit.>, the fused LASSO <cit.>, likelihood ratio selection <cit.>, and screening and ranking algorithm <cit.>. Some other change-point detection tools such as wild binary segmentation <cit.> and simultaneous multiscale changepoint estimator <cit.> can be also applied to CNV data. See <cit.> for a recent review on modern change-point analysis techniques. A majority of existing methods are based on Gaussian assumption, although quantile normalization <cit.> or local median transformation <cit.> can be used for normalization. The computational complexity is also of concern for some of the existing methods as the modern technologies produce extraordinarily big data. In spite of some fast algorithms <cit.>, few algorithms are known to possess both computational efficiency and solid theoretical foundation. Moreover, with only a few exceptions <cit.>, the existing methods focus on detection whereas not offering statistical inference. For these reasons, we develop a fast nonparametric method for CNV detection with theoretical foundation and the opportunity of conducting statistical inference.

In this paper, we model the CNVs as short segments with nonzero height parameters, which are sparsely hidden in a long sequence. The goal is to identify those segments with high probability and, moreover, to assess the significance levels for detected segments. In particular, we propose a scalable nonparametric algorithm for short segment detection. It depends on only the ranks of the absolute values of the measurements and hence requires minimal assumptions on the noise distribution. A short segment may be present when there are a large enough number of observations exceeding a certain threshold on a short segment; for instance, 8 on a segment of 10 observations are larger than the 99th percentile of the data. Following this idea, we implement a super scalable short segment (4S) detection algorithm to cluster the points to form a segment when such a phenomenon occurs.
The advantages of our method are fourfold. First, this nonparametric method requires minimal assumption on the noise distribution. Second, it is super fast as the core algorithm requires only O(n) operations to analyze a sequence of n measurements. In particular, it takes less than 2 seconds for our R codes to analyze 272 sequences with a range of about 34,000 measurements. Third, we establish a non-asymptotic theory to ensure the detection of all signal segments. Last but not least, our method can compute the significance level for each detected segment and offer a convenient approach to statistical inference.



§ METHOD



 §.§ Notations and the main idea

Let {X_j}_j=1^n be a sequence of random variables such that

    X_j=μ_j+_̌j,

where the height parameter vector =(μ_1,...,μ_n)^⊤ is sparse and {_̌j}_j=1^n are random noises with median 0.
Moreover, we assume that the nonzero entries of  are supported in the union of disjoint intervals =⋃_k=1^K I_k so that

    μ_j={[                   ν_k0, if j∈ I_k⊂ for some  k;;                      0,                 if  j∉. ].

Here we assume that I_k={ℓ_k,ℓ_k+1,...,r_k} and r_k<ℓ_k+1-1 for all k. Note that such a representation of  is unique and used throughout this paper for  or its estimator . For convenience and without confusion, we use the interval [ℓ_k,r_k] to imply the set of integers {ℓ_k,ℓ_k+1,...,r_k} when referred to a set of locations. We call those intervals segments. In particular, a signal segment is a segment where the height parameter is a nonzero constant. Let _I denote a subvector of  restricted to I⊂{1,...,n}. We denote by || the cardinality of a set . In particular, |I|=r-ℓ+1 for I=[ℓ,r].  and  denote vectors, (0,...,0)^⊤ and (1,...,1)^⊤, respectively.

Naturally, a primary goal for model (<ref>) is to identify the set of signal segments { I_k}_k=1^K. Moreover, while rarely done, it is useful to assign a significance level for each of the detected segments. In this paper, we will study both estimation and related inference problems on segment detection. Our strategy is to cluster “putative locations” using spatial information. In particular, we consider the set of positions _c={j:|X_j|>c, 1≤ j≤ n}, where the observations exceed a threshold c>0. Intuitively, for a properly chosen c and some segment Î, if |Î∩_c| is big enough compared with |Î|, it is likely that _Î.

To illustrate our idea, we consider a game of ball painting. Suppose that we start with n white balls in a row, and paint the jth ball with black color if |X_j|>c. Let m be the total number of black balls, which is much smaller than n. If we observe that there are a few black balls crowded in a short segment, e.g., `segment 1' illustrated in Figure <ref>, it is plausible that the height parameter is not zero in the segment. Our proposed algorithm can easily identify those segments. Also it may happen that in a neighborhood there is only a single black ball, e.g., `segment 2' in Figure <ref>. Then, we may not have strong evidence against =0. To put this intuition into a sound theoretical framework, it is imperative to evaluate the significance for each pattern. In fact, given the numbers of white and black balls in a short segment, we may calculate how likely a certain pattern appears in a sequence of length n with m black balls, when white and black balls are actually randomly placed. We will develop a framework of inference based on this idea in Section 2.3.







 §.§ Algorithm

To estimate , we propose a super scalable short segment (4S) detection algorithm which is described as follows.

Step 1: thresholding. Define _c={j:|X_j|>c}. That is, we collect the positions where the observations exceed a threshold.

Step 2: completion. Construct the completion set _c,d by the criterion that j∈_c if and only if there exist j_1, j_2∈_c, j_1≤ j_2≤ j_1+d+1 such that j_1≤ j≤ j_2. That is, we add the whole segment [j_1,j_2] into the completion set if the gap between j_1, j_2∈_c is small enough.


Write _c,d=⋃_k=1^K̃Ĩ_k where Ĩ_k =[ℓ̃_k, r̃_k] with r̃_k<ℓ̃_k+1-1. Note that this decomposition is unique.

Step 3: clean up.  We delete Ĩ_k from _c,d if |Ĩ_k|=r̃_k-ℓ̃_k+1≤ h, and obtain our final estimator _c,d,h=⋃_k=1^K̂Î_k. That is, we ignore the segments that are too short to be considered.

The whole procedure depends on three parameters c, d, and h. The choice of c is crucial and depends on applications. d and h are relatively more flexible as we can screen false positives using significance levels defined later. We may ignore the subscripts and simply refer to ,  and  when the sets obtained from the three steps above corresponding to c, d and h are clear in the context. Figure <ref> illustrates our procedure, where the location set obtained in each step is indicated by the positions of black balls.





 §.§ Theory: consistency and inference

Our goal is to identify the set of signal segments =⋃_k=1^K I_k with a false positive control. Here we say that I_k∈ is identified by an estimator =⋃_k=1^K̂Î_k if there is a unique Î_k'∈, such that Î_k'∩ I_k∅, and Î_k'∩ I_j=∅ for all I_j∈ and j k. Such an Î_k' is a true positive.
We define that  is identified by an estimator  if every I_k∈ is identified by . That is, there is a one-to-one correspondence between  and a subset of , and the K pairs under this correspondence are the only pairs with nonempty interaction among all segments in  and . See Figure <ref> for an illustration of the definition.



In our three-step procedure, the first two steps establish an estimator and the last one aims to delete obvious false positives. Our theory proceeds in two main steps. First, we characterize the non-asymptotic probability that the first two steps produce an estimator which successfully identifies . In order to identify I_k∈, we should ensure two conditions. Condition one is that, after step 1, the black balls are dense enough on I_k so that they do not split into two or more segments in step 2. Condition two is that, in the gap between I_k and I_k+1, the black balls are sparse enough so that the black balls on I_k and I_k+1 do not connect to a big segment. Theorem 1 addresses how to bound the probabilities of these two conditions for all k.

Second, we develop a framework of inference to control false positives. As a rough control, Theorem 2 gives an upper bound for the expected number of false positives if all segments of length one are deleted in step 3. In general, after step 2, it is not optimal to decide the likelihood of a detected segment being a false positive only by its length. Therefore, for each segment in , we check its original color pattern back in step 1 and calculate a p-value of this pattern under null hypothesis =. This assigns a significance level for each detected segment which helps control false positive. It is difficult to find the exact p-values. Lemma 2 offers a reasonable approximation.

To facilitate theoretical analysis, we assume that, in this subsection, {_̌j}_j=1^n are independent and identically distributed (IID) noises with median 0. Moreover, _̌j has a continuous density function f that is symmetric with respect to 0. Under this assumption, the black balls are randomly distributed for arbitrary threshold c when =.


Now we investigate when a signal segment can be detected by our algorithm. Let F be the cumulative distribution function of the noise density f. We use f_α=F^-1(α) to denote the α-percentile. Suppose that there is a segment I such that |I|=L, _I=ν and _I^c∩ H=, where H is a segment containing I such that I^c∩ H is the union of two segments both of which are of length D. Without loss of generality, we assume that ν=f_α>0 i.e. α>0.5. For a threshold 0<c=f_β≤ν, let us continue our game of ball painting and focus on this segment and its neighborhood. Recall that we paint the ball at position j with black color if and only if |X_j|>c=f_β. The following two events together can ensure that the segment I is identified by our method.

    =   {On I, there does not exist a sub-segment of min{d,L} consecutive white balls}
    _D    =   {On both sides of I, within distance  D,  there are  d  consecutive white balls}

Event  ensures that this segment can be detected as a whole segment while event _D controls the total length of the detected segment and makes sure that the detected segments are separated from each other.  and _D together guarantee that our algorithm identifies a segment Î such that Î∩ I∅, and Î∩ I'=∅ for any other signal segment I'. The following lemma gives non-asymptotic bounds for P() and P(_D).


Let H⊃ I be two segments such that _I=ν and _I^c∩ H=, |I|=L. I^c∩ H is the union of two segments both of which are of length D. Let β'=2β-1.
For 0<c=f_β≤ν and d>0, after the thresholding and completion steps, we have
                                   
    P()≥{[ 1-(L-d+2)2^-d-1,          if  d<L;          1-2^-L,         if  d≥ L ].
 and 
    P(_D)≥ 1-2(1-β'^d)^⌊D/d⌋.




Let ν_min=min_k|ν_k| be the minimal signal strength among all I_k's, L_min=min_k |I_k|, L_max=max_k |I_k| be the minimal and maximal lengths of signal segments, respectively, and D_min=min_k(ℓ_k+1-r_k-1) be the minimal gap between two signal segments. Define β_min=F(ν_min) so ν_min=f_β_min. Let β'_min=2β_min-1. Taking into account all signal segments in , the theorem below gives a lower probability bound for identifying  after first two steps.


With c=ν_min=f_β_min, d>0 and h=0, _c,d,h can identify all signal segments in  with probability at least

    1-Kmax{1/2(L_max-d+2),1}2^-min{d,L_min}-(K-1)(1-β_min'^d)^⌊D_min/d⌋.




The probability (<ref>) goes to 1 asymptotically if log K+log L_max≪min{d,L_min}→∞ and log K≪ D_min/d→∞ as n→∞.



Although Theorem 1 gives a theoretical guarantee to recover all signal segments with a large probability, there are some false positives. As an ad-hoc way, we may take h=2 or 3 to eliminate some obvious false positives. This clean up step is simple and helpful to delete isolated black balls. The Theorem below gives an upper bound on the number of false positives with a conservative choice h=1.


Assume =0 and |_c|=m. Then E|_c,d,h|≤ m (1-∏_k=1^dn-m+1-k/n-k) with h=1.



The expected number of false positive segments can be well controlled if both m/n and d are small. Next, we illustrate how to access the significance levels for the detected segments by our method, which is helpful to control false positives. Recall our estimator =⋃_k=1^K̂Î_k where Î_k=[ℓ̂_k,r̂_k]. For each Î_k, let s_k=|Î_k|, t_k=|Î_k ∩_c|, and m=|_c|. Now consider n balls in a row with m black and n-m white balls. Let _n,m,s_k,t_k be an event that there exists a segment of length s_k where at least t_k balls are black, in a sequence of n balls with m blacks ones. The p-value of Î_k can be defined as the probability of _n,m,s_k,t_k if the balls are randomly placed. This p-value can effectively control the false positives. However, it is challenging to find the exact formula to calculate the p-value. The lemma below gives an upper bound of P(_n,m,s_k,t_k).


P(_n,m,s_k,t_k)≤P̂(_n,m,s_k,t_k)=mP(Y≥ t_k-1) where Y follows a hypergeometric distribution with total population size n-1, number of success states m-1, and number of draws s_k-1.


This approximated p-value is useful to eliminate false positives.



 §.§ Implementation

Our proposed method is nonparametric and depends on only the rank of absolute measurements {|X_j|}_j=1^n. For a fixed triplet (c,d,h), it typically needs less than 3n operations to determine  when |_c|/n is small, say, less than 0.1. We need 2n operations to compare each measurement with the threshold to determine _c. Let =(w_1,...,w_m)^⊤ be a vector of locations in _c in an ascending order. In the completion step, we compare Δ=(w_2-w_1,...,w_m-w_m-1)^⊤ to a threshold. In particular, we declare that w_i and w_i+1 belong to different segments if and only if w_i+1-w_i>d. Let i_1,...,i_K̃ be those indices such as w_i_k+1-w_i_k>d.  consists of segments [w_1,w_i_1],...,[w_i_K̃,w_m]. We record the start and end points of each segment only. In the deletion step, we delete [w_i_k,w_i_k+1-1] if its length is not greater than h. The total operations can be controlled within 2n+10m.

The choice of threshold c is crucial to the 4S algorithm, and may need to be determined on a case-by-case basis. Here we offer a general guideline for parameter selection. Recall that for the Gaussian model, the signal strength of a segment with length L and height ν=δσ is usually measured by S=δ^2L; see, e.g. Table 1 in <cit.>. If there are two segments with the same overall signal strength S, however, one with large δ and small L (say, type A), and another one with small δ and large L (say, type B), then it is usually not equally easy to detect both of them by an algorithm of complex O(n). Indeed, for many segment detection algorithms, it is tricky to balance the powers to detect these two types of segments. Intuitively, the threshold parameter c controls this tradeoff in our methods. A higher threshold may be more powerful in detecting type A segments but less powerful in detecting type B segments; and vice versa. In practice, we may choose the threshold as a certain sample percentile of the absolute values of the observations based on a pre-specified preference. For example, if we know the signal segments have relatively large height parameters but can be as short as 5 data points, then with a fixed n, we can find largest m such as P̂(_n,m,5,5)≤ 0.05 and the percentile is chosen as α=1-m/n. That is, we want to guarantee that a segment of 5 consecutive black balls is significant enough to stand out. In another scenario, our preference might be longer segments with possibly lower heights. Then we may choose a threshold to include segments of length 10 with at least 6 black balls. In general, such m (or α) can be easily determined given n, s, t and p, by solving P̂(_n,m,s,t)≤ p. We illustrate in Figure <ref> the relationship between log n and selected percentile α for (s,t)=(5,5), (10,6), p=0.05 and 0.1. Because our main goal in this paper is to identify short segments, we prefer a large threshold such as 95th sample percentile. An even larger threshold can be used to identify shorter segments, and a smaller threshold can be used for detecting longer segments with lower heights.





§ NUMERICAL STUDIES



 §.§ Simulated data

We use simulation studies to evaluate the performance of our method in terms of the average number of true positives (TP) and false positives (FP) for identifying signal segments. Recall that in our definition, a detected segment Î∈ is a true positive, if it interacts with only one signal segment I∈, and it is the only one in  that interacts I.

In Example 1, we show the effectiveness of our inference framework on the false positive control of the 4S algorithm by a null model. As suggested by Figure <ref>, we choose the 95th percentile of absolute values of the observations as the threshold c. We set d=9 and h=3. We use various p-value thresholds for false positive control and compare them with a vanilla version of 4S, which is the one without p-value control.
[Null Model]
We generate a sequence based on model (<ref>) with n=10,000 and =. We consider three scenarios for the error distributions. In the first two scenarios, we consider {_̌i}_i=1^n which are IID from N(0,1) and t_3, respectively. In the last scenario, we consider {_̌i}_i=1^n which are marginally N(0,1) and jointly from an autoregressive (AR) model with autocorrelation 0.2.


As =∅ in this example, all detected segments are FPs. We report the average FPs for three versions of 4S (Vanilla, p=0.05 and p=0.1) based on 100 replicates in Table <ref>. We see that our inference framework can effectively control the number of FPs.



In Example 2, we compare 4S with three algorithms CBS <cit.>, LRS <cit.> and WBS <cit.>. The CBS and WBS methods, implemented by  packages  and  respectively, give a segmentation of the sequence which consists of a set of all segments rather than only the signal segments. In order to include their results for comparison, we ignore the long segments (with length greater than 100) detected by CBS or WBS, which decreases their false positives. For LRS, we set the maximum length of signal segments as 50.


We generate a sequence based on model (<ref>) with n=10,000. There are 5 signal segments with lengths 8, 16, 24, 32 and 40 respectively. We use the same error distributions as in Example 1. We consider two levels of height parameter for different signal strengths. In particular, we set height ν as the 99-, and 97-th percentiles of the marginal error distribution in two scenarios, labeled by S1 and S2. For the standard normal error, the height values are 2.326 and 1.881, respectively.


The threshold c we used for 4S is the 95-th sample percentile of absolute values, that is around the 97.5-th percentile of the error distribution, e.g., around 1.96 for the Gaussian case. Therefore, the true height is greater than c in S1, but lower in S2. Average numbers of TPs and FPs are reported in Tables <ref> and <ref>.












Overall, CBS performs the best for the IID Gaussian case, but suffers from a low power in the heavy-tail case, and high FPs in the correlated case. LRS and WBS perform reasonable well with slightly high FPs in the heavy-tail case. 4S methods are more robust against the error type. When the noise is Gaussian and the signal strength is weak, it is slightly less powerful than the methods based on the Gaussian assumption. In terms of computation time (Table <ref>), 4S is about 100 times faster than all other methods.





 §.§ Real data example


We applied the 4S method to the 272 individuals from HapMap project. In particular, we tried 4S (with p-value thresholds 0.05 and 0.5) to the LRR sequence of chromosome 1, which consists of 33991 measurements for each subject. We compared 4S with CBS, which has been a benchmark method in CNV detection. Note that CBS produces a segmentation of the sequence rather than the CNV segments directly. Therefore, we focused on only the short (less than 100 data points) segments detected by CBS because the long segments had means close to zero and are not likely to be CNVs. We found that most of these short segments are separated. But a very small portion of them are connected as CBS sometimes tends to over segment the sequence. Therefore, we merged two short segments detected by CBS if they are next to each other.

The 4S algorithm is extremely fast. It took less than 2 seconds (on a desktop with CPU 3.6 GHz Intel Core i7 and 16GB memory) to complete 272 sequences with p-values calculated for all detected segments. CBS algorithm is reasonably fast, but much slower than our algorithm. In Table <ref> we list the total number of detected CNVs, average length of CNVs and computation time for all methods.



Overall, the segment detection results were very similar. We further compared the segments detected by two algorithms, i.e., 4S with threshold p=0.05 and CBS. We found that 2753 segments are in common. Here by a common segment we mean a pair of segments, one detected from each algorithm, such that they overlap to each other but do not overlap with other detected segments. Among these common segments, we calculated a similarity measure, called affinity in <cit.>, defined as follows.

    ρ(I,I')=|I∩ I'|/√(|I|·|I'|).

ρ(I,I')=1 if two segments I and I' are the same and ρ(I,I')=0 if they do not overlap. We found that the average value of this similarity measure is 0.9290 among 2753 pairs.

Figure <ref> presents the histogram of affinity among 2753 pairs of commonly detected segments by 4S and CBS. We can see that 87.76% of those pairs have affinity values larger than 0.8. We further divided the detected segments into three groups: those detected by both methods (group 1); those detected by only 4S method with p=0.05 (group 2); those detected by only CBS (group 3). For each detected segment, we calculated its length and the sample mean of the measurements on the segment. Figure <ref> displays the scatter plots of sample means versus lengths for all the segments in three groups. Most segments in group 1 carries relatively strong signals. So it is not surprised that they were detected by both algorithms. The groups 2 and 3 have much smaller sizes than group 1. In particular, we found that most segments in group 3 (i.e. those detected by only CBS) are very short, consisting of only 2 or 3 data points. Those segments are not significant in our inference framework unless we set a very high threshold c in step 1. Some segments in group 2 (i.e. those detected by only the 4S method) have relatively small sample mean values, which explains why they were not detected by CBS. Some of these segments might be true positives with the sample mean affected by outliers. Overall, the 4S and CBS methods gave similar results. The segments detected by only one method may be prone to false positives, or true positives have weak signal strengths.








§ DISCUSSION

We proposed a scalable nonparametric algorithm for segment detection, and applied it to real data for CNV detection. Two main advantages of the 4S algorithm are its computational efficiency and independence of the normal error assumption. We introduced an inference framework to assign significance levels to all detected segments. Our numerical studies demonstrated that our algorithm was much faster than CBS and performed similarly to CBS under the normality assumption and better when the normality assumption was violated. Although our inference framework depended on the assumption of IID noise, our numerical experiments suggested that our algorithm worked well under weakly correlated noises. Hence, the proposed method is faster and more robust against non-normal noises than CBS. Overall, the 4S algorithm is a safe and fast alternative to CBS, which has been a benchmark method in CNV studies.

In the literature, there are two popular classes of change-point models used to study CNV related problems. The first one assumes only a piecewise constant median/mean structure. The second one assumes, in addition, a baseline, which reflects the background information or normal status of the data. Quite often, it assumes that the abnormal part, called signal segments in our paper, are sparse. For the first approach, the goal is to identify the change points. In contrast, the second approach emphasizes more on segment detection rather than change-point detection. The difference is subtle for estimation but might become remarkable for inference. For example, it is technically difficult to define `true positive' in the context of change-point detection <cit.>. But it is easier to define related concepts for segment detection as we did in this paper. Roughly speaking, the first approach is more general, and the second one is more specific and suitable to model certain CNV data, e.g., SNP array data. In particular, the 4S algorithm aims to solve change-point models in the second class. It can be applied to any data sequence when there is a baseline. When the baseline mean/median is unknown, we suggest that the data should be centered first by the estimated mean/median. Our method can not be applied to data when a baseline does not exist. Besides change-point models, there are other approaches to study CNV such as hidden Markov model <cit.>. Due to the space limit, we restricted our comparison to the methods based on change-point models and implemented by R packages.


Most segment detection algorithms involve one or more tuning parameters, whose values are critical to the results. In the study of segment detection, there are two trade-offs that researchers should consider in choosing algorithms as well as their parameters. The first one is the usual type I/type II errors trade-off, which might be tricky sometimes but well-known. The second one is more delicate and quite unique. For a signal segment, both its height and length determine the signal strength. Therefore, segments with weak but detectable signals can be roughly divided into two categories, the ones with small length (say, type A) and the ones with small height (say, type B).
Typically a method may detect type A segments more powerfully, but type B segments less powerfully, than the other method. For the proposed 4S algorithm, a choice of a larger threshold parameter in step 1 makes the algorithm more powerful in detecting short and high signal segments (type A), and vice versa. The 4S algorithm can be easily tuned to maximize the power in detecting of a certain type of the signal segments. We may also try different thresholding levels in data analysis in order to detect different types of segments. In general, the choice of the parameters depends on the research goals and balance of two trade-offs mentioned above.

There are various platforms and technologies which produce data for CNV detection. Besides the SNP array data studied in this work, read depth data from next generation sequencing (NGS) technologies are often used in CNV studies. As one referee pointed out, the speed of 4S algorithm would be an advantage when applied to read depth data from whole genome sequencing. This is a wonderful research direction that we will investigate next.




An R package  implementing our proposed method can be download via

.





§ ACKNOWLEDGEMENT
 The authors are partially supported by National Science Foundation, National Institutes of Health, Simons Foundation and University of Arizona faculty seed grant.



§ APPENDIX


Proof of Lemma 1.

Let I be the interval of integers [ℓ,r] with L=r+1-ℓ. For each X_i, i∈ I, the probability that the ball at i is white is π=P(|X_i|≤ c)≤ P(X_i≤ c)≤1/2 as ν≥ c and f is symmetric. It is trivial to bound P() for the case d≥ L as P(^c)=π^L≤ 2^-L.  Now let us consider the case d<L. Let _i, i∈ I be the event that the first segment of d consecutive white balls starts from position i. Then

    P(_i){[         =π^d,      if  i=ℓ;;   ≤ (1-π)π^d, if ℓ<i≤ r+1-d;           =0, if  i≥ r+2-d. ].

Therefore, P(^c)= ∑_i∈ I P(_i) ≤π^d+ (L-d)(1-π)π^d ≤ (1+1/2(L-d))(1/2)^d and P()≥ 1- (L-d+2)2^-d-1. Note that π is a constant depending on f, ν and c, so a sharper bound than 1/2 for π may be used to bound P() if more information is available.

Let us consider the segment [r+1,r+D] on the right hand side of I. _D^c implies that there is at least one black ball in each of the segments [r+1,r+d], [r+d+1,r+2d], etc. Note that X_i∼ F on these segments so the probability of white ball at i is P(|X_i|≤ c=f_β)=β'=2β-1. Consider all ⌊D/d⌋ segments of length d on the right side of I. The probability that all these segments contain at least one black ball is (1-β'^d)^⌊D/d⌋. Therefore, P(_D)≥ 1-2(1-β'^d)^⌊D/d⌋. □

Proof of Theorem 1.

For =∪_k=1^K I_k, let L_k=|I_k| and D_k be the gap between I_k and I_k+1. Define

    _k    =   {on I_k, there does not exist a sub-segment of min{d,L_k} consecutive white balls}
    _k    =   {there are  d  consecutive white balls on  D_k }.

Note that all segments in  are identified under event (⋂_k=1^K_k)⋂(⋂_k=1^K-1_k). By Lemma <ref>, P(_k^c)≤ 2^-L_k or 1/2(L_k-d+2)2^-d, which can be bounded by max{1/2(L_max-d+2),1}2^-min{d,L_min}. Moreover, P(_k^c)≤ (1-β_min'^d)^⌊D/d⌋≤ (1-β_min'^d)^⌊D_min/d⌋. The conclusion follows Bonferroni inequality. □

Proof of Theorem 2.

Let _c={j_1,...,j_m} be the locations of black balls after step 1. Note that j_i and j_i+1 will be connected in step 2 if and only if j_i+1-j_i≤ d. We aims to count the number of segments with at least 2 consecutive black balls after step 2, as all isolated black balls will be eliminated in step 3. Such a segment starts at j_i only if j_i+1-j_i≤ d. So the total number of such segments is at most |{i: j_i+1-j_i≤ d}|. Let Z_i follow Bernoulli distribution with Z_i=1 if and only if j_i+1-j_i≤ d for i=1,...,m. When =, all black balls are randomly distributed. P(Z_i=0), i.e., the probability that all balls are white in next d positions following j_i is n-1-d m-1/n-1 m-1=∏_k=1^dn-m+1-k/n-k. Therefore, E|_c,d,h|≤∑_i=1^m Z_i≤ m (1-∏_k=1^dn-m+1-k/n-k). □


Proof of Lemma 2.

We drop the subscript k in P(_n,m,s_k,t_k) as it is irrelevant in our derivation below. Under the assumption that m black balls are randomly assigned in n position, at a position of black ball, we calculate the probability that there are at least t-1 black balls in next s-1 positions. Let Y be the count of black balls in those s-1 positions. Y follows a hypergeometric distribution with total population size n-1, number of success states m-1, and number of draws s-1. Therefore, P(_n,m,s_k,t_k)≤ mP(Y≥ t-1) as there are m black balls. □


biometrika



