Detecting Metachanges in Data Streams from the Viewpoint of the MDL Principle

Abstract

This paper addresses the issue of how we can detect changes of changes, which we call metachanges, in data streams. A metachange refers to a change in patterns of when and how changes occur, referred to as “metachanges along time” and “metachanges along state”, respectively. Metachanges along time mean that the intervals between change points significantly vary, whereas metachanges along state mean that the magnitude of changes varies. It is practically important to detect metachanges because they may be early warning signals of important events. This paper introduces a novel notion of metachange statistics as a measure of the degree of a metachange. The key idea is to integrate metachanges along both time and state in terms of “code length” according to the minimum description length (MDL) principle. We develop an online metachange detection algorithm (MCD) based on the statistics to apply it to a data stream. With synthetic datasets, we demonstrated that MCD detects metachanges earlier and more accurately than existing methods. With real datasets, we demonstrated that MCD can lead to the discovery of important events that might be overlooked by conventional change detection methods.

1. Introduction

1.1. Purpose of This Paper

In this study, we are concerned with detecting changes in data streams. The goal of change detection is to detect the time points at which the nature of the data-generating mechanism significantly changes.

Thus far, many algorithms have been proposed to detect change points in data streams (e.g., [1,2,3,4,5,6,7,8,9,10,11]), and several studies addressed or have been related to the issue of changes of changes [12,13,14,15,16,17,18]. In this paper, we refer to the changes of changes as metachanges. A metachange refers to a change in the pattern of when or how changes occur. It is practically important to detect metachanges because they may be early warning signals of important events [12,13]. Metachanges have been treated from a viewpoint of metachanges along time. Metachanges along time indicate that the interval significantly varies between the change points. Such metachanges were called burstiness [12] and volatility [13] in previous studies. The detection of metachanges along time provides users with useful information from data streams. For example, in a machine in a manufacturing factory, a decrease in the interval between change points might be a sign of a serious failure.

There is also another type of metachange: metachanges along state. Here, “state” refers to the parameter value of the probability density function of a distribution. We consider a situation where change points $t_1,\dots$ are detected for a data stream $y_1,y_2,\dots$, and $y_t$ is drawn from $p_y(y_t;\eta )$. Here, $p_y$ is a probability density function of distributions, and $\eta$ is the associated parameter. Note that $\eta$ is called state in this paper, and it varies before and after a change point. A metachange along state means a change of how significantly $\eta$ varies before and after a change point. Metachanges along state might provide information such as changes of magnitude and velocity, which indicate an important change in the underlying data-generating mechanism. For example, in a machine in a manufacturing factory, a shift to an abrupt (sudden) change from a gradual (incremental) change [19], or its inverse shift, might be a sign of serious events.

A conceptual illustration of metachanges is shown in Figure 1, where the upper graph shows a data stream $y_1,\dots$ and change points ${\left\{t_i\right\}}_{i=1}^8$ on the horizontal axis. The lower left graph shows intervals between change points $\Delta t=t_i-t_{i-1}$ on the vertical axis. Metachanges along time occur at $t_4,t_5,t_6,t_7$: for example, $t_4-t_3$ is different from $t_3-t_2$ and $t_2-t_1$. The lower right graph shows the states estimated piecewisely between the change points. Here, we assume $y_t$ is drawn from the univariate normal distribution $p_y(y_t;\mu ,\sigma )$, where $\mu$ is the mean and $\sigma$ is the standard deviation. In this case, $(\mu ,\sigma )$ is a state. In Figure 1, because there is no significant change in the magnitude of state change between $t_1$ and $t_2$, a metachange along state does not occur at $t_2$. However, there is a significant change in the magnitude of change of $\mu$ between $t_2$ and $t_3$: thus, a metachange along state occurs at $t_3$. Because the magnitudes of the changes of $\mu$ and $\sigma$ are almost the same between $t_3$ and $t_4$, a metachange along state does not occur at $t_4$. Using the same procedure, we conclude that metachanges along state occur at $t_3$ and $t_7$ with respect to $\mu$. Moreover, metachanges along state occur at $t_6$ and $t_8$ with respect to $\sigma$: the magnitude of the change of standard deviations around $t_6$ ($t_8$) is greater than those around $t_5$ ($t_7$). As a result, metachanges along state occur at $t_3,t_6,t_7,$ and $t_8$. We can infer that metachanges along both time and state occur at $t_6$ and $t_8$, by combining the metachanges along time and state.

Metachanges along time have been investigated in previous studies [12,13], and, although there have been several studies related to metachanges along state [14,15,16,17,18], the focus of these studies was not on metachanges along state in particular. The purpose of this paper is to propose a framework and an approach to detect metachanges along time and state from a unified view with the minimum description length (MDL) [20]. Therefore, our framework and approach not only include previous notions such as burstiness [12] and volatility [13] but also extend these notions to metachanges along state. MDL asserts that the best statistical decision strategy is the one that compresses the data best. Description and coding with MDL are suitable for quantifying changes, and they enable us to easily integrate the code lengths of time and state.

1.2. Related Work

Change detection has been extensively explored in the area of data mining. Thus far, several methods have been proposed to detect metachanges in data streams [12,13], and there have been several studies related to metachanges along state [14,15,16,17,18].

Kleinberg [12] and Huang et al. [13] proposed algorithms for detecting metachanges along time. Kleinberg [12] proposed an algorithm to detect bursts in a time series. This algorithm assumes that intervals between successive events are drawn from an exponential distribution. The discretized values of the parameters of the exponential distribution are regarded as states. For intervals between successive events, states are estimated with dynamic programming. Changes of state indicate changes of intervals between the successive events. Huang et al. [13] proposed an algorithm, called the volatility detector, which detects changes of rates of change. The volatility detector prepares two buckets, called the buffer and the reservoir, to store intervals between change points. The intervals are put into the buffer sequentially. When the buffer is full, an interval is dropped from the buffer and moved to the reservoir in a first-in-first-out fashion. The reservoir stores the dropped interval by randomly replacing one of its stored intervals. If the ratio of variances of the buffer and the reservoir is over or under the specified threshold, the algorithm judges that the intervals change between change points. The authors called this event volatility shift. Both the burst detector and volatility detector are assumed to be used in two steps. That is, change points are detected with other change detection algorithms, and then changes of intervals between the change points are detected. While the burst detector works in an offline fashion, the volatility detector works in an online fashion.

Moreover, there have been several studies related to metachanges along state [14,15,16,17,18]. Aggarwal [15] introduced velocity density estimation to understand, visualize, and determine trends in the evolution of fast data streams. Spiliopoulou et al. [16,17] proposed an algorithm, called MONIC, to model and track cluster transitions. Ntoutsi et al. [18] proposed an algorithm, called FINGERPRINT, to summarize cluster evolution. Huang et al. [14] proposed a change type detector, intended to categorize change types into three relative types, some of which correspond to concept drifts proposed in [19]. Although their algorithms [14,15,16,17,18] are related to metachanges along state, they are not intended to characterize and detect metachanges directly. In addition, many change detection algorithms have been proposed based on detecting changes of state (e.g., [6,7,8,9,21,22]). The dynamic model selection [6,7] is the seminal work to apply MDL to the task of dynamic model selection and change detection. The MDL change statistics [8], SCAW [9], and STREAMKRIMP [22] are change detection algorithms with MDL. However, these algorithms are not intended to characterize and detect metachanges along state directly.

1.3. Significance of This Paper

In the context of Section 1.1 and Section 1.2, the contributions of this paper are summarized in the following subsections.

1.3.1. Proposal of Concept of Metachange

To detect changes of changes in data streams, we define a concept of metachanges along both time and state. Previous studies [12,13] considered metachanges along time only. In this paper, we deal with metachanges along both time and state. Metachanges along time include the notions proposed in previous studies such as burstiness [12] and volatility [13]. Metachange along state could capture changes of changes of the parameters of distribution between change points.

Our concept of metachange can detect the potential change of changes in data streams, which was overlooked by previous studies.

1.3.2. Novel Algorithm for Detection of Metachanges

We define metachange statistics along both time and state. There is a challenge to combining the metachange statistics along time and those along state. In this paper, these statistics are defined based on the MDL principle. Metachange statistics along time (MCAT) is defined as the code length of an interval between the change points, whereas metachange statistics along state (MCAS) is defined as the difference between the predictive code length and the normalized maximum likelihood (NML) code length [23] after a change. It is possible to simply add these statistics because they are defined as code lengths, which enables us to detect metachanges along both time and state in a unified manner.

