Detection of an anomalous cluster in a network

Introduction

We discuss the problem of detecting whether or not, in a given network, there is a cluster of nodes which exhibit an “unusual behavior.” Suppose that we are given a set of nodes with a random variable attached to each node. We observe a realization of this process and would like to tell whether all the variables at the nodes have the same behavior, in the sense that they are all sampled from a common distribution, or whether there is a cluster of nodes at which the variables have a different distribution.

The task of detection in networks is critical for an increasing number of applications, for example, in surveillance and environment monitoring. We describe a few of these applications below.

The advent of sensor networks overviewsensor , 1024422 , zhao2004wireless has multiplied the amount of data and the variety of applications where the task of detection is central. Surveillance and environment monitoring are prime areas of application for sensor networks. Take, for example, the transport of hazardous materials. Currently, some major traffic bottlenecks (e.g., airports, subways and borders) use portal monitoring systems fitch2003tcr , 1352095 . Sensor networks offer a more flexible, decentralized alternative and are considered for the detection of radioactive, biological or chemical materials hills2001sd , radiation , cui2001nnh . Sensor networks are also extensively used in other target tracking settings linesand , targetsensor .

Detection in digital signals and images

A digital camera may be seen as a sensor network, with CCD or CMOS pixel sensors. As imaging systems have been available for quite some time, the literature on detection in images is quite extensive, spanning several decades, particularly in satellite imagery roadtracking , artificialnatural , shipdetection , firedetection , computer vision eyedetection , objecttracking and medical imaging medicalsurvey , breasttumor , braintumor , braams1995detection .

Disease outbreak detection

The presence of a biological or chemical material in a given geographical region may also be detected indirectly through its impact on human health. In this context, early detection of the disease outbreak is crucial in order to minimize the severity of the epidemic. For that purpose, some specific information networks are used, with surveillance systems now incorporating data from hospital emergency visits, ambulance dispatch calls and pharmacy sales of over-the-counter drugs rotz2004advances , heffernan2004ssp , wagner2001esv .

Virus detection in a computer network

Diseases affect computers as well, in the form of viruses and worms spreading from host to host in a computer network szor2005acv . Affected machines may exhibit slightly anomalous behavior (e.g., a loss of performance or violations of specific rules) which may be hard to detect on an individual machine.

Detection from field measurements

In waterquality , the water quality in a network of streams in Pennsylvania is assessed by field biologists performing a variety of analyses at various locations along the streams; the objective is to determine whether there are regions of low biological integrity based on the collected data, and to identify these regions. Other field measurements include census data and surveys involving geographical location.

Detection is, of course, closely related to estimation (i.e., the localization or extraction of the anomalous cluster of nodes), but different. This distinction is rarely made clear, however. Indeed, reliable detection is possible at lower signal-to-noise ratios than reliable estimation and it may be important to detect the presence of signals from noisy data without being able to estimate them. For example, one could imagine developing a surveillance system performing detection at relatively low energy/bandwidth costs, yet efficient at low signal-to-noise ratios, and then switching to estimation mode whenever the presence of a signal is detected. Another example would be a low cost preliminary survey involving fewer field measurements, with findings subsequently confirmed by a larger, more expensive survey.

2 Mathematical framework

The random variables are assumed to be independent. For concreteness, we consider a normal location model, popular in signal and image processing, to model the noise. Our analysis, however, generalizes to any exponential family under some conditions on the sizes of the anomalous clusters, such as Bernoulli models which arise in sensor arrays where each sensor collects one bit (i.e., makes a binary decision) or Poisson models which come up with count data, for instance, arising in infectious disease surveillance systems Kul . The extension to exponential families is detailed in Section 4.1.

The situation where no signal is present, that is, “business as usual,” is modeled as

where $\mu_{K}>0$. We choose to decompose $\mu_{K}$ as $\mu_{K}=|K|^{-1/2}\Lambda_{K}$, where $|K|$ denotes the number of nodes in $K$ and $\Lambda_{K}$ is the signal strength. Indeed, with this normalization, for any cluster $K$,

Figures 1–4 illustrate the setting for various types of clusters.

