The Annals of Statistics 2013, Vol. 41, No. 2, 670–692 DOI: 10.1214/13-AOS1094 Institute of Mathematical Statistics, 2013 SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION1 BY YAO XIE AND DAVID SIEGMUND Duke University and Stanford University We develop a mixture procedure to monitor parallel streams of data for a change-point that affects only a subset of them, without assuming a spatial structure relating the data streams to one another. Observations are assumed initially to be independent standard normal random variables. After a change- point the observations in a subset of the streams of data have nonzero mean values. The subset and the post-change means are unknown. The procedure we study uses stream specific generalized likelihood ratio statistics, which are combined to form an overall detection statistic in a mixture model that hypothesizes an assumed fraction p0 of affected data streams. An analytic expression is obtained for the average run length (ARL) when there is no change and is shown by simulations to be very accurate. Similarly, an ap- proximation for the expected detection delay (EDD) after a change-point is also obtained. Numerical examples are given to compare the suggested proce- dure to other procedures for unstructured problems and in one case where the problem is assumed to have a well-defined geometric structure. Finally we discuss sensitivity of the procedure to the assumed value of p0 and suggest a generalization. 1. Introduction. Single sequence problems of change-point detection have a long history in industrial quality control, where an observed process is assumed initially to be in control and at a change-point becomes out of control. It is desired to detect the change-point with as little delay as possible, subject to the constraint that false detections occurring before the true change-point are very rare. Outstand- ing early contributions are due to Page [7, 8], Shiryaev [18] and Lorden [5]. We assume there are parallel streams of data subject to change-points. More pre- cisely suppose that for each n = 1, . . . ,N , we make observations yn,t , t = 1,2, . . . . The observations are mutually independent within and across data streams. At a certain time κ , there are changes in the distributions of observations made at a subset N {1, . . . ,N} of cardinality |N | ≤ N . Also denote by N c the set of un- affected data streams. The change-point κ , the subset N and its size, and the size of the changes are all unknown. As in the case of a single sequence, N = 1, the goal is to detect the change-point as soon as possible after it occurs, while keeping the frequency of false alarms as low as possible. In the change-point detection litera- ture, a surrogate for the frequency of false alarms is the average-run-length (ARL), Received June 2012; revised November 2012. 1Supported in part by NSF Grant 1043204 and Stanford General Yao-Wu Wang Graduate Fellowship. MSC2010 subject classifications. 62L10. Key words and phrases. Change-point detection, multi-sensor. 670 SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 671 defined to be the expected time before incorrectly announcing a change of distri- bution when none has occurred. It may be convenient to imagine that the data streams represent observations at a collection of N sensors and that the change-point is the onset of a localized sig- nal that can be detected by sensors in the neighborhood of the signal. In this paper we assume for the most part that the problem is unstructured in the sense that we do not assume a model that relates the changes seen at the different sensors. An example of an unstructured problem is the model for anomaly detection in com- puter networks developed in [4]. For other discussions of unstructured problems and applications, see [2, 6, 9, 19]. At the other extreme are structured problems where there exists a profile deter- mining the relative magnitudes of the changes observed by different sensors, say, according to their distance from the location of a signal (e.g., [12, 16]). A prob- lem that is potentially structured, may behave more like an unstructured problem if the number of sensors is small and/or they are irregularly placed, so their dis- tances from one another are large compared to the point spread function of the signals. Alternatively, local signals may be collected at the relatively widely dis- persed hubs of small, star-shaped subnetworks, then condensed and transmitted to a central processor, thus in effect removing the local structure. The detection problem of particular interest in this paper involves the case that N is large and |N | is relatively small. To achieve efficient detection, the detection procedure should use insofar as possible only information from affected sensors and suppress noise from the unaffected sensors. In analogy to the well-known CUSUM statistic (e.g., Page [7, 8], Lorden [5]), Mei [6] recently proposed a multi-sensor procedure based on sums of the CUSUM statistic from individual sensors. He then compares the sum with a suitable thresh- old to determine a stopping rule. While the distributions of the data, both before and after the change-point, are completely general, they are also assumed to be completely known. The method is shown to minimize asymptotically the expected detection delay (EDD) for a given false alarm rate, when the threshold value (and hence the constraint imposed by the ARL) becomes infinitely large. The procedure fails to be asymptotically optimal when the specified distributions are incorrect. Tartakovsky and Veeravalli proposed a different procedure [19] that sums the local likelihood ratio statistic before forming CUSUM statistics. They also assume the post-change distributions are completely prescribed. Moreover, both procedures assume the change-point is observed by all sensors. When only a subset of sen- sors observe the change-point, these procedures include noise from the unaffected sensors in the detection statistic, which may lead to long detection delays. In this paper, we develop a mixture procedure that achieves good detection per- formance in the case of an unknown subset of affected sensors and incompletely specified post-change distributions. The key feature of the proposed procedure is that it incorporates an assumption about the fraction of affected sensors when com- puting the detection statistic. We assume that the individual observations are inde- pendent and normally distributed with unit variance, and that the changes occur 672 Y. XIE AND D. SIEGMUND in their mean values. At the t th vector of observations (yn,t , n = 1, . . . ,N), the mixture procedure first computes a generalized likelihood ratio (GLR) statistic for each individual sensor under the assumption that a change-point has occurred at k ≤ t . The local GLR statistics are combined via a mixture model that has the ef- fect of soft thresholding the local statistics according to an hypothesized fraction of affected sensors, p0. The resulting local statistics are summed and compared with a detection threshold. To characterize the performance of our proposed procedure, we derive analytic approximations for its ARL and EDD, which are evaluated by comparing the approximations to simulations. Since simulation of the ARL is quite time consuming, the analytic approximation to the ARL proves very useful in de- termining a suitable detection threshold. The proposed procedure is then compared numerically to competing procedures and is shown to be very competitive. It is also shown to be reasonably robust to the choice of p0, and methods are suggested to increase the robustness to mis-specification of p0. Although we assume throughout that the observations are normally distributed, the model can be generalized to an exponential family of distributions satisfying some additional regularity conditions. The remainder of the paper is organized as follows. In Section 2 we establish our notation and formulate the problem more precisely. In Section 3 we review several detection procedures and introduce the proposed mixture procedure. In Section 4 we derive approximations to the ARL and EDD of the mixture procedure, and we demonstrate the accuracy of these approximations numerically. Section 5 demon- strates by numerical examples that the mixture procedure performs well compared to other procedures in the unstructured problem. In Section 6 we suggest a “par- allel” procedure to increase robustness regarding the hypothesized fraction of af- fected data streams, p0. In Section 7 we also compare the mixture procedure to that suggested in [16] for a structured problem, under the assumption that the assumed structure is correct. Finally Section 8 concludes the paper with some discussion. 2. Assumptions and formulation. Given N sensors, for each n = 1,2, . . . , N , the observations from the nth sensor are given by yn,t , t = 1,2, . . . . Assume that different observations are mutually independent and normally distributed with unit variances. Under the hypothesis of no change, they have zero means. Probabil- ity and expectation in this case are denoted by P∞ and E∞, respectively. Alterna- tively, there exists a change-point κ , 0 ≤ κ < ∞, and a subset N of {1,2, . . . ,N}, having cardinality |N |, of observations affected by the change-point. For each n ∈ N , the observations yn,t have means equal to μn > 0 for all t > κ , while observations from the unaffected sensors keep the same standard normal distri- bution. The probability and expectation in this case are denoted by Pκ and Eκ , respectively. In particular, κ = 0 denotes an immediate change. Note that proba- bilities and expectations depend on N and the values of μn, although this depen- dence is suppressed in the notation. The fraction of affected sensors is given by p = |N |/N . SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 673 Our goal is to define a stopping rule T such that for all sufficiently large pre- scribed constants c > 0, E∞{T } ≥ c, while asymptotically Eκ{T κ|T > κ} is a minimum. Ideally, the minimization would hold uniformly in the various unknown parameters: κ , N and the μn. Since this is clearly impossible, in Section 5 we will compare different procedures through numerical examples computed under vari- ous hypothetical conditions. 3. Detection procedures. Since the observations are independent, for an as- sumed value of the change-point κ = k and sensor n ∈ N , the log-likelihood of observations accumulated by time t > k is given by n(t, k,μn) = t∑ i=k+1 ( μnyn,i μ2n/2 ) .(3.1) We assume that each sensor is affected by the change with probability p0 (inde- pendently from one sensor to the next). The global log likelihood of all N sensors is N∑ n=1 log ( 1 p0 + p0 exp[n(t, k,μn)]).(3.2) Expression (3.2) suggests several change-point detection rules. One possibility is to set μn equal to a nominal change, say δ > 0, which would be important to detect, and define the stopping rule T1 = inf { t : max 0≤k≤t N∑ n=1 log ( 1 p0 + p0 exp[+n (t, k, δ)])≥ b } ,(3.3) where x+ denotes the positive part of x. Here thresholding by the positive part plays the role of dimension reduction by limiting the current considerations only to sequences that appear to be affected by the change-point. Another possibility is to replace μn by its maximum likelihood estimator, as follows. The maximum likelihood estimate of the post-change mean as a function of the current number of observations t and putative change-point location k is given by μ n = ( t∑ i=k+1 yn,i )+/ (t k).(3.4) Substitution into (3.1) gives the log generalized likelihood ratio (GLR) statistic. Putting Sn,t = t∑ i=1 yn,i, (3.5) Un,k,t = (t k) 1/2(Sn,t Sn,k), 674 Y. XIE AND D. SIEGMUND we can write the log GLR as n(t, k, μ n) = (U+n,k,t )2/2.(3.6) We define the stopping rule T2 = inf { t : max 0≤kN∑ n=1 log ( 1 p0 + p0 exp[(U+n,k,t )2/2])≥ b } .(3.7) REMARK. In what follows we use a window limited version of (3.7), where the maximum is restricted to m0 ≤ t k < m1 for suitable m0 < m1. The role of m1 is two-fold. On the one hand it reduces the memory requirements to implement the stopping rule, and on the other it effectively establishes a minimum level of change that we want to detect. For asymptotic theory given below, we assume that b → ∞, with m1/b also diverging. More specific guidelines in selecting m1 are discussed in [3]. In the numerical examples that follow, we take m0 = 1. In practice a slightly larger value can be used to provide protection against outliers in the data, although it may delay detection in cases involving very large changes. The detection rule (3.7) is motivated by the suggestion of [17] for a similar fixed sample change-point detection problem. For the special case p0 = 1, (3.7) becomes the (global) GLR procedure, which for N = 1 was studied by [14]. It is expected to be efficient if the change-point affects a large fraction of the sensors. At the other extreme, if only one or a very small number of sensors is affected by the change-point, a reasonable procedure would be Tmax = inf { t : max 0≤k( U+n,k,t )2 /2 ≥ b } .(3.8) The stopping rule Tmax can also be window limited. Still other possibilities are suggested by the observation that a function of y of the form log[1 p0 + p0 exp(y)] is large only if y is large, and then this function is approximately equal to [y + log(p0)]+. This suggests the stopping rules T3 = inf { t : max 0≤kN∑ n=1 [ n(t, k, δ) + log(p0)]+ ≥ b } (3.9) and T4 = inf { t : max 0≤kN∑ n=1 [( U+n,k,t )2 /2 + log(p0)]+ ≥ b } ,(3.10) or a suitably window limited version. Mei [6] suggests the stopping rule TMei = inf { t : N∑ n=1 max 0≤k} ,(3.11) SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 675 which simply adds the classical CUSUM statistics for the different sensors. Note that this procedure does not involve the assumption that all distributions affected by the change-point change simultaneously. As we shall see below, this has a negative impact on the efficiency of the procedure in our formulation, although it might prove beneficial in differently formulated problems. For example, there may be a time delay before the signal is perceived at different sensors, or there may be different signals occurring at different times in the proximity of different sensors. In these problems, Mei’s procedure, which allows changes to occur at different times, could be useful. The procedure suggested by Tartakovsky and Veeravalli [19] is defined by the stopping rule TTV inf { t : max 0≤kN∑ n=1 n(t, k, δ) ≥ b } .(3.12) This stopping rule resembles T3(p0) with p0 = 1, but with one important differ- ence. After a change-point the statistics of the unaffected sensors have negative drifts that tend to cancel the positive drifts from the affected sensors. This can lead to a large EDD. Use of the positive part, [n(t, k, δ)]+, in the definitions of our stopping rules is designed to avoid this problem. Different thresholds b are required for each of these detection procedures to meet the ARL requirement. 4. Properties of the detection procedures. In this section we develop theo- retical properties of the detection procedures T1 to T4, with emphasis on T2 and the closely related T4. We use two standard performance metrics: (i) the expected value of the stopping time when there is no change, the average run length or ARL; (ii) the expected detection delay (EDD), defined to be the expected stopping time in the extreme case where a change occurs immediately at κ = 0. The EDD pro- vides an upper bound on the expected delay after a change-point until detection occurs when the change occurs later in the sequence of observations. The approx- imation to the ARL will be shown below to be very accurate, which is fortunate since its simulation can be quite time consuming, especially for large N . Accu- racy of our approximation for the EDD is variable, but fortunately this parameter is usually easily simulated. 4.1. Average run length when there is no change. The ARL is the expected value of the stopping time T when there is no change-point. It will be convenient to use the following notation. Let g(x) denote a twice continuously differentiable increasing function that is bounded below at ∞ and grows sub-exponentially at +∞. In what follows we consider explicitly g(u). With an additional argument discussed below the results also apply to g(u+). Let ψ(θ) = logE{exp[θg(U)]},(4.1) 676 Y. XIE AND D. SIEGMUND where U has a standard normal distribution. Also let γ (θ
) = 12θ2E {[ g˙(U) ]2 exp [ θg(U) ψ(θ)]},(4.2) where the dot denotes differentiation. Let H(N, θ) = θ [2πψ¨(θ)] 1/2 γ (θ)N1/2 exp { N [ θψ˙(θ) ψ(θ)]}.(4.3) Denote the standard normal density function by φ(x) and its distribution function by (x). Also let ν(x) = 2x 2 exp[ 2∑∞1 n 1 ( |x|n1/2/2)]; cf. [13], page 82. For numerical purposes a simple, accurate approximation is given by (cf. [15]) ν(x) ≈ (2/x)[ (x/2) 0.5] (x/2) (x/2) + φ(x/2) . THEOREM 1. Assume that N → ∞ and b → ∞ with b/N fixed. Let θ be defined by ψ˙(θ) = b/N . For the window limited stopping rule T = inf { t : max 0≤kN∑ n=1 g(Un,k,t ) ≥ b } (4.4) with m1 = o(br) for some positive integer r , we have E ∞{T } ～ H(N, θ) /∫ [2Nγ (θ)/m0]1/2 [2Nγ (θ)/m1]1/2 yν2(y) dy.(4.5) REMARK. The integrand in the approximation is integrable at both 0 and ∞ by virtue of the relations ν(y) → 1 as y → 0, and ν(y) ～ 2/y2 as y → ∞. The following calculations illustrate the essential features of approxima- tion (4.5). For detailed proofs in similar problems, see [14] (where additional complications arise because the stopping rule there is not window limited) or [16]. From arguments similar to those used in [17], we can show that P ∞{T ≤ m} = P∞ { max t≤m,m0≤t k≤m1 N∑ n=1 g(Un,k,t ) ≥ b } (4.6) ～ N2e N[θψ˙(θ) ψ(θ)][2πNψ¨(θ)] 1/2|θ | 1γ 2(θ) × ∫ m1/m m0/m ν2 ([ 2Nγ (θ)/(mt) ]1/2) (1 t) dt/t2. SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 677 Here it is assumed that m is large, but small enough that the right-hand side of (4.6) converges to 0 when b → ∞. Changing variables in the integrand and using the notation (4.3), we can re-write this approximation as P ∞{T ≤ m} ～ m ∫ [2Nγ (θ)/m0]1/2 [2Nγ (θ)/m1]1/2 yν2(y) dy/H(N, θ).(4.7) From the arguments in [14] or [16] (see also [1]), we see that T is asymptotically exponentially distributed and is uniformly integrable. Hence if λ denotes the factor multiplying m on the right-hand side of (4.7), then for still larger m, in the range where mλ is bounded away from 0 and ∞, P∞{T ≤ m} [1 exp( λm)] → 0. Consequently E∞{T } ～ λ 1, which is equivalent to (4.5). REMARKS. (i) The result we have used from [17] was motivated by a prob- lem involving changes that could be positive, or negative, or both; and in that paper it was assumed that the function g(u) is twice continuously differentiable. The required smoothness is not satisfied by the composite functions of principal interest here, of the form g(u+). However, (i) the required smoothness is required only in the derivation of some of the constant factors, not the exponentially small factor, and (ii) the second derivative that appears in the derivation in [17] can be eliminated from the final approximation by an integration by parts. As a conse- quence, we can approximate the indicator of u > 0 by (ru) and use in place of u+ the smooth function ∫ u ∞ (rv)dv = u (ru) + r 1φ(ru), which converges uni- formly to u+ as r → ∞. Letting r → ∞ and interchanging limits produce (4.6). An alternative approach would be simply to define g(u) to be appropriate for a one-sided change while having the required smoothness in u. An example is g(u) = log[1 p0 +p0 exp(u2 (ru)/2)], which sidesteps the technical issue, but seems less easily motivated. (ii) The fact that all the stopping times studied in this paper are asymptotically exponentially distributed when there is no change can be very useful. (A simulation illustrating this property in the case of T2 is given in Section 4.3.) To simulate the ARL, it is not necessary to simulate the process until the stopping time T , which can be computationally time consuming, but only until a time m when we are able to estimate P∞{T ≤ m} with a reasonably small percentage error. For the numerical examples given later, we have occasionally used this shortcut with the value of m that makes this probability 0.1 or 0.05. (iii) Although the mathematical assumptions involve large values of N , some numerical experimentation for T2(p0) shows that (4.5) gives roughly the correct values even for N = 1 or 2. For p0 = 1 (4.5) provides numerical results similar to those given for the generalized likelihood ratio statistic in [14]. (iv) Theorem 1 allows us to approximate the ARL for T2 and T4. The stop- ping rule Tmax is straightforward to handle, since the minimum of N independent exponentially distributed random variables is itself exponentially distributed. The 678 Y. XIE AND D. SIEGMUND stopping rules T1 and T3, where g is composed with +t,k,δ , can be handled by a similar argument with one important difference. Now the cumulant generating function ψ(θ) depends on w = t k, so the equation defining θ must be solved for each value of w, and the resulting approximation summed over possible values of w. Fortunately only a few terms make a substantial contribution to the sum, except when δ is very small. For the results reported below, the additional amount of computation is negligible. 4.2. Expected detection delay. After a change-point occurs, we are interested in the expected number of additional observations required for detection. For the detection rules considered in this paper, the maximum expected detection delay over κ ≥ 0 is attained at κ = 0. Hence we consider this case. Here we are unable to consider stopping times defined by a general function g, so we consider the specific functions involved in the definitions of T2 and T4. Let g(u,p0) = log(1 p0 +p0 exp[(u+)2/2]) or [(u+)2/2+ log(p0)]+, and let U denote a standard normal random variable. Recall that N denotes the set of sensors at which there is a change, |N | is the cardinality of this set and p = |N |/N is the true fraction of sensors that are affected by the change. For each n ∈ N the mean value changes from 0 to μn > 0, and for n ∈ N c the distribution remains the same as before the change-point. Let
= (∑ n∈N μ2n )1/2 .(4.8) Note that the Kullback–Leibler divergence of a vector of observations after the change-point from a vector of observations before the change-point is
2/2, which determines the asymptotic rate of growth of the detection statistic after the change- point. Using Wald’s identity [13], we see to a first-order approximation that the expected detection delay is 2b/
2, provided that the maximum window size, m1, is large compared to this quantity. In the following derivation we assume m1 2b/
2. In addition, let S t t∑ i=1 zi(4.9) be a random walk where the increments zi are independent and identically dis- tributed with mean
2/2 and variance
2. Let τ = min{t : S t > 0}. Our approxi- mation to the expected detection delay given below depends on two related quan- tities. The first is ρ(
) = 12E { S 2τ } /E{S τ }(4.10) SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 679 for which exact computational expressions and useful approximations are available in [13]. In particular, ρ(
) = E{z21}/(2E{z1}) ∞∑ i=1 i 1E { S i }=
2/4 + 1 ∞∑ i=1 i 1E { S i } ,(4.11) where (x) = min{x,0}. The second quantity is E{mint≥0 S t }, which according to (Problem 8.14 in [13]) is given by E { min t≥0 S t } = ρ(
) 1
2/4.(4.12) The following approximation refines the first-order result for the expected de- tection delay. Recall that E0 denotes expectation when the change-point κ = 0. THEOREM 2. Suppose b → ∞, with other parameters held fixed. Then for T = T2 or T4, E 0{T } = 2
2 [ b + ρ(
) |N | logp0 |N |/2 + E { min t≥0 S t } (4.13) (N |N |)E{g(U,p0)}+ o(1)]. The following calculation provides the ingredients for a proof of (4.13). For de- tails in similar problems involving a single sequence, see [10] and [14]. For conve- nience we assume that T = T2, but there is almost no difference in the calculations when T = T4. Let k0 = b1/2. For k < T k0, we can write the detection statistic at the stopping time T as follows, up to a term that tends to zero exponentially fast in probability: Zk,T = N∑ n=1 g(Un,k,T ,p0) = ∑ n∈N g(Un,k,T ,p0) + ∑ n∈N c g(Un,k,T ,p0) = ∑ n∈N log ( p0 exp {( U+n,k,T )2 /2 }[ 1 + 1 p0 p0 exp { (U+n,k,T )2/2} ]) + ∑ n∈N c g(Un,k,T ,p0) = ∑ n∈N [ logp0 + (U+n,k,T )2/2]+ ∑ n∈N c g(Un,k,T ,p0)(4.14) + ∑ n∈N log ( 1 + 1 p0 p0 exp { (U+n,k,T )2/2} ) 680 Y. XIE AND D. SIEGMUND = |N | logp0 + ∑ n∈N ( U+n,k,T )2 /2 + ∑ n∈N c g(Un,k,T ,p0) + o(1) = |N | logp0 + ∑ n∈N [ (Sn,T Sn,k)+]2/2(T k) + ∑ n∈N c g(Un,k,T ,p0) + o(1). The residual term ∑ n∈N log(1 + (1 p0) exp{ (U+n,k,T )2/2}/p0) tends to zero exponentially fast when b → ∞ because when b → ∞, T → b/
, and n ∈ N , (U+n,k,T )2 grows on the order of μ2n(T k) > μ2nk0 = μ2n √ b. We then use the following simple identity to decompose the second term in (4.14) for the affected sensors into two parts: ( S+n,t )2 /2t = S2n,t /2t ( S n,t )2 /2t (4.15) = μn(Sn,t μnt/2) + (Sn,t μnt)2/2t (S n,t )2/2t. From the preceding discussion, we see that max0≤kwhile maxT k0≤kprobability the max over all k is attained for k < T k0, so from (4.15) and (4.14) we have max 0≤k= max 0≤kN∑ n=1 g(Un,k,T ,p0) + o(1) = |N | logp0 + max 0≤k[∑ n∈N μn [ (Sn,T Sn,k) (T k)μn/2] + ∑ n∈N [ (Sn,T Sn,k) (T k)μn]2/[2(T k)] [(Sn,T Sn,k) ]2/2(T k) (4.16) + ∑ n∈N c g(Un,k,T ,p0) ] + o(1) = |N | logp0 + ∑ n∈N μn(Sn,T T μn/2) + max 0≤k[ ∑ n∈N μn(Sn,k kμn/2) SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 681 + ∑ n∈N [ (Sn,T Sn,k) (T k)μn]2/[2(T k)] ∑ n∈N [ (Sn,T Sn,k) ]2/[2(T k)] + ∑ n∈N c g(Un,k,T ,p0) ] + o(1). The following lemma forms the basis for the rest of the derivation. The proof is omitted here; for details see [20] (or [14] for the special case N = 1). LEMMA 4.1. For k0 = b1/2, asymptotically as b → ∞ max 0≤k[ ∑ n∈N μn(Sn,k kμn/2) + ∑ n∈N [(Sn,T Sn,k) (T k)μn]2 2(T k) ∑ n∈N [(Sn,T Sn,k) ]2 2(T k) + ∑ n∈N c g(Un,k,T ,p0) ] = ∑ n∈N (Sn,T T μn)2/2T + ∑ n∈N c g(Un,0,T ,p0) + max 0≤k[ ∑ n∈N μn(Sn,k kμn/2) ] + op(1), where op(1) converges to 0 in probability. By taking expectations in (4.16), letting b → ∞ and using Lemma 4.1, we have E 0 { max 0≤kN∑ n=1 g(Un,k,T ,p0) } = E0 { |N | logp0 + ∑ n∈N μn(Sn,T T μn/2) + ∑ n∈N (Sn,T T μn)2 2T (4.17) + ∑ n∈N c g(Un,0,T ,p0) + max 0≤k[ ∑ n∈N μn(Sn,k kμn/2) ]} + o(1). We will compute each term on the right-hand side of (4.17) separately. We will need the lemma due to Anscombe and Doeblin (see Theorem 2.40 in [13]), which states that the standardized randomly stopped sum of random variables are asymp- totically normally distributed under quite general conditions. 682 Y. XIE AND D. SIEGMUND (i) By Wald’s identity [13], E 0 {∑ n∈N μn(Sn,T T μn/2) } = E0{T }
2/2.(4.18) (ii) By the Anscombe–Doeblin lemma, (Sn,T T μn)/T 1/2 is asymptotically normally distributed with zero mean and unit variance. Hence ∑ n∈N (Sn,T T μn) 2/T is asymptotically a sum of independent χ21 random variables, so E 0 {∑ n∈N (Sn,T T μn)2/2T } = |N |/2 + o(1).(4.19) (iii) Similarly, E 0 { ∑ n∈N c g(Un,0,T ,p0) } → (N |N |)E0{g(U,p0)}.(4.20) (iv) The term ∑n∈N μn(Sn,k μnk/2) (k ≥ 0) is a random walk with neg- ative drift
2/2 and variance
2. Hence E0{max0≤k∑ n∈N μn(Sn,k kμn/2)]} converges to the expected minimum of this random walk, which has the same distribution as mint≥0 S t defined above. Having evaluated the right-hand side of (4.17), we now consider the left-hand side, to which we will apply a nonlinear renewal theorem. This requires that we write the process of interest as a random walk and a relatively slowly varying remainder, and follows standard lines by using a Taylor series approximation to show that for large values of t and bounded values of k (cf. [10, 14], and the argument already given above) the asymptotic growth of ∑Nn=1 g(Un,k,t , p0) for t > κ is governed by the random walk ∑ n∈N μn(Sn,t tμn/2), which has mean value t
2/2 and variance t
2. By writing E 0 { max 0≤kN∑ n=1 g(Un,k,T ,p0) } = b + E0 { max 0≤kN∑ n=1 g(Un,k,T ,p0) b } ,(4.21) and using nonlinear renewal theory to evaluate the expected overshoot of the pro- cess of (4.9) over the boundary ([13], Chapter IX), we obtain E 0 { max 0≤kN∑ n=1 g(Un,k,T ,p0) b } → ρ(
).(4.22) REMARKS. (i) Although the proof of Theorem 2 follows the pattern of argu- ments given previously in the case N = 1, unlike that case where the asymptotic approximation is surprisingly accurate even when the EDD is relatively small, here the accuracy is quite variable. The key element in the derivation is the asymptotic linearization of g(U+n,k,t , p0) for each n ∈ N into a term involving a random walk SEQUENTIAL MULTI-SENSOR CHANGE-POINT DETECTION 683 and a remainder. A simple test for conditions when the approximation will be rea- sonably accurate is to compare the exact value of E0{Z0,t }, which is easily evalu- ated by numerical integration, to the expectation of the linearized approximation, then take t large enough to make these two expectations approximately equal. If such a value of t makes the expectations less than or equal to b, the approximation of the theorem will be reasonably accurate. Indeed the preceding argument is sim- ply an elaboration of these equalities at t = T combined with Wald’s identity to extract E0{T } from the random walk, and numerous technical steps to approximate the nonnegligible terms in the remainders. For a crude, but quite reliable approx- imation that has no mathematically precise justification that we can see, choose t to satisfy E0{Z0,t } = b. Fortunately the EDD is easily simulated when it is small, which is where problems with the analytic approximation arise. (ii) In principle the same method can be used to approximate the expected de- tection delay of T1 or T3. In some places the analysis is substantially simpler, but in one important respect it is more complicated. In the preceding argument, for n ∈ N c the term involving the expected value of g(Un,0,T ,p0) is very simple, since U2n,0,T has asymptotically a χ2 distribution. For the stopping rules T1 and T3, the term g(n,0,T ,p0) does not have a limiting distribution, and in fact for n ∈ N c it converges to 0 as b → ∞. However, there are typically a large number of these terms, and in many cases T is relatively small, nowhere near its asymp- totic limit. Hence it would be unwise simply to replace this expectation by 0. To a crude first-order approximation T ～ b/[δ0(∑n∈N μn δ0/2)] = t0, say. Although it is not correct mathematically speaking, an often reasonable approximation can be obtained by using the term (N |N |)E0{n,0,t0} to account for the statis- tics associated with sequences unaffected by the change-point. Some examples are included in the numerical examples in Table 5. 4.3. Accuracy of the approximations. We start with examining the accuracy of our approximations for the ARL and the EDD in (4.5) and (4.13). For a Monte Carlo experiment we use N = 100 sensors, m1 = 200 and μn = 1 for all affected data streams. The comparisons for different values of p0 between the theoretical and Monte Carlo ARLs obtained from 500 Monte Carlo trials are given in Tables 1 and 2, which show that the approximation in (4.5) is quite accurate. Figure 1 illustrates the fact that T2(0.1) is approximately exponentially dis- tributed. Results for the EDD obtained from 500 Monte Carlo trials are given in Table 3. Although the approximation seems reasonable, it does not appear to be as accurate as the approximation for the ARL. Since the EDD requires considerably less com- putational effort to simulate and needs to be known only roughly when we choose design parameters for a particular problem, there is less value to an accurate ana- lytic approximation. We have performed considerably more extensive simulations that yield results consistent with the small experiments reported in Tables 1, 2 and 3. Since the 684 Y. XIE AND D. SIEGMUND TABLE 1 ARL of T2(p0),m1 = 200 p0 b Theory Monte Carlo 0.3 31.2 5001 5504 0.3 32.3 10,002 10,221 0.1 19.5 5000 4968 0.1 20.4 10,001 10,093 0.03 12.7 5001 4830 0.03 13.5 10,001 9948 TABLE 2 ARL of T4(p0), m1 = 200 p0 b Theory Monte Carlo 0.3 24.0 5000 5514 0.1 15.1 5000 5062 0.03 10.8 5000 5600 parameter p0 defining T2 must be chosen subjectively, it is interesting to observe that Table 3 suggests these procedures are reasonab