2. Theoretical Background of Metachange Statistics

In this section, we consider how to encode both intervals between change points and states around the change points. We assume that for a data stream $y_1,y_2,\dots$ change points $t_1,\dots$ are detected and that the intervals between change points $x_i=t_i-t_{i-1}$ and $y_t$ are drawn, respectively, from

$$x_i\sim p_x(x_i;\xi ),y_t\sim p_y(y_t;\eta ),$$

where $p_x$ and $p_y$ are probability density functions of distributions and $\xi$ and $\eta$ are the associated parameters. Finally, $\eta$ is the state whose metachanges are addressed in this paper.

2.1. Definitions of Metachanges

In this subsection, we give definitions of metachanges.

Definition 1.

(Metachange along time) For intervals between change points $x_1,x_2,\dots$, we say that a metachange along time occurs at a change point $t_i$ for a threshold parameter ${\delta }_{\mathrm{t}}>0$ if and only if

$$\begin{array}{c}{q}_1\to q_2\mathrm{at}t=t_i,\hfill \\ d(q_1,q_2)>{\delta }_{\mathrm{t}},q_1,q_2\in {\mathcal{F}}_{\mathrm{t}},{\mathcal{F}}_{\mathrm{t}}=\left\{p_x(x;\xi )\right\},\hfill \end{array}$$

(1)

where $q_1$ and $q_2$ are distributions of intervals. $q_1\to q_2$ means that $x_t\sim q_1$ at $t=t_{i-1}$ and $x_t\sim q_2$ at $t=t_i$. d is a distance function between the probability density functions.

Definition 2.

(Metachange along state) For a data stream $y_1,y_2,\dots$, we say that a metachange along state occurs at a change point $t_i$ for a threshold parameter ${\delta }_{\mathrm{s}}>0$ if and only if

$$\begin{array}{c}{q}_1\to q_2\mathrm{at}t=t_{i-1},q_2\to q_3\mathrm{at}t=t_i,\hfill \\ |d(q_2,q_3)-d(q_1,q_2)|>{\delta }_{\mathrm{s}},q_1,q_2,q_3\in {\mathcal{F}}_{\mathrm{s}},{\mathcal{F}}_{\mathrm{s}}=\left\{p_y(y;\eta )\right\},\hfill \end{array}$$

(2)

where $q_1$, $q_2$, and $q_3$ are distributions of values of the data stream. Equation (2) means that $y_t\sim q_1$ at $t=t_{i-2},\dots ,t_{i-1}-1$, $y_t\sim q_2$ at $t=t_{i-1},\dots ,t_i-1$, and $y_t\sim q_3$ at $t=t_i\dots ,t_{i+1}-1$. Here, d is the same as that in Definition 1.

Definition 3.

(Integrated metachange) For a change point $t_i$, we say that an integrated metachange occurs at $t_i$ if and only if Equation (1) or Equation (2) holds.

2.2. Problem Setting

In this subsection, we consider a situation where $(m+1)$ change points $t_1,\dots ,t_{m+1}$ are given. We consider how to encode $x_i$ and $y_t$ as shortly as possible. The ideal code length required for encoding $x_i$ is given by what we call the predictive code length, which is the sum of the negative logarithm of its predictive density $p_x$ at each time point, defined as follows:

$$\underset{{\left\{{\widehat{\xi }}_{x^{i-1}}\right\}}_{i=1}^m}{min}\sum _{i=1}^m-log{p}_x(x_i;{\widehat{\xi }}_{x^{i-1}}),$$

(3)

where ${\widehat{\xi }}_{x^{i-1}}$ are estimated at each change point. Similarly, the ideal code length required for encoding $y_t$ around change points is given by the predictive code length as follows:

$$\underset{{\left\{{\widehat{\eta }}_{y^{t-1}}\right|t\in \mathrm{neighbor}\left(t_i\right)\}}_{i=1}^m}{min}\sum _{i=1}^m\sum _{t\in \mathrm{Neighbor}\left(t_i\right)}-log{p}_y(y_t;{\widehat{\eta }}_{y^{t-1}}),$$

(4)

where $\mathrm{Neighbor}\left(t_i\right)$ indicates the neighborhood of a change point $t_i$. In practice, as explained in Section 3.3, $\mathrm{Neighbor}\left(t_i\right)=[t_i-h,t_i+h],h\in \mathbb{N}$. ${\widehat{\xi }}_{x^{i-1}}$ and ${\widehat{\eta }}_{y^{t-1}}$ are estimated using $x^{i-1}=x_1\dots x_{i-1}$ and $y^{t-1}=y_1\dots y_{t-1}$, respectively. A change of ${\widehat{\eta }}_{y^{t-1}}$ indicates a change of state. Detection of a metachange along time is asserted as a problem of detection of a change of ${\widehat{\xi }}_{x^{i-1}}$ in Equation (3). On the other hand, detection of a metachange along state is asserted as a problem of detection of a change of how ${\widehat{\eta }}_{y^{t-1}}$ in Equation (4) changes around a change point between change points.

3. Metachange Detection Algorithm

In this section, we present our online algorithm called metachange detection algorithm (MCD) for detecting metachanges along both time and state. We consider how to achieve Equations (3) and (4) in an online fashion. A schematic description of MCD is shown in Figure 2.

First, we detect change points from data stream (A). Next, we concurrently detect metachanges along time (B) and along state (C). We introduce metachange statistics to quantify these metachanges. Finally, we integrate the metachange statistics along time and state into a statistics (D).

The key challenge of detecting metachanges along time and state is how to describe and integrate them. Our approach describes both metachanges as code lengths with MDL; therefore, it is easy to combine them.

3.1. Detecting Change Points

First, we detect change points $t_1,t_2,\dots$. As our proposed algorithm MCD works in an online fashion, it is necessary for the change detection algorithm to work in an online fashion (e.g., [1,2,3,4,8,9]). In general, MCD is prone to errors by the change detection algorithm and its threshold parameter. We empirically investigate and discuss this point in detail in Section 4.

3.2. Detecting Metachanges along Time

For the detected change points $t_1\dots$, let us consider intervals between the successive change points $I_i=[t_{i-1},t_i-1]$, with length $x_i=t_i-t_{i-1}$. For an interval sequence $x^i=x_1\dots x_i$, we consider how to achieve Equation (3) in an online fashion. We define metachange along time (MCAT) $a_{t_i}$ as the predictive code length

$$a_{t_i}\stackrel{\mathrm{def}}{=}-log{p}_x(x_i;{\widehat{\xi }}_{x^{i-1}}),$$

(5)

where $p_x\in {\mathcal{F}}_{\mathrm{t}}$, ${\mathcal{F}}_{\mathrm{t}}=\left\{p_x(x;\xi )\right\}$ is a parametric class of probability distribution, and ${\widehat{\xi }}_{x^{i-1}}$ is estimated using $x^{i-1}=x_1\dots x_{i-1}$. For example, we can estimate ${\widehat{\xi }}_{x^{i-1}}$ as the maximum likelihood estimator. To deal with nonstationary data streams, we use the online discounting maximum likelihood estimator [24]

$$\begin{array}{c}{\widehat{\xi }}_{x^{i-1}}=\underset{\xi }{argmax}\sum _{t=1}^{i-1}r{(1-r)}^{i-1-t}log{p}_x(x_t;\xi ),\end{array}$$

(6)

where $0<r<1$ is a discounting parameter. An increase in r has a greater effect on forgetting past data.

In this paper, we introduce a parametric class of the exponential distribution

$${\mathcal{F}}_{\mathrm{t}}=\left\{p_x(x;\xi )=\xi exp(-\xi x),\xi >0\right\}.$$

(7)

By substituting Equation (7) into Equation (6), we get

$$\begin{array}{cc}\hfill {\widehat{\xi }}_{x^{i-1}}& =\underset{\xi }{argmax}\sum _{t=1}^{i-1}r{(1-r)}^{i-1-t}log\left(\xi exp(-\xi x_t)\right)\hfill \\ & =\underset{\xi }{argmax}\sum _{t=1}^{i-1}r{(1-r)}^{i-1-t}\left(log\xi -\xi x_t\right).\hfill \end{array}$$

(8)

The inside of argmax in the right-hand side of Equation (8) is expanded as