The situation we just described is purely spatial and relevant in some applications not involving time. Such situations are common in image processing. In other applications, especially in surveillance, time is an intrinsic part of the setting. In the following section, we modify the model above to incorporate time.

2.2 Spatio-temporal model

3 Structured multiple hypothesis testing

In stark contrast, in all of the applications described earlier, the set of alternatives is composite. Viewing each node as performing a test of hypotheses, which is common in the literature on sensor networks, our problem falls within the framework of multiple comparisons. Multiple hypothesis testing is a rich and active line of research which is receiving a considerable amount of attention within the statistical community at the moment; see MR2387976 and references therein. The vast majority of the papers assume that the tests are independent of each other, which is clearly not the case here since, in general, the class contains clusters that intersect. This is particularly true in engineering applications, although this assumption is often made MR2213561 , MR926001 , MR1951265 .

4 The scan statistic

We will focus on the test that rejects for large values of the following version of the scan statistic:

The chosen normalization is such that each term in the maximization is standard normal under the null and allows us to compare clusters of different sizes. It corresponds to the generalized likelihood ratio test in our context if $\Lambda_{K}$ is independent of $K\in\mathcal{K}_{m}$. The scan statistic was originally proposed in the context of cluster detection in point clouds glaz01 . This is the method of matched filters which is ubiquitous in problems of detection in a wide variety of fields, sometimes in the form of deformable templates in the engineering literature deformablereview , medicalsurvey or their nonparametric equivalent, active contours or snakes xu1998ssa . Note that the scan statistic is the prevalent method in disease outbreak detection, with many variations kulldorff2005stp , treebased , kulldorff2006ess , duczmal2006ess .

and will restrict the scanning to an $\varepsilon$-net of $\mathcal{K}_{m}$ with respect to $\delta$, that is, a subset $\{K_{j}\dvtx j\in J\}\subset\mathcal{K}_{m}$ with the property that for each $K\in\mathcal{K}_{m}$, there is a $j\in J$ such that $\delta(K,K_{j})\leq\varepsilon$. We will elaborate on this approach in the LABEL:supp. When $J$ is minimal, we call the resulting statistic an $\varepsilon$-scan statistic. The approximation precision $\varepsilon$ will be chosen appropriately, depending on the situation.

We focus on $\varepsilon$-scan statistics for two reasons. First, their performance is easier to analyze than that of the scan statistic itself; in fact, the main approach to analyzing the scan statistic, the chaining method of Dudley MR0220340 , MR2133757 , is via a properly chosen $\varepsilon$-scan statistic. Second, some of the classes we consider are rather large and we believe that it would be computationally impractical to scan through all of the clusters in the class; furthermore, our results show that, from an asymptotic standpoint, no substantial improvement would be gained by using the full scan statistic.

We also note that the tuning parameter $\varepsilon$ may be dispensed with if we scan over subsets of different sizes in a multiscale fashion and use a scale-dependent threshold.

5 Existing theoretical results

Most of the literature assumes that the class $\mathcal{A}$ is parametric, exemplified by deformable templates, for which theoretical results are available, especially in the case of the square lattice boxscan , MGD , morel , jiang , perone , boutsikas . In particular, with a normal location model, the scan statistic performs well, in the sense that it is asymptotically minimax; this is shown in MGD in a slightly different context tailored to image processing applications. We also mention the recent work HJ09 , which considers the detection of multiple clusters (intervals) of various amplitudes in the one-dimensional lattice. As for nonparametric classes of domains, MGD argues that the scan statistic is asymptotically minimax for the case of star-shaped clusters with smooth boundaries.

6 New theoretical results

We describe here in an informal way the results we obtain.

Note that the detection rate is the same as for the class of balls so that, perhaps surprisingly, scanning for the location (and not the shape) is what drives the minimax detection risk.

Hence, some form of scan statistics achieves a detection rate within a factor of $(\log m)^{3/2}$ from the minimax rate.

In Section 2.3, we consider the spatio-temporal setting. We first consider cluster sequences that admit a “thick” limit. Cellular automata, which have been used to model epidemics MR1683061 , satisfy this condition in some cases. In Proposition 5, we show that scanning over space–time cylinders, as done in disease outbreak detection, achieves the asymptotic minimax risk. We then consider cluster sequences with controlled space–time variations, which may be a relevant model for applications such as target tracking targetsensor . We consider a fairly general model in Proposition 7.

