Minimax-Optimal Algorithms for Detecting Changes in Statistically Periodic Random Processes

TL;DR

This paper introduces minimax-optimal algorithms for detecting changes in statistically periodic random processes, with applications to real-world datasets like social media activity during events.

Contribution

It proposes a new class of stochastic processes for modeling periodicity and develops algorithms that are asymptotically optimal for change detection in these processes.

Findings

Algorithms achieve minimax asymptotic optimality as false alarm rate approaches zero.

Developed methods handle unknown post-change distributions and multiple observation streams.

Applied algorithms successfully to detect a running event from Instagram data.

Abstract

Theory and algorithms are developed for detecting changes in the distribution of statistically periodic random processes. The statistical periodicity is modeled using independent and periodically identically distributed processes, a new class of stochastic processes proposed by us. An algorithm is developed that is minimax asymptotically optimal as the false alarm rate goes to zero. Algorithms are also developed for the cases when the post-change distribution is not known or when there are multiple streams of observations. The modeling is inspired by real datasets encountered in cyber-physical systems, biology, and medicine. The developed algorithms are applied to sequences of Instagram counts collected around a 5K run in New York City to detect the run.

Equations

$$f_{n+T}=f_n,\forall n\ge 1.$$

$$\begin{array}{rl}{X}_n& \stackrel{\text{ind}}{\sim }p(\cdot ;{\theta }_n),n\ge 1.\end{array}$$

$$\begin{array}{rl}{\theta }_n& ={\theta }_{n+T},n\ge 1.\end{array}$$

$$X_n\sim p\left(\cdot ,\theta \left(\frac{n\%T}{T}\right)\right)n\ge 1.$$

$$X_n=\theta \left(\frac{n\%T}{T}\right)+Z_{n}n\ge 1,$$

$$g_{n+T}=g_n,\forall n\ge 1.$$

$$g_i\ne f_i,\text{for some }i=1,2,\cdots ,T.$$