$$\begin{array}{cc}\hfill \sum _{t=1}^{i-1}r{(1-r)}^{i-1-t}(log\xi -\xi x_t)& =rlog\xi \sum _{t=1}^{i-1}{(1-r)}^{i-1-t}-r\xi \sum _{t=1}^{i-1}{(1-r)}^{i-1-t}{x}_t\hfill \\ \hfill & =rlog\xi \frac{1-{(1-r)}^{i-1}}{r}-r\xi \sum _{t=1}^{i-1}{(1-r)}^{i-1-t}{x}_t\hfill \\ \hfill & =log\xi (1-{(1-r)}^{i-1})-r\xi \sum _{t=1}^{i-1}{(1-r)}^{i-1-t}{x}_t.\hfill \end{array}$$

(9)

The right-hand side of Equation (9) is maximized by deriving it with respect to $\xi$. As a result, we obtain the following optimal solution:

$${\widehat{\xi }}_{x^{i-1}}=\frac{1-{(1-r)}^{i-1}}{r{\sum }_{t=1}^{i-1}{(1-r)}^{i-1-t}{x}_t}.$$

(10)

Thus, by substituting Equation (10) into Equation (5), MCAT at $t_i$ is

$$\begin{array}{cc}\hfill a_{t_i}& =-log{p}_x(x_i;{\widehat{\xi }}_{x^{i-1}})=-log{\widehat{\xi }}_{x^{i-1}}+{\widehat{\xi }}_{x^{i-1}}{x}_i.\hfill \end{array}$$

(11)

In practice, we judge that a metachange occurs along time when MCAT changes greatly between the change points. Technically, we use the change rate of MCAT: a metachange occurs along time when $|(a_{t_i}-a_{t_{i-1}})/a_{t_{i-1}}|>{ϵ}_{\mathrm{t}}$ holds, where ${ϵ}_{\mathrm{t}}>0$ is a threshold parameter. We call the algorithm described above as the metachange detection along time algorithm (MCD-T).

As for computational cost of MCAT, Equation (10) is written as

$${\widehat{\xi }}_{x^{i-1}}=\frac{1-{(1-r)}^{i-1}}{r{s}_{i-1}},$$

where

$$s_{i-1}\stackrel{\mathrm{def}}{=}\sum _{j=1}^{i-1}{(1-r)}^{i-1-j}{x}_j.$$

$s_i$ and $s_{i-1}$ satisfy the following relation:

$$s_i=(1-r)s_{i-1}+x_i.$$

Therefore, the computational cost of MCAT $a_{t_i}$ is $O\left(i\right)$.

Example:

We consider a data stream with a length of 200 time intervals between change points: $x_i=100(i=1,\dots ,100)$ and $x_i=500(i=101,\dots ,200)$. This means that there are 201 change points ${\left\{t_i\right\}}_{i=1}^{201}$. If we assume $t_1=100$, then $t_2=200,\dots ,t_{101}=$ 10,100, $t_{102}=$ 10,600, $t_{103}=$ 11,100, $\dots ,$ $t_{201}=$ 60,100. Then, $x_i$ is calculated as $x_i=t_i-t_{i-1}$. Figure 3 shows the time intervals at the change points (Figure 3, top), MCATs $a_{t_i}$ (Figure 3, second graph), the change rate of MCATs $|(a_{t_i}-a_{t_{i-1}})/a_{t_{i-1}}|$ (Figure 3, third graph), and ${\widehat{\xi }}_{x^{i-1}}$ (Figure 3, bottom). We observe in Figure 3 that we can detect the metachange along time when we choose a suitable threshold ${ϵ}_{\mathrm{t}}$. Here, the discounting parameter is set to $r=0.5$.

3.3. Detecting Metachanges Along State

For a change point $t_i$ detected in Section 3.1, we consider how to achieve Equation (4) in an online fashion. We consider a subset of time around $t_i$ for $\mathrm{Neighbor}\left(t_i\right)$ in Equation (4). The subset is denoted by $J_i=\left[t_i-h,t_i+h\right]$, where $h\in \mathbb{N}$ is a window size. Thus, we consider a sequence $y_{t_i-h}^{t_i+h}=y_{t_i-h}\dots y_{t_i+h}$, with length $n=2h+1$. We introduce a parametric class of probability distributions ${\mathcal{F}}_{\mathrm{s}}=\{p_y(Y;\eta );\eta \in H\}$. Here, Y is a random variable and $\eta$ is a real-valued parameter. H is the associated parameter space.

Next, we define metachange statistics along state (MCAS) at change point $t_i$. First, two statistics, $b_{t_i}^{+}$ and $b_{t_i}^{-}$, are introduced. These are defined as the difference between two code lengths for $y_{t_i+1}^{t_i+h}$: one is the “expected” code length, estimated using the parameter change at $t_{i-1}$ and the estimated parameter with $y_{t_i-h}^{t_i-1}$. The other is the code length with the parameter estimated in terms of $y_{t_i+1}^{t_i+h}$. Formally, $b_{t_i}^{\pm }$ is defined as the difference between the predictive code length and the NML code length [20] after the change point. The former is calculated as the predictive code length, which is the total code length for encoding $y_{t_i+1}^{t_i+h}$ in a predictive way, using the estimated parameter ${\eta }^{\pm }$ as follows:

$$\frac{1}{h}\sum _{t=t_i+1}^{t_i+h}-log{p}_y(y_t;{\widehat{\eta }}^{\pm }),$$

(12)

where ${\widehat{\eta }}^{\pm }$ is defined as

$${\widehat{\eta }}^{\pm }\stackrel{\mathrm{def}}{=}{\widehat{\eta }}_{y_{t_i-h}^{t_i-1}}\pm \left({\widehat{\eta }}_{y_{t_{i-1}+1}^{t_{i-1}+h}}-{\widehat{\eta }}_{y_{t_{i-1}-h}^{t_{i-1}-1}}\right),$$

(13)

which indicates the parameter change to the same side and the opposite side in the same way as the previous change point $t_{i-1}$. Here, ${\widehat{\eta }}_{y_{{\tau }_1}^{{\tau }_2}}$ means the maximum likelihood estimator of $\eta$ using $y_{{\tau }_1}^{{\tau }_2}=y_{{\tau }_1}\dots y_{{\tau }_2}$.

The latter is calculated as the NML code length, which is defined as the negative logarithm of the NML distribution [20]:

$$\frac{1}{h}\left(\sum _{t=t_i+1}^{t_i+h}-log{p}_y(y_t;{\widehat{\eta }}_{y_{t_i+1}^{t_i+h}})+log{C}_h\right).$$

(14)

The difference between Equation (12) and Equation (14) is given by

$$\begin{array}{c}{b}_{t_i}^{\pm }\stackrel{\mathrm{def}}{=}\frac{1}{h}\left\{\sum _{t=t_i+1}^{t_i+h}\left(-log{p}_y(y_t;{\widehat{\eta }}^{\pm })+log{p}_y(y_t;{\widehat{\eta }}_{y_{t_i+1}^{t_i+h}})\right)-log{C}_h\right\},\end{array}$$

(15)

where $C_h={\sum }_{z_{t_i+1}^{t_i+h}}{max}_{\eta }{p}_y(z_{t_i+1}^{t_i+h};\eta )$ in Equation (15) is computed using Rissanen’s approximation formula under some regularity conditions [23]:

$$log{C}_h\approx \frac{k}{2}log\frac{h}{2\pi }+log\int \sqrt{\left|I\right(\theta \left)\right|}\mathrm{d}\theta ,$$

where k is the dimension of H and $I\left(\theta \right)\stackrel{\mathrm{def}}{=}{\mathbb{E}}_{\eta }[-{\partial }^{2}log{p}_y(Y;\eta )/\partial {\eta }_i\partial {\eta }_j]$ is the Fisher information matrix at the parameter value $\eta$. Intuitively, Equation (15) quantifies the redundant code length for coding $y_{t_i+1}^{t_i+h}$ with the parameters estimated in terms of the parameter change at $t_{i-1}$ and the parameter values in the former part of $t_i$.

Finally, we define MCAS as

$$b_{t_i}\stackrel{\mathrm{def}}{=}min(b_{t_i}^{+},b_{t_i}^{-}),$$

(16)

which means that metachanges along state are quantified by the relative magnitude of changes in the parameters in this paper. The computational cost of MCAS is $O\left(h\right)=O\left(1\right)$. We judge that a metachange along state occurs at $t_i$ when $b_{t_i}>{ϵ}_{\mathrm{s}}$ holds, where ${ϵ}_{\mathrm{s}}>0$ is a threshold parameter. We call the algorithm described above as the metachange detection along state algorithm (MCD-S).