Therefore, some form of scan statistic is again within a factor of $(\log m)^{3/2}$ from optimal. In Section 3.2, we consider arbitrary connected components, constraining only the size; see Figure 4. In Proposition 9, we obtain a sharp detection rate for clusters of very small size.

7 Structure of the paper

We have just described the contents of Sections 2 and 3. Section 4 is our discussion section. We extend the results obtained for the normal location model to any exponential family in Section 4.1. Other extensions are described in Section 4.2. We state some open problems in Section 4.5. In Section 4.4, we briefly discuss the challenge of computing the scan statistic. The technical arguments are gathered in the LABEL:supp.

8 Notation

Clusters as geometric shapes in Euclidean space

In this section, we consider clusters as in (4), where $\mathcal{A}$ is a class of bi-Lipschitz deformations of the unit $d$-dimensional ball. This includes the vast majority of all the parametric clusters considered in the literature, such as hyperrectangles and ellipsoids, as long as the shape is not too narrow. Note that a slightly less general situation is briefly mentioned in MGD .

We start with a lower bound on the minimax detection rate for discrete balls of a given radius.

Consider $\lambda_{m}\to 0$ such that $\lambda_{m}\geq r_{m}^{*}$ and let $\mathcal{K}_{m}$ be the class of all discrete balls of radius $\lambda_{m}$, that is,

Note that $\lambda_{f}$ is intimately related to the size of $\operatorname{im}(f)$ and therefore of $K_{f}$. Indeed, a simple application of (6) implies that, for any $f\in\mathcal{F}_{d,d}(\kappa)$,

This implies that sets of the form $\operatorname{im}(f)$, with $f\in\mathcal{F}_{d,d}(\kappa)$, are “thick,” in the sense that the smallest ball(s) containing $\operatorname{im}(f)$ and the largest ball(s) included in $\operatorname{im}(f)$ are of comparable sizes.

Consider $\lambda_{m}\to 0$ such that $\lambda_{m}\geq r_{m}^{*}$ and define

where $\eta_{m}=\varepsilon_{m}^{2}\sqrt{2d\log(1/\lambda_{m})}$. Moreover, if $r_{m}^{*}\asymp m^{-1/d}$ and

where $\eta_{m}=\varepsilon_{m}^{2}\sqrt{2\log m}$.

Therefore, on a larger class of mild deformations of the unit ball, some form of scan statistic achieves essentially the same detection rate as for the class of balls stated in Proposition 1.

We note that the lower bound on $\underline{\Lambda}_{m}$ is driven by the smaller clusters in the class and that the performance guarantee is subject to a proper choice of $\varepsilon_{m}$. A simple fix for both issues is to combine the tests for different cluster sizes with an appropriate correction for multiple testing. We summarize the consequence of Proposition 1 and Theorem 1 with this observation in the following result, inspired by boxscan .

Consider $\lambda_{m}\to 0$ such that $\lambda_{m}\geq r_{m}^{*}$ and define

The same procedure, that is, combining $\varepsilon$-scan statistics at different (dyadic) scales, may be implemented in any of the settings we consider in this paper to obtain a test that does not depend on a tuning parameter like $\varepsilon_{m}$ and achieves the same optimal rate at every size. This is simply due to the fact that we only need to consider the order of $\log m$ scales and the fast decaying tails of the scan statistics under the null.

In a number of situations, the signal to be detected may be composed of several clusters. Our results extend readily to this case. Let $j_{m}$ be a positive integer and consider sets of the form $\bigcup_{j=1}^{j_{m}}K_{f_{j}}$, where the union is over some $f_{j}\in\mathcal{F}_{d,d}(\kappa)$ such that, for $j,j^{\prime}$, $\lambda_{f_{j}}\leq C\lambda_{f_{j^{\prime}}}$ and