$$X_n\sim \{\begin{array}{ll}{f}_n,& \forall n<\nu ,\\ g_n& \forall n\ge \nu .\end{array}$$

$$\begin{array}{rl}\underset{\tau }{\min }& \underset{\nu \ge 1}{\sup }{\mathsf{E}}_{\nu }[\tau -\nu |\tau \ge \nu ]\\ \text{subj. to}& {\mathsf{E}}_{\infty }[\tau ]\ge \beta ,\end{array}$$

$$\begin{array}{rl}\underset{\tau }{\min }& \underset{\nu \ge 1}{\sup }\text{ess}\sup {\mathsf{E}}_{\nu }[\tau -\nu |X_1,\cdots ,X_{\nu -1}]\\ \text{subj. to}& {\mathsf{E}}_{\infty }[\tau ]\ge \beta ,\end{array}$$

$$W_{n+1}=\underset{1\le k\le n+1}{\max }\sum _{i=k}^{n+1}\log \frac{g_i(X_i)}{f_i(X_i)}$$

$${\tau }_c=\inf \{n\ge 1:W_n>A\}.$$

$$W_{n+1}=W_n^{+}+\log \frac{g_{n+1}(X_{n+1})}{f_{n+1}(X_{n+1})},$$

$$I=\frac{1}{T}\sum _{i=1}^{T}D(g_i\Vert f_i),$$

$$D(g_i\Vert f_i)<\infty ,\forall i=1,2,\cdots ,T,$$

$$0<D(g_i\Vert f_i),\text{ for some }i=1,2,\cdots ,T.$$

$$\begin{array}{rl}\underset{\nu \ge 1}{\sup }\text{ess}& \sup {\mathsf{E}}_{\nu }[\tau -\nu |X_1,\cdots ,X_{\nu -1}]\\ & \ge \underset{\nu \ge 1}{\sup }{\mathsf{E}}_{\nu }[\tau -\nu |\tau \ge \nu ]\ge \frac{\log \beta }{I}(1+o(1)),\end{array}$$

$${\mathsf{E}}_{\infty }[{\tau }_c]\ge \beta ,$$

$$\begin{array}{rl}\underset{\nu \ge 1}{\sup }& {\mathsf{E}}_{\nu }[{\tau }_c-\nu |{\tau }_c\ge \nu ]\\ & \le \underset{\nu \ge 1}{\sup }\text{ess}\sup {\mathsf{E}}_{\nu }[{\tau }_c-\nu |X_1,\cdots ,X_{\nu -1}]\\ & \le \frac{A}{I}(1+o(1))=\frac{\log \beta }{I}(1+o(1)).\end{array}$$

$$W_{n+1}^{(\ell )}={\left(W_n^{(\ell )}\right)}^{+}+\log \frac{g_{n+1}^{(\ell )}(X_{n+1})}{f_{n+1}(X_{n+1})},$$

$${\tau }_{c\ell }=\inf \left\{n\ge 1:W_n^{(\ell )}\ge \log (\beta M)\right\},$$

$${\tau }_{cm}=\inf \left\{n\ge 1:\underset{1\le \ell \le M}{\max }{W}_n^{(\ell )}\ge \log (\beta M)\right\}.$$

$${\tau }_{cm}=\underset{\ell =1,\cdots ,M}{\min }{\tau }_{c\ell }.$$

$$R_n=\sum _{\ell =1}^M\sum _{k=1}^n\prod _{i=k}^n\frac{g_i^{(\ell )}(X_i)}{f_i(X_i)}$$

$${\tau }_{sr}=\inf \left\{n\ge 1:R_n\ge \beta M\right\}.$$

$${\tau }_{sr}\le {\tau }_{cm}.$$

$${\mathsf{E}}_{\infty }[{\tau }_{cm}]\ge {\mathsf{E}}_{\infty }[{\tau }_{sr}]\ge \beta .$$

$$I_{\ell }=\frac{1}{T}\sum _{i=1}^{T}D(g_i^{(\ell )}\Vert f_i),$$

$$\begin{array}{rl}\underset{\nu \ge 1}{\sup }& {\mathsf{E}}_{\nu }[{\tau }_{cm}-\nu |{\tau }_{cm}\ge \nu ]\\ & \le \underset{\nu \ge 1}{\sup }\text{ess}\sup {\mathsf{E}}_{\nu }[{\tau }_{cm}-\nu |X_1,\cdots ,X_{\nu -1}]\\ & \le \frac{\log \beta }{I_{\ell }}(1+o(1)).\end{array}$$

$$D_{n+1}^{(\ell )}={\left(D_n^{(\ell )}\right)}^{+}+\log \frac{g_{n+1}^{(\ell )}(X_{n+1,\ell })}{f_{n+1}^{(\ell )}(X_{n+1,\ell })},$$

Full text

000Taposh Banerjee is with the Department of Electrical and Computer Engineering, University of Texas at San Antonio, San Antonio, TX. Prudhvi Gurram is with Army Research Laboratory, Adelphi, MD, and with Booz Allen Hamilton, McLean, VA. Gene Whipps is with the Army Research Laboratory, Adelphi, MD. (e-mail: taposh.banerjee@utsa.edu; gene.t.whipps.civ@mail.mil; gurram_prudhvi@bah.com) The work of Taposh Banerjee was partially supported by a grant from the U.S. Army Research Laboratory, W911NF1820295, and by the Cloud Technology Endowed Professorship IV. A preliminary version of this paper has been presented at IEEE ICASSP 2019.

I Introduction

The problem of detecting an abrupt change in the statistical properties of a measurement process has many applications in engineering and sciences including statistical process control [1, 2, 3, 4], sensor networks [5, 6, 7, 8, 9, 10], computer networks [11, 12], public health [13, 14, 15], and telecommunication engineering [16, 17]. In many of these applications, algorithms are needed to detect the changes in real time, that is the changes have to detected as soon as they occur. Mathematical problems for real-time detection are studied in an area of statistics called sequential analysis [18, 19, 20]. Specifically, in the problem of quickest change detection (QCD), theory and algorithms are developed to detect a change in the distribution of a sequence of random variables with the minimum possible delay subject to a contraint on the rate of false alarms [21, 22, 23, 24, 25, 26, 27, 28]. See [29, 30, 31] for surveys.

In this paper, we study the problem of change detection for applications where the observation process exhibits statistically periodic behavior. Below we give three such applications:

Cyber-physical systems: In Fig. 1, we have plotted the average number of Instagram messages posted near three CCTV cameras in New York City (NYC). These counts were extracted from multi-modal traffic data (CCTV images, Twitter and Instagram messages) that we collected in NYC [32], [33]. We collected data during a 5K run that occurred on September 24, 2017. We also collected data on two Sundays before and one Sunday after the event day. The first two plots in Fig. 1 are for two CCTV cameras that were not on the path of the 5K run while the third figure contains data corresponding to a camera on the path of the run. As can be observed from the figures, the average message counts show similar statistical behavior on the four Sundays for the off-path cameras. For the on-path camera, the average message count increases due to the event. This periodic pattern was also observed in other data that we collected in [32] and [33]. Thus, the problem of anomaly detection or change detection in traffic data can be posed as the problem of detecting deviations away from statistically periodic behavior. 2. 2.

Neuroscience: In Fig.2, we have plotted the neural spike data (binned) collected from a brain-computer interface (BCI) study on mice [34]. In this BCI experiment, it is of interest to investigate if an observer mouse learns to associate a cue to a shock given to a target mouse. In the first $15$ trials (below dashed horizontal line), no shocks are given. Starting trial number $16$, a cue is followed by a shock. The time of the cue is indicated by the dashed vertical line in the figure. The change point problem of interest here is to detect a change in the neural firing patterns caused by the behavioral learning starting trial number $16$. If a change is detected, then further shock trials can be avoided. The neural spike data shows statistical periodicity across trials as the observer mouse has to be part of the same experiment in every trial. For further discussion, we refer the readers to [35].

Medicine: Signals collected in electrocardiography (ECG) show regular patterns of P waves, QRS complexes, and ST segments; see Fig. 3. It is well-known that arrhythmia can change the shape of this pattern. Typically, a 12-lead ECG is used to record a patient’s ECG data and a human expert reads a chart to detect the anomaly. This process can be automated and can be done in real time using a change-point detection algorithm.

The classical QCD literature can be broadly classified into two categories: results for independent and identically distributed (i.i.d.) processes with algorithms that can be computed recursively and enjoy strongly optimality properties [24], and results for non-i.i.d. data with algorithms that are hard to compute but are asymptotically optimal [25, 26, 27, 28]. As can be seen in the figures (Fig. 1–Fig. 3), the data encountered in applications of interest cannot be modeled as i.i.d. processes. To model statistically periodic data, such as shown in these figures, we use a new class of stochastic processes called independent and periodically identically distributed processes (i.p.i.d.). This is a class of processes introduced by us in [37]. For this class of non-i.i.d. processes, we will show in this paper that the optimal algorithms can be computed recursively and using finite memory. This makes the developed algorithms amenable to implementation in various applications in cyber-physical systems, medicine, and biology.

In [37], we provided a Bayesian optimality theory of change detection in i.p.i.d. processes where it was assumed that a prior distribution of the change point variable is available to the decision maker. Since such an assumption is typically not satisfied in practice, in this paper we take a non-Bayesian approach. Specifically, we obtain algorithms that are asymptotically optimal for two minimax stochastic optimization formulations [22], [23]. We also develop algorithms when the post-change distribution is unknown or when there are multiple streams of observations. Finally, we apply the developed algorithms to the NYC traffic data to detect the 5K run.

II Model and Problem Formulation

In this section, we discuss the concept of an i.p.i.d. process which we will use to model statistically periodic data. We also discuss a change point model and two minimax problem formulations.

II-A Statistical Model

An i.i.d. process is a sequence of random variables that are independent and have the same distribution. To model the periodic statistical behavior of data, we use a new class of stochastic processes called i.p.i.d. processes which are defined as follows.

Definition 1**.**

Let $\{X_{n}\}$ be a sequence of random variables such that the variable $X_{n}$ has density $f_{n}$. The stochastic process $\{X_{n}\}$ is called independent and periodically identically distributed (i.p.i.d) if $X_{n}$ are independent and there is a positive integer $T$ such that the sequence of densities $\{f_{n}\}$ is periodic with period $T$:

[TABLE]

We say that the process is i.p.i.d. with the law $(f_{1},\cdots,f_{T})$.

If $T=1$, then the i.p.i.d. process is an i.i.d. process. The law of an i.p.i.d. process is completely characterized by the finite-dimensional product distribution involving $(f_{1},\cdots,f_{T})$. We assume that in a normal regime, the data can be modeled as an i.p.i.d. process. At some point in time, due to an anomaly, the distribution of the i.p.i.d. process deviates from $(f_{1},\cdots,f_{T})$. Our objective in this paper is to develop algorithms to process $\{X_{n}\}$ in real time and detect changes in the distribution as quickly as possible, subject to a constraint on the rate of false alarms.

In practice, the following parametric form of an i.p.i.d. process can be used. Let $\{X_{n}\}$ is an independent sequence of random variables with distribution in a parametric family with parameters $\{\theta_{n}\}$:

[TABLE]

The process $\{X_{n}\}$ is an i.p.i.d. process (a parametric i.p.i.d. process) if there is an integer $T>0$ such that the parameter sequence $\{\theta_{n}\}$ is periodic with period $T$:

[TABLE]

Note that the statistical model in (2) has only $T$ parameters $\theta_{1},\cdots,\theta_{T}$. The change detection problem in this case reduces to detecting changes in these $T$ parameters. Given the parameters $\theta_{1},\cdots,\theta_{T}$, we can use an algorithm to observe the process $\{X_{n}\}$ sequentially over time $n$ and detect any changes in the values of any of the parameters. The baseline parameters in the problem, the period $T$ and the parameters within a period $\theta_{1},\cdots,\theta_{T}$, can be learned from the training data. A special and important case of a parametric i.p.i.d. process is when we have a smooth function $\theta(t),t\in[0,1],$ and

[TABLE]

where $n\%T$ represents $n$ modulo $T$. An example is a regression set up:

[TABLE]

where $\{Z_{n}\}$ is a zero-mean i.i.d. sequence. The change detection problem for the regression setup is the problem of detecting changes in the regression function $\theta(t)$.

Note that the sequence model in (2) is different from the sequence model studied in [38] and [39]. In the model studied in [38] and [39], the random variables $\{X_{n}-\theta_{n}\}$ are modeled as Gaussian random variables and the parameters $\{\theta_{n}\}$ are not periodic. Furthermore, the problem there is of simultaneous estimation of all the different parameters $\{\theta_{n}\}$ given all the observations $\{X_{n}\}$. That is, the problem is not sequential in nature. It is also not a change point problem.

An i.p.i.d. process is a cyclostationary process [40]. Although a cyclostationary process can also be used to model statistically periodic behavior, the i.p.i.d. definition captures sample level behavior and the independence assumption across time allows for the development of a strong change point detection theory.

II-B Change Point Model

To define a change point model, consider another periodic sequence of densities $\{g_{n}\}$ such that

[TABLE]

Thus, we essentially have $T$ distinct set of densities $(g_{1},\cdots,g_{T})$. We assume that at some point in time $\nu$, called the change point in the following, the law of the i.p.i.d. process is governed not by the densities $(f_{1},\cdots,f_{T})$, but by the new set of densities $(g_{1},\cdots,g_{T})$. These densities need not be all different from the set of densities $(f_{1},\cdots,f_{T})$, but we assume that there exists at least an $i$ such that they are different:

[TABLE]

The change point model is as follows. At a time point $\nu$, the distribution of the random variable changes from $\{f_{n}\}$ to $\{g_{n}\}$:

[TABLE]

We emphasize that the densities $\{f_{n}\}$ and $\{g_{n}\}$ are periodic. This model is equivalent to saying that we have two i.p.i.d. processes, one governed by the densities $(f_{1},\cdots,f_{T})$ and another governed by the densities $(g_{1},\cdots,g_{T})$, and at the change point $\nu$, the observation process switches from one i.p.i.d. process to another.

II-C Problem Formulation

We want to detect the change described in (6) as quickly as possible, subject to a constraint on the rate of false alarms. We are looking for a stopping time $\tau$ for the process to minimize a metric on the delay $\tau-\nu$ and to avoid the event of false alarm $\{\tau<\nu\}$. Specifically, we are interested in the popular false alarm and delay metrics of Pollak [23] and Lorden [22]. Let $\mathsf{P}_{\nu}$ denote the probability law of the process $\{X_{n}\}$ when the change occurs at time $\nu$ and let $\mathsf{E}_{\nu}$ denote the corresponding expectation. When there is no change, we use the notation $\mathsf{E}_{\infty}$. The quickest change detection problem formulation of Pollak [23] is defined as

[TABLE]

where $\beta$ is a given constraint on the mean time to false alarm. Thus, the objective is to find a stopping time $\tau$ that minimizes the worst case conditional average detection delay subject to a constraint on the mean time to false alarm. A popular alternative is the worst-worst case delay metric of Lorden [22]:

[TABLE]

where $\text{ess}\sup$ is used to denote the supremum of the random variable $\mathsf{E}_{\nu}[\tau-\nu|X_{1},\cdots,X_{\nu-1}]$ outside a set of measure zero. Further motivation and comparison of these and other problem formulations for change point detection can be found in the literature [29], [31], [30], [25].

In Section III, we develop algorithms and optimality theory for detecting changes in i.p.i.d. processes. In Section IV, we extend the results to the case when the post-change i.p.i.d. law $(g_{1},\cdots,g_{T})$ is unknown. In Section V, we study the distributed case where there are multiple parallel streams of i.p.i.d. processes and the change can occur in any one of them.

III Change Detection in a Single Sequence with Known Post-Change Law

We now propose a CUSUM-type scheme to detect the above change. This algorithm belongs to the class of generalized CUSUM schemes discussed in the literature [25]. We compute the sequence of statistics

[TABLE]

and raise an alarm as soon as the statistic is above a threshold $A$:

[TABLE]

We show below that this scheme is asymptotically optimal for the minimax problem formulations in (7) and (8). We first show that the statistic $W_{n}$ can be computed recursively and using finite memory.

Lemma 1**.**

The statistic sequence $\{W_{n}\}$ can be recursively computed as

[TABLE]

where $(x)^{+}=\max\{x,0\}$. Further, since the set of pre- and post-change densities $(f_{1},\cdots,f_{T})$ and $(g_{1},\cdots,g_{T})$ are finite, the recursion (11) can be computed using finite memory needed to store these $2T$ densities, the current observation, and the past statistic.

Proof.

The proof is provided in the appendix. ∎

In the rest of the paper, we refer to (11) to as the Periodic-CUSUM algorithm. Towards proving the optimality of the Periodic-CUSUM scheme, we obtain a universal lower bound on the performance of any stopping time for detecting changes in i.p.i.d. processes. Define

[TABLE]

where $D(g_{i}\;\|\;f_{i})$ is the Kullback-Leibler divergence between the densities $g_{i}$ and $f_{i}$. We assume that

[TABLE]

and

[TABLE]

Theorem 1**.**

Let the information number $I$ as defined in (12) satisfy $0<I<\infty$. Then, for any stopping time $\tau$ satisfying the false alarm constraint $\mathsf{E}_{\infty}[\tau]\geq\beta$, we have as $\beta\to\infty$

[TABLE]

where an $o(1)$ term is one that goes to zero in the limit as $\beta\to\infty$.

Proof.

The proof is provided in the appendix. ∎

We now show that the Periodic-CUSUM scheme (9)–(11) is asymptotically optimal for both the formulations (7) and (8).

Theorem 2**.**

Let the information number $I$ as defined in (12) satisfy $0<I<\infty$. Then, the Periodic-CUSUM stopping time $\tau_{c}$ (9)–(11) with $A=\log\beta$ satisfies the false alarm constraint, i.e.,

[TABLE]

and as $\beta\to\infty$,

[TABLE]

Proof.

The proof is provided in the appendix. ∎

IV Change Detection With Unknown Post-Change I.P.I.D. Law

In the previous section, we assumed that the post-change law $(g_{1},\cdots,g_{T})$ is known to the decision maker. This information was used to design the Periodic-CUSUM algorithm (11). In practice, this information may not be available. We now show that if the post-change law belongs to a finite set of $M$ possible distributions, $(g_{1}^{(1)},\cdots,g_{T}^{(1)}),\cdots,(g_{1}^{(M)},\cdots,g_{T}^{(M)})$, then an asymptotically optimal test can be designed.

For $\ell\in\{1,\cdots,M\}$, define the statistic

[TABLE]

and the stopping rule

[TABLE]

which is the Periodic-CUSUM stopping rule for the $\ell$th post-change law $(g_{1}^{(\ell)},\cdots,g_{T}^{(\ell)})$. Now, define

[TABLE]

Then, note that

[TABLE]

The stopping rule $\tau_{cm}$ is the stopping rule under which we stop the first time any of the $\ell$ Periodic-CUSUMs is above the threshold $\log(\beta M)$.

We now show that this stopping rule is optimal for both Lorden’s and Pollak’s criteria. Towards this end, we define a Shiryaev-Roberts-type statistic

[TABLE]

and a Shiryaev-Roberts-type stopping rule

[TABLE]

Note that

[TABLE]

We have the following theorem.

Theorem 3**.**

The process $\{R_{n}-nM\}$ is a $\mathsf{P}_{\infty}$ martingale. If $\mathsf{E}_{\infty}[R_{\tau_{sr}}]<\infty,$ then

[TABLE]

Further, if $(g_{1}^{(\ell)},\cdots,g_{T}^{(\ell)})$ is the true post-change i.p.i.d. law and

[TABLE]

then as $\beta\to\infty$,

[TABLE]

Proof.

The proof is provided in the appendix. ∎

Given the lower bound in Theorem 1, the stopping rule $\tau_{cm}$ is thus asymptotically optimal with respect to the criteria of Lorden and Pollak, uniformly over each possible post-change hypothesis $(g_{1}^{(\ell)},\cdots,g_{T}^{(\ell)})$, $\ell=1,\cdots,M$.

V Change Detection in a Distributed I.P.I.D. Setting

In the previous sections, we assumed that there is a single sequence of random variable $\{X_{n}\}$. In many applications, as in the problem of event detection in NYC data, the sensors are distributed. An event can occur near any one of the sensors. Thus, it is of interest to develop algorithms to detect event using multiple streams of data. In this section, we obtain optimal algorithms for detecting changes when the observation process in each stream is an i.p.i.d. process.

Let there be $M$ independent streams of data and let $\{X_{n,\ell}\}_{n=1}^{\infty}$ be the i.p.i.d. observation process of the $\ell$th stream with law $(f_{1}^{(\ell)},\cdots,f_{T}^{(\ell)})$. At the change point $\nu$, the law of the i.p.i.d. process in one of the streams changes from $(f_{1}^{(\ell)},\cdots,f_{T}^{(\ell)})$ to $(g_{1}^{(\ell)},\cdots,g_{T}^{(\ell)})$. The objective is to detect this change in distribution.

For $\ell\in\{1,\cdots,M\}$, define the statistic

[TABLE]

and

[TABLE]

Note that this rule is different from that discussed in (15) because here the statistic $D_{n+1}^{(\ell)}$ utilizes a different stream of observations and a different pre-change distribution for each stream index $\ell$.

Now, define the Shiryaev-Roberts statistic

[TABLE]

and the Shiryaev-Roberts stopping rule

[TABLE]

We have the following theorem.

Theorem 4**.**

If $\mathsf{E}_{\infty}[S_{\tau_{srd}}]<\infty,$ then

[TABLE]

Further, if the change occurs in the $\ell$th stream and

[TABLE]

then as $\beta\to\infty$,

[TABLE]

Proof.

The proof is similar to the proof of Theorem 3 and is skipped. ∎

In view of the lower bound obtained in Theorem 1, the above theorem shows that the stopping rule $\tau_{dm}$ is asymptotically optimal for each post-change scenario.

VI Learning a Parametric I.P.I.D. Model

In practice, learning pre- and post-change laws $(f_{1},\cdots,f_{T})$ and $(g_{1},\cdots,g_{T})$ can be hard. Thus, data will be typically modeled using a parametric i.p.i.d. model. We recall the definition of a parametric i.p.i.d. process from (1):

Definition 2**.**

Let $\{X_{n}\}$ be an independent sequence of random variables with distribution in a parametric family with parameters $\{\theta_{n}\}$:

[TABLE]

The process $\{X_{n}\}$ is called a parametric i.p.i.d. process if there is an integer $T>0$ such that the parameter sequence $\{\theta_{n}\}$ is periodic with period $T$:

[TABLE]

If the statistical properties of the observation process repeats every day, and if the data is collected once per hour, then in the above model, the period $T$ would correspond to $T=24$ hours in a day, and the variables $X_{1},\cdots,X_{T}$ would correspond to the data collected in the first day. In many applications, the data is often collected more frequently, at the rate of many samples per second. In such applications, $T$ could be, for example, equal to $24\times 60\times 60\times m$, where $m$ is the number of samples collected per second. In practice, it may be hard to learn a large number of parameters, and detect changes in them, especially if the post-change parameters are not known. The learning process can be made simpler by making additional assumptions on the model. We discuss one such simplification now.

In order to control the complexity of the problem, we may assume that the parameters are divided into batches and parameters in each batch are approximately constant. For example, a batch may correspond to data collected in an hour and the average count of objects may not change in an hour. Mathematically, we assume that in each cycle or period of length $T$, the vector of parameters $\{\theta_{k}\}_{k=1}^{T}$ is partitioned into $E$ batches or episodes. Specifically, for $N_{0}=0$ and positive integers $\{N_{e}\}_{e=1}^{E}$ we define $B_{e}=\{N_{e-1}+1,\cdots,N_{e}\}$ such that $\{1,\cdots,T\}=\cup_{e=1}^{E}B_{e},\quad B_{e}\cap B_{f}=\emptyset,\text{ for }e\neq f.$ For $e\in\{1,\cdots,E\}$, we define $\theta_{B_{e}}=(\theta_{N_{e-1}+1},\cdots,\theta_{N_{e}}).$ Thus, $\{\theta_{k}\}_{k=1}^{T}$ is partitioned as

[TABLE]

Note that we have $T=\sum_{e=1}^{E}|B_{e}|$.

We further assume a step model for parameters. Under this assumption, the parameters remain constant within a batch resulting in the step-wise constant sequence model

[TABLE]

That is $\theta^{(1)}=\theta_{1}=\cdots=\theta_{N_{1}}$, $\theta^{(2)}=\theta_{N_{1}+1}=\cdots=\theta_{N_{2}}$, and so on. Thus, if the batch sizes are large, there are only $E\ll T$ parameters to learn from the data. Also, we have $|B_{e}|$ samples for batch $e$. The objective is then to observe the process $\{X_{n}\}$ over time and detect any changes in the parameters $\theta^{(1)},\cdots,\theta^{(E)}$.

Note that for the regression model, the batch assumption above is equivalent to approximating the smooth function $\theta(t)$ by a step function.

VII Numerical Results

We now apply the stopping rule $\tau_{c}$ in (10) to the NYC Instagram count data. In Fig. 4, we have plotted the evolution of the test statistic $W_{n}$ applied to the Instagram counts collected near nine CCTV cameras. Two of these CCTV cameras are near or on the path of the run (called on-path cameras) and the rest are outside the race path (called off-path cameras). Out of the two on-path cameras, one is near the north end of the run and another is at the south end of the run. To obtain Fig. fig:AllTestConcatenatedCount, we arranged the Instagram counts near each CCTV camera in a concatenated fashion, with labeled segments separated via red vertical lines. Each day has $6598$ samples. To compute the statistic, we divided the data for each day into four batches, with the first three batches being of length $1500$. We modeled the data as a sequence of Poisson random variables. We used the count data from Sept. 10 (one of the non-event days) to learn the averages of these Poisson random variables for each of the four batches. We assumed that there is only one post-change parameter per batch with value equal to three times the value of the normal parameter for that batch. We then applied the test $\tau_{c}$ to all the four days of data. As observed in the top figure in Fig. 4, the test statistic stays close to zero for all the cameras on the non-event days. On the event day (Sept. 24th), the test statistic grows for the on-path cameras indicating that a deviation from the baseline was detected. The most significant change was observed near the CCTV camera near the north end of the run. In the bottom figure of Fig. 4, we have replotted the test statistic applied to the on-path camera near the south end of the run.

In Fig. 5, we have plotted the delay $\mathsf{E}_{1}[\tau_{c}]$ versus log of mean time to false alarm $\log\mathsf{E}_{\infty}[\tau_{c}]$ for the periodic-CUSUM algorithm. The simulation plot was obtained for the following set of pre- and post-change parameters:

[TABLE]

To obtain each of the five points in the figure for simulations, the value of the threshold $A$ in (10) was set to values $3,4,5,5.5$ and $6$ and both delay and false alarm estimates were obtained using $5000$ sample paths. The analysis plot was obtained by dividing the threshold by the average KL-divergence between the densities. In Fig. 6, we have plotted a typical evolution of the algorithm applied to simulated data. This plot was obtained for the same set of parameters specified in (34).

VIII Conclusions

We developed a minimax asymptotic theory for quickest detection of changes in i.p.i.d. models. We also studied the cases where the post-change i.p.i.d. law is unknown or where there are multiple streams of data. The algorithms developed were applied to count data extracted from multimodal data collected around a 5K run from NYC to detect the 5K run. The theoretical results show that many of the results valid for i.i.d. data models are also true for i.p.i.d. models. An important question one can ask motivated from [24] is whether the periodic-CUSUM algorithm is exactly optimal for Lorden’s formulation [22]. The answer to this question, we believe, is negative. Our belief is based on the Bayesian analysis done in [37] where we observed that single-threshold policies are not strictly optimal. In fact, in the Bayesian setting, the optimal algorithm has periodic thresholds. We conjecture that even for the minimax settings an algorithm with periodic thresholds will be strictly optimal.

Proof of Lemma 1.

For any sequence $\{Z_{i}\}$ of random variables, we can write

[TABLE]

Substituting $Z_{i}=\log\frac{g_{i}(X_{i})}{f_{i}(X_{i})}$ into the above equation we get the desired recursion for $W_{n}$ in (9):

[TABLE]

Note that the increment term $\log\frac{g_{n+1}(X_{n+1})}{f_{n+1}(X_{n+1})}$ is only a function of the current observation $X_{n+1}$. Also, since the processes are i.p.i.d. with laws $(f_{1},\cdots,f_{T})$ and $(g_{1},\cdots,g_{T})$, the likelihood ratio functions $\frac{g_{n}(\cdot)}{f_{n}(\cdot)}$ are not all distinct, and there are only $T$ such functions $\frac{g_{1}(\cdot)}{f_{1}(\cdot)}$ to $\frac{g_{T}(\cdot)}{f_{T}(\cdot)}$. Thus, we need only a finite amount of memory to store the past statistic, current observation, and $2T$ densities to compute this statistic recursively. ∎

Proof of Theorem 1.

Let $Z_{i}=\log\frac{g_{i}(X_{i})}{f_{i}(X_{i})}$ be the log likelihood ratio at time $i$. We show that the sequence $\{Z_{i}\}$ satisfies the following statement: as $n\to\infty$,

[TABLE]

where $I$ is as defined in (12). The lower bound then follows from Theorem 1 in [25]. Towards proving (36), note that as $n\to\infty$

[TABLE]

The above display is true because of the i.p.i.d. nature of the observation processes. This implies that as $n\to\infty$

[TABLE]

To show this, note that

[TABLE]

For a fixed $\epsilon>0$, because of (37), the LHS in (38) is greater than $I(1-\epsilon)$ for $n$ large enough. Also, let the maximum on the LHS be achieved at a point $k_{n}$, then

[TABLE]

Now $k_{n}$ cannot be bounded because of the presence of $n$ in the denominator. This implies $k_{n}>i$, for any fixed $i$, and $k_{n}\to\infty$. Thus, $\frac{1}{k_{n}}\sum_{i=\nu}^{\nu+k_{n}}Z_{i}\to I$. Since $k_{n}/n\leq 1$, we have that the LHS in (38) is less than $I(1+\epsilon)$, for $n$ large enough. This proves (38). To prove (36), note that due to the i.p.i.d. nature of the process

[TABLE]

The right hand side goes to zero because of (38) and because the maximum on the right hand side in (40) is over only finitely many terms. ∎

Proof of Theorem 2.

Again with $Z_{i}=\log\frac{g_{i}(X_{i})}{f_{i}(X_{i})}$, we show that the sequence $\{Z_{i}\}$ satisfies the following statement:

[TABLE]

The upper bound then follows from Theorem 4 in [25]. To prove (41), note that due to the i.p.i.d nature of the process we have

[TABLE]

The right hand side of the above equation goes to zero for any $\delta$ because of (37) and also because of the finite number of maximizations. The false alarm result follows directly from [25] with $A=\log\beta$ because the likelihood ratios here also form a $\mathsf{P}_{\infty}$ martingale. ∎

Proof of Theorem 3.

That $\{R_{n}-nM\}$ is a $\mathsf{P}_{\infty}$ martingale can be proved by direct verification. For the false alarm proof, we assume that $\mathsf{E}_{\infty}[\tau_{sr}]<\infty$, otherwise the proof is trivial. Since $\mathsf{E}_{\infty}[R_{\tau_{sr}}]<\infty$, we have that $R_{\tau_{sr}}-\tau_{sr}M$ is integrable. Further, as $n\to\infty$,

[TABLE]

Thus, by the optional sampling theorem [20] and (21) we have

[TABLE]

The delay result is true because of (18). ∎