Example:

We generate a data stream with length 11,250:

$$y_t\sim \left\{\begin{array}{c}\mathcal{N}(0.0,0.05)(t=1,\dots ,1000)\hfill \\ \mathcal{N}(1.0,0.05)(t=1001,\dots ,2000)\hfill \\ \mathcal{N}(0.0,0.05)(t=2001,\dots ,3000)\hfill \\ \mathcal{N}(1.0,0.05)(t=3001,\dots ,4000)\hfill \\ \mathcal{N}\left(\right(t-4001)/1000,0.05)(t=4001,\dots ,5000)\hfill \\ \mathcal{N}(0.0,0.05)(t=5001,\dots ,6000)\hfill \\ \mathcal{N}\left(\right(t-6000)/1000,0.05)(t=6001,\dots ,7000)\hfill \\ \mathcal{N}(1.0,0.05)(t=7001,\dots ,8000)\hfill \\ \mathcal{N}(1-(t-8000)/250,0.05)(t=8001,\dots ,8250)\hfill \\ \mathcal{N}(0.0,0.1)(t=8251,\dots ,9250)\hfill \\ \mathcal{N}(1.0,0.1)(t=9251,\dots ,10,250)\hfill \\ \mathcal{N}(1.0,0.3)(t=10,251,\dots ,11,250)\hfill \end{array}\right.,$$

where $\mathcal{N}(\mu ,\sigma )$ denotes the probability density function of the univariate normal distribution with mean $\mu$ and standard deviation $\sigma$.

Figure 4 shows data stream $\left\{y_t\right\}$ (Figure 4, top) and statistics $\left\{b_{t_i}\right\}$ (Figure 4, bottom). The parameter is set to $h=200$. True change points occur at 1001, 2001, 3001, 4001, 5101, 6001, 7001, 8001, 8251, 9251, and 10,251. Figure 4 shows that the statistics $b_{t_i}$ increase when there is a change in how parameters behave around a change point between successive change points. At $t_2=2001$ and $t_3=3001$, $b_{t_i}$ are relatively small, which shows that parameter changes (i.e., their magnitudes) do not differ much between $t_1=1001$ and $t_2=2001$ and between $t_2=2001$ and $t_3=3001$. However, $b_{t_i}$ increases at $t_4=4001$ because the change shifts to a gradual change from an abrupt one. These results indicate that MCAS provides information regarding changes in the behavior around the change points.

3.4. Integrating Metachange Statistics

Finally, we consider how to integrate MCAT $a_{t_i}$ and MCAS $b_{t_i}$ at a change point $t_i$. Because $a_{t_i}$ and $b_{t_i}$ are code lengths, they can be summed. Therefore, we propose adding $a_{t_i}$ and $b_{t_i}$ with weighting. Integrated metachange (MCI) $s_{t_i}$ at $t_i$ is defined as

$$s_{t_i}\stackrel{\mathrm{def}}{=}{a}_{t_i}+\lambda b_{t_i},$$

(17)

where $\lambda \in \mathbb{R}$ is a hyperparameter. We should carefully choose $\lambda$ with data. In Section 4.3, $\lambda$ is determined using a grid search.

In practice, we judge that a metachange along both time and state occur at $t_i$ when MCI greatly changes between the change points. As in the case of metachanges along time in Section 3.2, we use the change rate of MCI: a metachange along both time and state occurs at $t_i$ if $|(s_{t_i}-s_{t_{i-1}})/s_{t_{i-1}}|>{ϵ}_{\mathrm{ts}}$, where ${ϵ}_{\mathrm{ts}}>0$ is a threshold parameter.

We call the overall algorithm described above MCD; it is summarized in Algorithm 1.

Algorithm 1 MCD.

Input:r: discounting parameter ($0<r<1$), h: window size, ${ϵ}_{\mathrm{ts}}$: threshold parameter Output:$a_{t_i}$: metachange statistics along time, $b_{t_i}$: metachange statistics along state, $s_{t_i}$: integrated metachange statistics. 1: $i=1$ 2: for$t=1,\dots$do 3:    Input $y_t$. 4:    Detect change point with a change detection algorithm. 5:    if t is change point then 6:         $t_i\leftarrow t$. 7:         $x_i\leftarrow t_i-t_{i-1}$. 8:         Calculate metachange statistics along time $a_{t_i}$ according to Equation (11). 9:         Calculate metachange statistics along state $b_{t_i}$ according to Equation (16). 10:      Calculate integrated metachange statistics $s_{t_i}$ according to Equation (17). 11:       Raise an alarm if and only if $|(s_{t_i}-s_{t_{i-1}})/s_{t_{i-1}}|>{ϵ}_{\mathrm{ts}}$. 12:       $i\leftarrow i+1$. 13:    end if 14: end for

4. Experiment

We conducted five experiments to confirm the effectiveness of the proposed algorithm MCD (https://github.com/s-fuku/metachange).

4.1. Synthetic Dataset 1 (Metachanges along Time)

We defined six levels of time intervals between change points referring to the work in [13,25]. The interval lengths were 100,000, 50,000, 10,000, 5000, 1000, and 500. The change points were set using a Bernoulli distribution oscillating between $\mu =0.2$ and $\mu =0.8$. For each combination of two intervals, we generated the streams based on the scheme above. Each stream contained 100 change points. In what follows, $L_1$ and $L_2$ indicate the first and second interval lengths, respectively.

We confirmed the effectiveness of MCD by comparing it with a volatility detector (VD) [13]. We used the SEED algorithm [13] and the sequential MDL-change statistics algorithm (SMDL) [8] for change detection. SEED was based on ADWIN2 [21] and its parameters were set to $\delta =0.05,\mathsf{\Gamma }=75,\widehat{ϵ}=0.025$, and $\alpha =0.025$, which are the same as those in [13]. The window size w of SMDL was set to $w=0.2L_1$, and the threshold parameter $ϵ$ was set to $ϵ=0.01$. For the Bernoulli distribution, the change score ${\mathsf{\Psi }}_t$ of SMDL at time t was calculated as

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_t& =-{\widehat{\mu }}_{0}log{\widehat{\mu }}_0-(1-{\widehat{\mu }}_0)log(1-{\widehat{\mu }}_0)\hfill \\ \hfill & -\frac{1}{2}\left(-{\widehat{\mu }}_{1}log{\widehat{\mu }}_1-(1-{\widehat{\mu }}_1)log(1-{\widehat{\mu }}_1)\right)-\frac{1}{2}\left(-{\widehat{\mu }}_{2}log{\widehat{\mu }}_2-(1-{\widehat{\mu }}_2)log(1-{\widehat{\mu }}_2)\right),\hfill \end{array}$$

where ${\widehat{\mu }}_0={\sum }_{i=t-w}^{t+w}{y}_i/(2w+1)$, ${\widehat{\mu }}_1={\sum }_{i=t-w}^{t-1}{y}_i/w$, and ${\widehat{\mu }}_2={\sum }_{i=t}^{t+w}{y}_i/(w+1)$. If ${\mathsf{\Psi }}_t>ϵ$, t is regarded as a change point. We determined that t was a change point if the change score ${\mathsf{\Psi }}_t$ was the maximum. The parameter of MCD-T was set to $r=0.2$. Below, we discuss the dependency of MCD-T on r in Figure 5. For VD, buffer size $B=32$ and reservoir size $R=32$, which were the same as in [13]. We also discuss the dependency of VD on B and R below in Figure 6. In running SEED [13], we used the Java source code provided by the authors (https://www.cs.auckland.ac.nz/research/groups/kmg/DavidHuang.html). We started to use change points when its number reached $B+R$ for MCD-T and VD because the buffer and the reservoir of VD are not full until $B+R$ intervals arrive.

We investigated the trade-off between detection delay and accuracy in terms of benefit and false alarm rate, defined as in [8,26]. For MCD-T, we first fixed the threshold parameter ${ϵ}_{\mathrm{t}}$ and converted MCAT $\left\{a_{t_i}\right\}$ in Equation (11) to binary alarms $\left\{{\alpha }_{t_i}\right\}$. That is, ${\alpha }_{t_i}=𝟙\left(\right|(a_{t_i}-a_{t_{i-1}})/a_{t_{i-1}}|>{ϵ}_{\mathrm{t}})$, where $𝟙\left(t\right)$ denotes the binary function that takes 1 if and only if t is true. We evaluated MCD-T by varying ${ϵ}_{\mathrm{t}}$. We let $\tau$ be a maximum tolerant delay of metachange detection. When the metachange really started from $t^{*}$, we defined the benefit of an alarm at time t as

$$b(t;t^{*})=\left\{\begin{array}{cc}1-\frac{|t-t^{*}|}{\tau }\hfill & (0\le |t-t^{*}|<\tau ),\hfill \\ 0\hfill & \left(\mathrm{otherwise}\right).\hfill \end{array}\right.$$

(18)

The number of false alarms was calculated as

$$n\left({\alpha }_1^m\right)\stackrel{\mathrm{def}}{=}\sum _{k=1}^m{\alpha }_{t_k}𝟙(b(t_k,t^{*})=0).$$

(19)

We visualized the performance by plotting the recall rate of the total benefit, b, against the false alarm rate, $n/{sup}_{{ϵ}_{\mathrm{t}}}n$, with ${ϵ}_{\mathrm{t}}$ varying. Likewise, for VD, ${\alpha }_{t_i}$ was calculated using the relative volatility between the variances of the buffer and the reservoir by varying the threshold parameter $\beta$. We evaluated all four combinations of change detectors SEED and SMDL and metachange detectors MCD-T and VD by calculating the average and standard deviation of the area under the curve (AUC) of the benefit vs. FAR curves. The AUC scores were calculated over 50 sequences. The delay parameter was set to $\tau =5L_2$.

Table 1. Average area under the curve (AUC) scores on Synthetic Dataset 1 ($r=0.2,\tau =5L_2$). Boldfaces describe best performances.

${\mathit{L}}_1$	${\mathit{L}}_2$	SEED		SMDL
		MCD-T	VD	MCD-T	VD
100,000	50,000	$\mathbf{0.603}\pm \mathbf{0.180}$	$0.458\pm 0.199$	$0.500\pm 0.060$	$0.313\pm 0.160$
100,000	10,000	$0.621\pm 0.147$	$0.310\pm 0.167$	$\mathbf{0.710}\pm \mathbf{0.254}$	$0.463\pm 0.113$
100,000	5000	$0.645\pm 0.129$	$0.328\pm 0.152$	$\mathbf{0.668}\pm \mathbf{0.223}$	$0.416\pm 0.164$
100,000	1000	$\mathbf{0.651}\pm \mathbf{0.110}$	$0.275\pm 0.135$	$0.512\pm 0.123$	$0.448\pm 0.107$
100,000	500	$\mathbf{0.697}\pm \mathbf{0.140}$	$0.336\pm 0.140$	$0.660\pm 0.111$	$0.506\pm 0.138$
50,000	100,000	$\mathbf{0.788}\pm \mathbf{0.093}$	$0.647\pm 0.107$	$0.729\pm 0.067$	$0.639\pm 0.107$
50,000	10,000	$\mathbf{0.671}\pm \mathbf{0.103}$	$0.280\pm 0.130$	$0.605\pm 0.171$	$0.556\pm 0.060$
50,000	5000	$\mathbf{0.708}\pm \mathbf{0.087}$	$0.293\pm 0.144$	$0.617\pm 0.183$	$0.546\pm 0.146$
50,000	1000	$\mathbf{0.718}\pm \mathbf{0.067}$	$0.294\pm 0.140$	$0.655\pm 0.161$	$0.501\pm 0.144$
50,000	500	$\mathbf{0.767}\pm \mathbf{0.110}$	$0.316\pm 0.133$	$0.686\pm 0.074$	$0.470\pm 0.157$
10,000	100,000	$0.863\pm 0.059$	$0.794\pm 0.058$	$\mathbf{0.877}\pm \mathbf{0.068}$	$0.791\pm 0.015$
10,000	50,000	$0.834\pm 0.050$	$0.735\pm 0.050$	$\mathbf{0.876}\pm \mathbf{0.066}$	$0.823\pm 0.026$
10,000	5000	$\mathbf{0.723}\pm \mathbf{0.040}$	$0.344\pm 0.250$	$0.658\pm 0.159$	$0.498\pm 0.084$
10,000	1000	$\mathbf{0.781}\pm \mathbf{0.014}$	$0.375\pm 0.260$	$0.689\pm 0.083$	$0.444\pm 0.077$
10,000	500	$\mathbf{0.809}\pm \mathbf{0.063}$	$0.391\pm 0.256$	$0.671\pm 0.163$	$0.520\pm 0.070$
5000	100,000	$\mathbf{0.856}\pm \mathbf{0.060}$	$0.796\pm 0.067$	$0.854\pm 0.071$	$0.798\pm 0.036$
5000	50,000	$0.825\pm 0.047$	$0.726\pm 0.062$	$\mathbf{0.875}\pm \mathbf{0.032}$	$0.708\pm 0.043$
5000	10,000	$\mathbf{0.777}\pm \mathbf{0.030}$	$0.575\pm 0.139$	$0.716\pm 0.060$	$0.630\pm 0.031$
5000	1000	$\mathbf{0.783}\pm \mathbf{0.009}$	$0.436\pm 0.257$	$0.709\pm 0.009$	$0.353\pm 0.098$
5000	500	$0.816\pm 0.054$	$0.493\pm 0.269$	$\mathbf{0.839}\pm \mathbf{0.191}$	$0.413\pm 0.097$
1000	100,000	$\mathbf{0.872}\pm \mathbf{0.072}$	$0.814\pm 0.072$	$0.836\pm 0.036$	$0.812\pm 0.036$
1000	50,000	$0.844\pm 0.059$	$0.754\pm 0.061$	$\mathbf{0.947}\pm \mathbf{0.037}$	$0.810\pm 0.027$
1000	10,000	$0.802\pm 0.022$	$0.668\pm 0.050$	$\mathbf{0.873}\pm \mathbf{0.049}$	$0.805\pm 0.023$
1000	5000	$0.801\pm 0.014$	$0.648\pm 0.064$	$\mathbf{0.895}\pm \mathbf{0.053}$	$0.812\pm 0.053$
1000	500	$\mathbf{0.816}\pm \mathbf{0.048}$	$0.560\pm 0.242$	$0.711\pm 0.141$	$0.409\pm 0.108$
500	100,000	$\mathbf{0.876}\pm \mathbf{0.068}$	$0.831\pm 0.063$	$0.830\pm 0.079$	$0.820\pm 0.023$
500	50,000	$\mathbf{0.845}\pm \mathbf{0.062}$	$0.767\pm 0.062$	$0.836\pm 0.044$	$0.818\pm 0.010$
500	10,000	$0.827\pm 0.051$	$0.676\pm 0.047$	$\mathbf{0.872}\pm \mathbf{0.023}$	$0.822\pm 0.016$
500	5000	$0.830\pm 0.047$	$0.663\pm 0.042$	$\mathbf{0.864}\pm \mathbf{0.047}$	$0.819\pm 0.017$
500	1000	$0.830\pm 0.050$	$0.612\pm 0.100$	$\mathbf{0.935}\pm \mathbf{0.022}$	$0.853\pm 0.095$

shows the average AUC scores. Table 1 shows that MCD-T with SEED or MCD-T with SMDL outperforms VD with SEED or VD with SMDL. This indicates the effectiveness of MCD-T.

Because MCD-T depends on discounting parameter r and the change detection algorithm used, we investigated these effects. First, we examined the dependency of AUC on r for all combinations of $L_1$ and $L_2$. We calculated AUC for 30 times with $r=0.01,0.05,0.1,0.2,0.3,0.4,$ and $0.5$. We used SEED [13] as the change detection algorithm, and its parameters were set to the same values as above. The dataset used was also the same as in the previous experiment. Figure 5 shows that, when $L_1$ is relatively small (e.g., $L_1=500,1000,5000,10,000$), AUC is not heavily dependent on r. When $L_1$ is larger, however, we observe that the larger r is, the smaller AUC is. This is because, with an increase of $L_1$, the number of false alarms of SEED also increases. In such situations, MCD-T is more prone to the false alarms when r is larger.

Figure 6 shows the dependency of AUC of VD on the buffer size B and the reservoir size R ($B=R$) for comparison. We calculated AUC for 50 times. We observe from Figure 6 that AUC decreases as B increases. In addition, we also see that MCD-T outperforms VD for various combinations of r and $B(=R)$ by comparing Figure 5 with Figure 6.

Next, we investigated the effect of the change detection algorithm used. We used SEED by changing the parameter $\widehat{ϵ}=0.0025,0.005$, and $0.0075$. Other conditions and the dataset were the same as in the previous experiment. Here, $\widehat{ϵ}$ is a hyperparameter that controls the threshold parameter [13]. Figure 7 shows that AUC does not heavily depend on $\widehat{ϵ}$ for all combinations of $L_1$ and $L_2$. In general, the threshold parameter of the change detection algorithm controls the performance of MCD-T. Hence, it should be carefully set.

4.2. Synthetic Dataset 2 (Metachanges along State)

We generated a data stream with length $24L$, where $L=500,1000,2000$. The generated data stream contained a metachange along state. In the former part, each datum was drawn from

$$y_t\sim \left\{\begin{array}{cc}\mathcal{N}(0,0.1)\hfill & (t=1,\dots ,L),\hfill \\ \mathcal{N}(0.5,0.1)\hfill & (t=L+1,\dots ,2L).\hfill \end{array}\right.$$

(20)

After we repeated the procedure 10 times, we obtained a subsequence with length $20L$. In the latter part, each datum was drawn from

$$y_t\sim \left\{\begin{array}{cc}\mathcal{N}\left(\right(t-20L)/2L,0.1)\hfill & (t=20L+1,\dots ,21L),\hfill \\ \mathcal{N}(0,0.1)\hfill & (t=21L+1,\dots ,24L).\hfill \end{array}\right.$$

A metachange along state occurred at $t=20L+1$. For change detection, we employed four algorithms for comparison: (1) SMDL [8], a semi-instant method with the MDL change statistics; (2) ChangeFinder (CF) [1,2,4], a state-of-the-art method of abrupt change detection; (3) Bayesian online change point detection (BOCPD) [3], a retrospective online change point detection with a Bayesian scheme; and (4) ADWIN2 [21], adaptive windowing methods. As we assumed a situation where change and metachange mechanisms do not vary significantly, we decided to choose the best combinations of parameters of each change detection algorithm by grid search, as in [8,27]. We generated 10 sequences with the scheme above and calculated the F-scores for each combination of the following parameters:

• SMDL: Window size $w=50,100$ ($L=500$), $w=100,200$ ($L=1000$), $w=200,400$ ($L=2000$). Threshold parameter $ϵ=0.1,0.2,0.3,0.4,0.5,0.6,0.7$.
• CF: Discounting rate $r=0.003,0.005,0.01,0.03,0.1$. Threshold parameter $\delta =0,0.5,1.0,1.5,2.0$ (regression orders $k_1,k_2=3$, smoothing parameters $T_1,T_2=5$).
• BOCPD: Parameter related to change intervals $\alpha =100,300,600$. Threshold parameter $ϵ=0.1,0.3$.
• ADWIN2: Confidence parameter $\delta =0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9$.

F-score is defined as the harmonic mean of precision and recall, which are calculated using the number of true positives (TP), false positives (FP), and false negatives (FN) as follows [9]: TP is the number of true change points that are $\tau$-neighbors of estimated change points. Thus, FP and FN are calculated as $\mathit{FP}=\ell -\mathit{TP}$ and $\mathit{FN}=m-\mathit{TP}$, where ℓ and m are calculated as $\mathit{FP}=\ell -\mathit{TP}$ and $\mathit{FN}=m-\mathit{TP}$, where ℓ and m denotes the total number of estimated and true change points, respectively. Finally, we calculated $\mathit{recall}=\mathit{TP}/(\mathit{TP}+\mathit{FN})$ and $\mathit{precision}=\mathit{TP}/(\mathit{TP}+\mathit{FP})$ for each method. In this experiment, we set $\tau$ to 100.

After optimizing the parameters of each change detection algorithm, we generated 30 data streams with the scheme above and detected change points and the metachange. In the metachange detection, we compared MCD-S with SMDL. We chose SMDL for comparison because it calculates a change score at each time based on changes of parameters with MDL. Hence, a change rate of scores between change points is regarded as the degree of metachange along state. Hereafter, we refer to SMDL for metachange detection as SMDL metachange (SMDL-MC) and the window parameter as $w_{\mathrm{mc}}$. We calculated MCAS in Equation (16) for MCD-S and the change rate $|({\mathsf{\Psi }}_{t_i}-{\mathsf{\Psi }}_{t_{i-1}})/{\mathsf{\Psi }}_{t_{i-1}}|$ for SMDL-MC. ${\mathsf{\Psi }}_t$ is the change score at time t for a univariate normal distribution [8]:

$${\mathsf{\Psi }}_t=\frac{1}{2}log\frac{{\widehat{\sigma }}_0^2}{{\widehat{\sigma }}_1{\widehat{\sigma }}_2}+log\frac{C_{2w_{\mathrm{mc}}}}{C_{w_{\mathrm{mc}}}^2},$$

where ${\widehat{\sigma }}_0$, ${\widehat{\sigma }}_1$, and ${\widehat{\sigma }}_2$ are the maximum likelihood estimators of standard deviations calculated for $y_{t-w_{\mathrm{mc}}+1}^{t+w_{\mathrm{mc}}}$, $y_{t-w_{\mathrm{mc}}+1}^{t-1}$ and $y_t^{t+w_{\mathrm{mc}}}$, respectively. $C_k$ is the normalizer of the normalized maximum likelihood code length [20]

$$log{C}_k=\frac{1}{2}log\frac{16{\mu }_{max}}{\pi {\sigma }_{min}^2}+\frac{k}{2}log\frac{k}{2e}-log\mathsf{\Gamma }\left(\frac{k-1}{2}\right),$$

where $\mathsf{\Gamma }$ is the gamma function. In this paper, ${\mu }_{max}=2$ and ${\sigma }_{min}=0.005$. The window parameters h of MCD-S and $w_{\mathrm{mc}}$ of SMDL-MC were set to $h,w_{\mathrm{mc}}=100$ ($L=500$), $h,w_{\mathrm{mc}}=200$ ($L=1000$), and $h,w_{\mathrm{mc}}=400$ ($L=2000$). In calculating the F-scores, the maximum tolerant delay was set to $\tau =0.5L$.

Table 2. Average AUC scores on Synthetic Dataset 2. The first and second headers represent change detection and metachange detection algorithms, respectively. Boldfaces describe best performances.

L	SMDL		CF		BOCPD		ADWIN2
	MCD-S	SMDL-MC	MCD-S	SMDL-MC	MCD-S	SMDL-MC	MCD-S	SMDL-MC
500	$\mathbf{0.887}\pm \mathbf{0.100}$	$0.795\pm 0.156$	$0.874\pm 0.111$	$0.851\pm 0.170$	$0.701\pm 0.318$	$0.572\pm 0.332$	$0.797\pm 0.186$	$0.853\pm 0.114$
1000	$\mathbf{0.921}\pm \mathbf{0.018}$	$0.905\pm 0.012$	$0.912\pm 0.042$	$0.830\pm 0.052$	$0.751\pm 0.323$	$0.743\pm 0.291$	$0.834\pm 0.094$	$0.847\pm 0.048$
2000	$\mathbf{0.970}\pm \mathbf{0.010}$	$0.953\pm 0.011$	$0.912\pm 0.033$	$0.843\pm 0.022$	$0.829\pm 0.124$	$0.821\pm 0.138$	$0.951\pm 0.032$	$0.887\pm 0.046$

shows the average AUC values of MCD-S and SMDL-MC for the detection of metachanges along state at $t=20L+1$. The first and second rows in the header represent change detection and metachange detection algorithms, respectively. The best parameters for each combination of change detection and metachange detection algorithms are $ϵ=0.7,w=100$ ($L=500$), $ϵ=0.7,w=200$ ($L=1000$), and $ϵ=0.7,w=400$ ($L=2000$). Table 2 shows that MCD-S outperforms SMDL-MC overall because MCD-S deals with metachanges along state directly in terms of MCAS, whereas SMDL-MC only quantifies the difference in code lengths between situations where there is a change and where there is no change.

We further investigated the effects of window size h and threshold parameters of the change detection algorithms. We chose SMDL [8] for change detection. Figure 8 shows the dependency of AUC on h and threshold parameter $ϵ$ of SMDL. The interval length was set to $L=500$, threshold parameter was set to $ϵ=0.1,0.2,0.3,0.4,0.5,0.6,0.7$, and $h=w=50,100,150$, where w is the window parameter of SMDL. Figure 8 (top and bottom) shows the dependency of AUC of MCD-S on the threshold parameter $ϵ$ of SMDL and the dependency of F-score of SMDL on $ϵ$, respectively. We observe in Figure 8 (top) that AUC of MCD-S decreases between $ϵ=0.2$ and $0.4$, but, when $ϵ$ exceeds $0.4$, AUC begins to increase for $h=50,100,150$. This reflects the fact that there are many local maximum points of the change scores of SMDL, leading to false alarms of change points around $ϵ=0.2$–$0.4$. It is noticeable that F-scores of SMDL decrease for $ϵ=0.1$$(h=100)$, and for $ϵ=0.2$$(h=150)$, but AUCs of MCD-S do not do so much. This is because SMDL detects many false positive change points, but it detects the metachange point accurately.

As for the dependency of AUC on window size h, we observe that AUC generally increases as h increases for the same $ϵ$.

4.3. Synthetic Dataset 3 (Metachanges Along Time and State)

We generated a data stream that contained metachanges along both time and state. The stream consisted of two subsequences. The former part repeated changes of mean. Each instance was drawn from Equation (20) with $L=L_1$. We repeated the procedure for 50 times and obtained a subsequence with length $100L_1$. The latter part comprised the following four parts, each with length $L_2$:

$$y_t\sim \left\{\begin{array}{cc}\mathcal{N}(0,0.1)\hfill & (t=100L_1+1,\dots ,100L_1+L_2),\hfill \\ \mathcal{N}(0.45,0.1)\hfill & (t=100L_1+L_2+1,\dots ,100L_1+4L_2).\hfill \end{array}\right.$$

In total, we obtained a data stream with length $100L_1+4L_2$. A metachange along both time and state occurred at $t=100L_1+L_2+1$. We chose lengths $L_1$ and $L_2$ among 400, 450, and 500.

We detected the metachange in the following three ways: we first detected change points with the same algorithms as in Section 4.2, and then detected the metachanges with MCD-T, MCD-S, and MCD. The parameters of the change detection algorithms were tuned as in Section 4.2. The ranges of parameters were the same as those in Section 4.2. except that, for SMDL, the threshold parameter $ϵ=0.05,0.1,0.15$ for all combinations of $L_1$ and $L_2$. The parameter of MCD-T was selected among $r=0.1,0.2,0.3$ and MCD-S was among $h=0.1L_1,0.2L_1$. The window size of SMDL were selected among $w=h$, and the maximum tolerant delay was $\tau =L_2$. We chose the weight parameter $\lambda$ in Equation (17) among $\lambda =0.001,0.01,0.1,1,5,10$. For VD, the buffer and reservoir sizes (B and R) were selected among $16,24,32$. All the parameters were selected with grid search for the AUCs of metachange detection to be maximum.

Table 3 shows the average AUC values. Table 3a–c show average AUC values with MCD-T, MCD-S, and MCD. Table 3a shows that MCD combined with SMDL as the change detection algorithm outperforms MCD-S and MCD-T.

Table 4. Best parameters for each combination of intervals.

${\mathit{L}}_1$	${\mathit{L}}_2$	r	w	h	$\mathit{\lambda }$
400	450	$0.2$	$0.2L_1$	$0.2L_1$	$0.1$
400	500	$0.3$	$0.2L_1$	$0.2L_1$	$0.01$
450	400	$0.1$	$0.2L_1$	$0.2L_1$	$0.1$
450	500	$0.2$	$0.2L_1$	$0.2L_1$	$0.1$
500	400	$0.3$	$0.2L_1$	$0.2L_1$	$0.01$
500	450	$0.1$	$0.2L_1$	$0.2L_1$	$0.1$

shows the best parameters for each combination of intervals. We observe that the more intensive a metachange along time is, the bigger r is and the less $\lambda$ becomes. These results reflect the fact that it is necessary to adapt to recent data, and MCAT increases in such a situation, leading to the decrease of $\lambda$.

4.4. Real Dataset: Human Action Recognition Data

We applied MCD to the detection of metachanges in human action recognition data called HASC-PAC2016 dataset [28] (HASC-PAC2016 dataset is publicly available at http://hub.hasc.jp/). The data were collected from the Human Activity Sensing Consortium (HASC, http://hasc.jp/). HASC-PAC2016 dataset contains sequences of acceleration data for three axes, and each sequence is segmented into one of six action labels: “stay”, “walk”, “jog”, “skip”, “stair up”, (go upstairs) and “stair down” (go downstairs). For this experiment, we aimed to evaluate the effectiveness of our proposed algorithm MCD by using a data stream with ground truth of “changes of action changes” and “changes of intervals of actions”. The former corresponds to metachanges along state, and the latter to metachanges along time. We combined each action into a data stream as follows: first, we repeated “stay” and “walk” alternately for 15 times; then “jog” and “skip” for 15 times; and, finally, “stair up” and “stair down” for 15 times. We repeated each pair of actions for 15 times because “stair up” and “stair down” have only 15 files, which are the fewest in all the six actions. We obtained a data stream of length 89,324.

Table 5. Files for generating a sequence of Person06023.

Action Label	Files
stay	HASC N-acc.csv (N = 0605581–0605595)
walk	HASC N-acc.csv (N = 0608420–0608434)
jog	HASC N-acc.csv (N = 0611173–0611187)
skip	HASC N-acc.csv (N = 0613411–0613425)
stair up	HASC N-acc.csv (N = 0615620–0615634)
stair down	HASC N-acc.csv (N = 0614162–0614166)

shows the files used for a participant named Person06023. We read the files sequentially in alphabetical order for each action. Figure 9 shows the data stream we obtained. Here, acc_X, acc_Y, and acc_Z represent accelerations for x-, y-, and z-axes, respectively.

First, we detected change points with SMDL [8]. It was a challenge to determine the hyperparameters of SMDL—window size w and threshold parameter $ϵ$—in an online change detection. We tuned w and $ϵ$ with the remaining dataset for Person06023, which alternated “stay” and “walk” four times, and “jog” and “skip” likewise. Although this dataset lacked “stair up” and “stair down”, we thought that it was enough to estimate the best configuration of w and $ϵ$. We calculated F-score as described in Section 4.2 for the change points between different action labels. We selected $w=900$ and $ϵ=0.75$ among $w\in \{500,600,700,800,900,1000\}$ and $ϵ\in \{0,0.25,0.5,0.75,1\}$. Figure 10 shows histograms of intervals for each action label. We observe in Figure 10 that most of the intervals are around 960–970 for “jog”, “walk”, and “skip”, whereas, for “stay”, “stair up”, and “stair down”, the intervals are around 1020. We can see that $w=900$ was enough to detect changes.

We applied SMDL to the stream and obtained the estimated change scores $\left\{{\mathsf{\Psi }}_t\right\}$ at each time point. We calculated ${\mathsf{\Psi }}_t$ with the multivariate normal distribution. Specifically, ${\mathsf{\Psi }}_t$ is calculated as

$$\begin{array}{cc}\hfill {\mathsf{\Psi }}_t& =\frac{1}{2}log\frac{|{\widehat{\Sigma }}_0{|}^2}{|{\widehat{\Sigma }}_1\left|\right|{\widehat{\Sigma }}_2|}+\frac{1}{2w}log\frac{C_{2w}}{C_w^2}\hfill \\ \hfill & +\frac{1}{2w}\left\{\sum _{i=t-w}^{t+w}{(y_i-{\widehat{\mu }}_0)}^{\top }{\widehat{\Sigma }}_0^{-1}(y_i-{\widehat{\mu }}_0)\right.\hfill \\ \hfill & \left.-\sum _{i=t-w}^{t-1}{(y_i-{\widehat{\mu }}_1)}^{\top }{\widehat{\Sigma }}_1^{-1}(y_i-{\widehat{\mu }}_1)-\sum _{i=t}^{t+w}{(y_i-{\widehat{\mu }}_2)}^{\top }{\widehat{\Sigma }}_2^{-1}(y_i-{\widehat{\mu }}_2)\right\},\hfill \end{array}$$

(21)

where ${\widehat{\mu }}_0=1/(2w+1){\sum }_{i=t-w}^{t+w}{y}_i$, ${\widehat{\mu }}_1=1/w{\sum }_{i=t-w}^{t-1}{y}_i$, and ${\widehat{\mu }}_2=1/(w+1){\sum }_{i=t+1}^{t+w}{y}_i$. ${\widehat{\Sigma }}_0=\frac{1}{2w}{\sum }_{i=t-w}^{t+w}(y_i-{\widehat{\mu }}_0){(y_i-{\widehat{\mu }}_0)}^{\top }$, ${\widehat{\Sigma }}_1=\frac{1}{w}{\sum }_{i=t-w}^{t-1}(y_i-{\widehat{\mu }}_1){(y_i-{\widehat{\mu }}_1)}^{\top }$, and ${\widehat{\Sigma }}_1=\frac{1}{w+1}{\sum }_{i=t}^{t+w}(y_i-{\widehat{\mu }}_2){(y_i-{\widehat{\mu }}_2)}^{\top }$.

Note that $C_w$ in Equation (21) is the normalizer of the NML code length [29,30]:

$$\begin{array}{cc}\hfill log{C}_w& =-(m+1)log\frac{m}{2}+\frac{m}{2}log{\mu }_{max}-\frac{m^2}{2}log{\sigma }_{min}+\frac{mw}{2}log\frac{w}{2e}\hfill \\ \hfill & -log\mathsf{\Gamma }\left(\frac{m}{2}\right)-log{\mathsf{\Gamma }}_m\left(\frac{w-1}{2}\right),\hfill \end{array}$$

(22)

where m is the dimension of the data stream, $\mathsf{\Gamma }$ is the gamma function, and ${\mathsf{\Gamma }}_m$ is calculated as

$${\mathsf{\Gamma }}_m\left(x\right)={\pi }^{\frac{m(m-1)}{4}}\prod _{j=1}^m\mathsf{\Gamma }\left(x+\frac{1-j}{2}\right).$$

We set ${\mu }_{max}=50$ and ${\sigma }_{min}=0.005$.

Next, we defined the ground truths for metachanges along state at two time points where the changes of action label changes occurred: $t=$ 29,752 from “jog” to “stair up”, and $t=$ 59,588 from “walk” to “skip”. Moreover, we also defined the ground truths for metachanges along time at time points where the changes of intervals occurred. We see in Figure 10 that the distributions are significantly different between four types of “changes of action changes”: from “stay” to “jog”, from “jog” to “stair up”, from “stair up” to “walk”, and from “skip” to “stair down”.

We detected metachanges along time with MCD-T and volatility detector (VD) [13], and compared them. Figure 11 shows the estimated MCAT with MCD-T and the relative volatility with VD. The parameter of MCD-T was set to $r=0.1,0.2,0.3$, whereas one of VD was $B=R=10,15,20$. Figure 11 shows the results.

We observe in Figure 11 that MCAT detects the metachanges along time between the four action pairs, respectively, for $r=0.1,0.2,$ and $0.3$. However, the relative volatility fails to detect some of these metachanges along time.

We detected metachanges along state with MCD-S and the change rate of the MDL change statistics [8]. Figure 12 shows the estimated MCAS with MCD-S and the MDL change statistics. We observe in Figure 12 that both MCD-S and the MDL change statistics detect a time point around $t=$ 29,752 from “jog” to “stair up”. However, the MDL change statistics do not change significantly at a time point around $t=$ 59,588, where a metachange along state happened from “walk” to “skip”. It indicates that the change rate of the MDL change statistics failed to detect the metachange along state around $t=$ 59,588, whereas MCD-S detected it successfully.

In summary, the proposed algorithm MCD detected metachanges along both time and state more accurately than other methods.

4.5. Real Dataset: Production Condition Data

We applied MCD to the detection of metachanges in the production condition data. The data were collected from a factory of a manufacturing company. Each datum comprised eight attributes, and the length of the stream was 26,450. The factory reported that important events occurred 10 times during the study period, at $t=$ 668, 2634, 2635, 9663, 13,230, 13,231, 17,372, 17,832, 20,131, and 25,441. Figure 13 shows the attributes from the stream. The dashed line indicates the time points where important events occurred. We investigated whether the detected metachanges were signs of important events, and we finally concluded that it might be true. The details are as follows.

Figure 13 shows that the scales of attributes were different. Hence, we normalized each attribute X to $(X-\mu )/\sigma$, where $\mu$ and $\sigma$ are the sample mean and standard deviation, respectively, which were calculated with the first 250 time points. First, we applied SMDL [8] to the stream and obtained the estimated change scores $\left\{{\mathsf{\Psi }}_t\right\}$ at each time. We calculated ${\mathsf{\Psi }}_t$ with the multivariate normal distribution in Equation (21). The window sizes w of SMDL and h of MCD were set to $w=h=250$ by field knowledge that it roughly represents a unit of production. Moreover, ${\mu }_{max}$ and ${\sigma }_{min}$ in Equation (22) were set to 60 and $0.001$, respectively. Next, we detected change points $t_1,t_2,\dots$ as time points where the change scores ${\mathsf{\Psi }}_{t_i}$ were locally maximum within an interval where ${\mathsf{\Psi }}_t>ϵ$. We set $ϵ=0.3$ when the total change points detected was less than 0.5% of the total length. It is a business demand by a factory, and so there were not many alarms. The number of detected change points was 97 (0.37%). Finally, we determined the discounting parameter r and the weight parameter $\lambda$ of MCD in Equation (17) with the first 5000 time points. We selected $r=0.1$ and $\lambda =0.2$ so that the AUC score at $t=2634$ and $t=2635$ would be the maximum. The AUC score was calculated using Equations (18) and (19).

Figure 14 shows the MDL change statistics $\left\{{\mathsf{\Psi }}_t\right\}$ calculated with SMDL [8] (Figure 14, top), the estimated MCAT $a_{t_i}$ (Figure 14, second), logarithm of the estimated MCAS ${log}_{10}{b}_{t_i}$ (Figure 14, third), and logarithm of the estimated MCI ${log}_{10}{s}_{t_i}$ (Figure 14, fourth). We also estimated the relative volatility with VD [13,25] (Figure 14, fifth) and the change rate of the MDL change statistics $|({\mathsf{\Psi }}_{t_i}-{\mathsf{\Psi }}_{t_{i-1}})/{\mathsf{\Psi }}_{t_{i-1}}|$ (Figure 14, bottom) for comparison in detecting metachanges along both time and state. For VD, the buffer size B and the reservoir size R were both set to 10. In Figure 14 (top), the red points indicate the detected change points.

We summarize what can be seen for metachange statistics in Figure 14 as follows:

• $t=9663$: The trend of MCI increases roughly after $t=5000$, which can be interpreted as a combination of MCAT and MCAS in Figure 14. The relative volatility and the change rate of the MDL change statistics do not show such a significant sign.
• $t=$ 13,230, 13,231, 17,372, 17,832: For time points between $t=$ 10,000 and $t=$ 15,000, the trend of MCI increases. It is also due to the combination of MCAT and MCAS, but is more influenced by MCAS. It might also be a sign of important events at $t=$ 17,372 and 17,832 as well as $t=$ 13,230 and $t=$ 13,231. The relative volatility increases after $t=$ 13,231, which might be a sign of the important event at $t=$ 17,372. However, the change rate of the MDL change statistics does not show such a significant sign.
• $t=$ 25,440: For time points between $t=$ 20,000 and $t=$ 25,000, the trend of MCI increases with large fluctuations. It is also more influenced by MCAS. It might also be a sign of important events at $t=$ 25,440. The relative volatility increases for the time points, but the change rate of the MDL change statistics does not show such a significant sign.

In summary, we can observe a sign of metachange for each important event. We therefore infer that there might have been some symptoms that should be analyzed using field knowledge.

5. Conclusions

We propose the concept of metachanges along time and state in data streams, and we introduce metachange statistics to quantify metachanges from a unified view with MDL. The key idea of our proposed method is to encode the time intervals and change of states with code lengths in the same fashion. Next, we introduce the novel methodology of MCD. Using synthetic datasets, we empirically demonstrated that the proposed algorithm was highly effective in detecting metachanges along time and state. Using a real dataset, we demonstrated that the proposed algorithm could detect metachanges in both time and state, some of which were overlooked by VD [13] and the MDL change statistics [8]. The estimated metachange statistics might have been a sign of important events.

Future work will be directed toward the theoretical guarantee of metachange statistics, especially integrated metachange statistics. We will also consider how to adapt to a non-stationary data stream by updating the weight parameter $\lambda$ in Equation (17). Other research directions might lie in the extension of metachange statistics to transient periods between change points. Furthermore, metachange detection of model structure change and its change sign is another interesting line of research.