so that the sets $\operatorname{im}(f_{j})$ and $\operatorname{im}(f_{j^{\prime}})$ are of comparable sizes and not too far from each other. In that case, Theorem 1 applies unchanged, as long as the number of clusters is not too large, specifically if $j_{m}=o(\log(1/\lambda_{m}))^{1/d}$. (This can be improved if the $K_{f_{j}}$’s do not overlap too much.) If the proximity constraint is dropped, then the term $\log(1/\lambda_{m})$ in Theorem 1 is replaced by $j_{m}\log(1/\lambda_{m})$.

2 Thin clusters

In this section, we consider clusters that are built from smooth embeddings in $\Omega_{d}$ of the unit $p$-dimensional ball, where $p<d$. The special case of curves ($p=1$) is, for example, relevant in road tracking roadtracking and in modeling blood vessels in medical imaging bloodvessel04 . As in the previous section, the results we obtain below are valid for (some) unions of such subsets and, in particular, for submanifolds with a wide array of topologies.

Again, $\lambda_{f}$ is intimately related to the size of $B(\operatorname{im}(f),r)$ and $K_{f,r}$. This relationship is made explicit in the LABEL:supp. We consider classes of clusters of the form $\{K_{f,r}\dvtx f\in\mathcal{F}\}$, where $\mathcal{F}$ is a subclass of $\mathcal{F}_{p,d}(\kappa)$.

Let $C$ be the constant defined in Lemma B.2 in the \setattributereffmtLABEL:supp\setattributereffmt. Consider $\lambda_{m},r_{m}\to 0$ such that $C^{-1}\lambda_{m}\geq r_{m}\geq r_{m}^{*}$ and let $\mathcal{F}$ be a subclass of $\mathcal{F}_{p,d}(\kappa)$. Define

Just as in Theorem 1, if $r_{m}^{*}\asymp m^{-1/d}$, we can dispense with the restriction $r_{m}\geq r_{m}^{*}$ and replace the factor $\log(1/\lambda_{m}^{d})$ by $\log m$ in the bound.

For a typical parametric class $\mathcal{F}$, $\log N_{\varepsilon}(\mathcal{F})\sim a(\mathcal{F})\log(1/\varepsilon)$, so the scan statistic (over an appropriate net) is accurate if

On the other hand, $\log N_{\varepsilon}(\mathcal{F})\asymp(1/\varepsilon)^{a(\mathcal{F})}$ for a typical nonparametric class $\mathcal{F}$ MR0124720 , so the scan statistic (over an appropriate net) is accurate if

Finding a sharp lower bound for the minimax detection rate is more challenging for thin clusters compared to thick clusters. By considering disjoint tubes around $p$-dimensional hyperrectangles, we obtain a lower bound that matches, in order of magnitude, the rate achieved by the scan statistic when the class $\mathcal{F}$ is parametric, displayed in (8).

The proof is parallel to that of Proposition 1 and is therefore omitted.

For at least one family of nonparametric curves ($p=1$), we show that the rate displayed at (9) matches the minimax rate, except for a logarithmic factor. For concreteness, we assume that $\Omega_{d}=^{d}$. Let $\mathcal{H}(\alpha,\kappa)$ be the Hölder class of functions $g\dvtx\to$ satisfying

Let $r_{m}\to 0$ with $r_{m}\geq r_{m}^{*}$. Let $\mathcal{F}$ be the class of functions of the form $f(x)=(x,g_{1}(x),\ldots,g_{d-1}(x))$, where $g_{j}\in\mathcal{H}(\alpha,\kappa)$, with $\alpha\geq 2$. Define

Thus, for the detection of curves with Hölder regularity, a scan statistic achieves the minimax rate within a poly-logarithmic factor. We prove Proposition 3 by reducing the problem of detecting a band in a graph so that we can use results from Section 3.1. We do not know how to generalize this approach to higher-dimensional surfaces (i.e., $p\geq 2$).

3 The spatio-temporal setting

In this section, we consider the spatio-temporal setting described in Section 1.2.2. This is a special case of the spatial setting we have considered thus far, with time playing the role of an additional dimension. For their relevance in applications and concreteness of exposition, we focus on two specific models. In Section 2.3.1, we consider cluster sequences with a limit; as we shall see, this assumption is implicit in some popular models for epidemics. In Section 2.3.2, we consider cluster sequences of bounded variations.

Regarding the actual evolution of the cluster in time, a number of growth models have been suggested, for example, cellular automata MR2367301 , MR1849342 and their random equivalent, threshold growth automata MR2206346 , MR1698409 , which have been used to model epidemics MR1683061 . The latter includes the well-known Richardson model Richardson1973jk . Under some conditions, these models develop an asymptotic shape (with probability one), a convex polygon in the case of threshold growth automata. Less relevant for modeling epidemics, internal diffusion limited aggregation is another growth model with a limiting shape Lawler1992jk .

The simplest cluster sequences with limiting shape are space–time cylinders, for which we have the equivalent of Proposition 1. (The proof is completely parallel and we omit details.)

If the starting time is uniformly bounded away from $t_{m}$ and the convergence to the thick spatial cluster [in the sense of (10)] occurs at a uniform speed, then all of the cluster sequences in the class have sufficient time to develop into their “limiting” shapes. The space–time cylinders over which we scan are based on an $\varepsilon$-net for the possible limiting shapes, that is, the class of thick clusters.

Scanning over space–time cylinders (with balls as bases) is advocated in the disease outbreak detection literature diseaseoutbreak . Although seemingly naive, this approach achieves, in our asymptotic setting, the minimax detection rate if the cluster sequences develop into balls and, in general, falls short by a constant factor.

We mention that the equivalent of Corollary 1 holds here as well, in that we can combine the different scans at different space–time scales to obtain a test that does not depend on a tuning parameter (implicit here) and which achieves the same rate for the cluster class defined as above, but with $\lambda_{m}\geq\lambda_{f}\geq r_{m}^{*}$, which is the class that appears in Corollary 1.

3.2 Cluster sequences of bounded variation

In target tracking linesand , targetsensor , the target is usually assumed to be limited in its movements due to maximum speed and maneuverability. With this example in mind, we consider classes of cluster sequences of bounded variation, meaning that the cluster is limited in the amount it can change in a given period of time. As the rates we obtain in this subsection are the same with or without the condition $K_{t_{m}}\neq\varnothing$, we do not make that assumption. Let $t_{K}^{+}=\max\{t\dvtx K_{t}\neq\varnothing\}$.

We consider space–time tubes around Hölder space–time curves. For $\alpha\in(0,1]$ and $\kappa>0$, let $\mathcal{H}_{\infty}(\alpha,\kappa)$ be the Hölder class of functions $g\dvtx[0,\infty)\to$ satisfying

The following is the equivalent of Proposition 3.

For simplicity, assume that $t_{m}$ is a power of $m$. If $\xi_{m}r_{m}^{1/\alpha}=O(1)$, then the detection threshold is roughly of order $t_{m}^{1/2}$, while if $\xi_{m}r_{m}^{1/\alpha}$ is large, yet small enough that $t_{m}/(\xi_{m}r_{m}^{1/\alpha})$ is still a power of $m$, then the detection threshold is roughly of order $(t_{m}/(\xi_{m}r_{m}^{1/\alpha}))^{1/2}$.

A form of scan statistic is actually able to attain the same detection rate when the radius is unknown, but restricted to $r\geq r_{m}$. In fact, another form of scan statistic achieves a slightly different rate over a much larger class of cluster sequences with bounded variations. Let $\mathcal{S}(r,\kappa)$ be the set of subsets $S\subset\Omega_{d}$ such that $B(x,r)\subset S\subset B(x,\kappa r)$ for some $x\in\Omega_{d}$.

for a function $\nu\dvtx[0,\infty)\to[0,\sqrt{2}]$. Then, (12) is satisfied with $\eta=\nu(1)$ and the same $\xi_{m}$. The requirement in Proposition 7 is that $\nu(1)$ be small enough. In particular, the cluster sequences considered in Proposition 6 satisfy, for some constant $C>0$,

This comes from Lemma C.1 in the LABEL:supp and (11). Therefore, assuming $\xi_{m}\ll(r_{m}^{*})^{-1/\alpha}$, (13) is satisfied with $\nu(u)=u^{\alpha/2}$ and $\xi_{m}$ replaced by $\xi_{m}r_{m}^{1/\alpha}$. In that case, the detection rates obtained by the scan statistics of Propositions 6 and 7 are of comparable orders of magnitude.

Clusters as connected components in a graph

coming slightly closer to the minimax rate.

In fact, even when the band has unknown length, width and starting location, and when the path is not restricted to be nondecreasing, a form of scan statistic achieves the same rate, except for a logarithmic factor.

2 Arbitrary connected components

We consider here classes of connected components with a constraint on their sizes. Arbitrary connected components in the square lattice are sometimes called animals or polyominoes (polycubes in dimension $d\geq 3$), which are well-studied objects in combinatorics, where the goal is to count the number of polyominoes Klarner1967fp . We mention in passing the results in MR1825148 which provide a law of large numbers for the scan statistic under the null. Otherwise, such objects are fairly new to statistics. Detecting animals is, of course, harder than detecting paths since paths are themselves animals. The result below offers a sharp detection threshold for connected components of sufficiently small size.

Discussion

Under the alternative, the variables at the nodes belonging to the anomalous cluster $K\in\mathcal{K}_{m}$ have distribution $F_{\theta_{K}}$ with $\theta_{K}:=\sigma\Lambda_{K}|K|^{-1/2}$, that is,

As before, the variables are assumed to be independent.

If the clusters in the class are sufficient large, then the results presented for the normal location family hold unchanged. Intuitively, large enough clusters allow for the sums over them to be approximately normally distributed. Details are provided in the LABEL:supp. For example, we have the following equivalent of Corollary 1 in the context of thick clusters as in Section 2.1. Consider $\lambda_{m}\geq r_{m}\geq r_{m}^{*}$ with $mr_{m}^{d}(\log 1/r_{m})^{-3}\to\infty$ (which guarantees that the clusters in the class are large enough) and define the class

In this setting, under the Bernoulli model, the detection threshold is at

under the Poisson model, the detection threshold is at

Note that without a lower bound on the minimum size of the anomalous clusters, the general analysis breaks down and the results depend on the specific exponential model. For example, unless $\min\{|K|\dvtx K\in\mathcal{K}_{m}\}\to\infty$ fast enough, detection is impossible in the Bernoulli model, even if the anomalous nodes have value 1 under the alternative.

2 Other extensions

The array of possible models is as wide as the breadth of real-world applications. We mention a few possible variations below.

Using an exponential family of distributions allows us to obtain sharp detection lower bounds. Otherwise, similar results, although not as sharp, may be obtained for essentially any family of distribution $F_{\theta}$, where the distance between the null $\theta=0$ and an alternative $\theta$ is in terms of the chi-square distance between $F_{0}$ and $F_{\theta}$; see maze , Section 5.

Different means at the nodes

We could consider a situation where the mean varies over the nodes of the anomalous cluster. This situation is considered in HJ09 for the case of intervals, and the constant in the detection rate is indeed different. We implicitly considered a worst case scenario where the mean is bounded below over the anomalous cluster and subsequently assumed it was equal to that lower bound everywhere over the anomalous cluster. However, our results hold unchanged if we allow $X_{v}$ to have any mean above $\theta_{K}$, for every $v\in K$, $K$ being the anomalous cluster.

Dependencies

Also of interest is the case where the variables are dependent. In the spatial setting, the same paper HJ09 solves this problem for the case of the one-dimensional lattice, with the correlation between $X_{v}$ and $X_{w}$ decaying as a function of distance between $v$ and $w$. We postulate that the same result holds in higher dimensions. In the spatio-temporal setting, variables could be dependent across time as well, involving a higher degree of sophistication. We plan on pursuing these generalizations in future publications.

Unknown variance or other parameters

We assumed throughout that the variance was known (and equal to 1 after normalization). This is, in fact, a mild assumption, as one can consistently estimate the variance using a robust estimator, say the median absolute deviation (MAD), with the usual $\sqrt{m}$-convergence rate, assuming that the anomalous cluster corresponds to a small part of the entire network. When dealing with one-parameter families such as Bernoulli or Poisson, the issue is to estimate the parameter under the null and a robust version of the maximum likelihood (e.g., trimmed mean for these two examples) can be used for that purpose.

3 Energy, bandwidth and other constraints

We assume throughout that a central processor has access to all the information measured at the nodes and, based on that, makes a decision as to whether there is an anomalous cluster of nodes in the network or not. This assumption is reasonable in, for example, the context of image processing or syndromic surveillance. However, real-world sensor networks of the wireless type are often constrained by energy and/or bandwidth considerations. A growing body of literature targetsensor is dedicated to designing efficient (e.g., decentralized) communication protocols for sensor networks under such constraints. As mentioned in Section 1.3, the papers we are aware of consider very simplistic detection settings. In the context of the present paper, it would be interesting to study how the detection rates change when different communication protocols are used.

We also assume that we have infinite computational power. However, all real-world systems operate under finite energy and processing resources. In the same way, it would be interesting to know what detection rates are achievable under such computational constraints.

4 On computing the scan statistic

In all of the settings we consider in this paper, the scan statistic comes close to achieving the minimax detection rate. Turning to computational issues, however, it is very demanding, even when scanning for simple parametric clusters such as rectangles. For general shapes, Duczmal, Kulldorff and Huang duczmal2006ess suggests a simulated annealing algorithm, which, from a theoretical point of view, is extremely difficult to analyze. For parametric shapes and blobs, Arias-Castro, Donoho and Huo MGD advocates the use of $\varepsilon_{m}$-scan statistics based on multiscale nets built out of unions of dyadic hypercubes; similar ideas appear in boxscan . Partial results suggest that this approach yields, in theory, a near-optimal algorithm for detecting the more general thick clusters considered in Section 2.1.

For the thin clusters of Section 2.2, or for the bands of Section 3.1, the situation is quite different. Take the latter. After pre-processing the data by performing a moving average with an appropriate radius, it remains to find the maximum over a restricted, yet exponentially large, set of paths. Without further restriction, this problem, known as the “bank robber problem” or “reward budget problem” DasGuptaHespanhaRiehlSontag06 , is NP-hard. Note that DasGupta et al. DasGuptaHespanhaRiehlSontag06 suggests a polynomial time approximation that deserves further investigation. The case of thin clusters is even harder. In the context of point clouds, Arias-Castro, Efros and Levi AriasCastro2009 introduces multiscale nets that could be adapted to the setting of a network. It remains to compute the scan statistic over this net, which seems particularly challenging for surfaces of dimension $p\geq 2$, which no longer correspond to paths. In the spatio-temporal setting of Section 2.3, dynamic programming ideas could be used, as done in MSDFS in the context of point clouds and in chirpletpursuit in the context of a harmonic analysis decomposition of chirps.

5 Open theoretical problems

The paper leaves two main theoretical problems unresolved. The first one concerns obtaining sharper bounds for the detection of thin clusters. This is in the context of Section 2.2. For parametric classes, the challenge is to match constants in the rate, while, for nonparametric classes, the challenge is to obtain sharper lower bounds, perhaps closer to what a scan statistic is shown to achieve in Theorem 1. We were only able to do the latter for curves; see Proposition 3.

The second one concerns comparing the detection rates for arbitrary connected components and for paths. At a given size, the thicker the band (relative to its length), the easier it is to detect it; see Theorem 3. It seems, therefore, that the most difficult connected components to detect are paths or unions of paths. But is this true? In other words, are the minimax detection rates for arbitrary connected components and paths of a similar order of magnitude?

Acknowledgments

The authors are grateful to the anonymous referees for suggesting an expansion of the discussion section, for encouraging them to obtain sharper bounds and for alerting them of the possibility of improving on the performance of the scan statistic by using a different threshold for each scale, which resulted in Corollary 1.

[id=supp] \snameSupplement \stitleTechnical Arguments \slink[doi]10.1214/10-AOS839SUPP \sdatatype.pdf \sfilenamecluster-suppl.pdf \sdescriptionIn the supplementary file clustersuppl , we prove the results stated here. It is divided into three sections. In the first section, we state and prove general lower bounds on the minimax rate and upper bounds on the detection rate achieved by an $\varepsilon$-scan statistic. We do this for the normal location model first and extend these results to a general one-parameter exponential family. In the second section, we gather a number of results on volumes and node counts. In the third and last section, we prove the main results.

References