\section{Introduction}

Since the seminal paper by Barab\'asi~\cite{Barabasi2005Origin}, the origin of bursty temporal patterns in human dynamics has been extensively debated for the last few decades~(see the review book~\cite{Karsai2018Bursty} and references therein). Here the bursts indicate the rapidly occurring events in short-time periods that are separated by long periods of low activity. In his paper, Barab\'asi argued that non-Poissonian or bursty patterns observed in the email communication dataset can be understood by considering how individuals manage the incoming emails according to their priorities. In contrast, Malmgren et al.~\cite{Malmgren2008Poissonian} suggested that bursty patterns in the same email dataset can be largely explained by cyclic behaviors such as daily and weekly cycles of humans. Later Jo et al.~\cite{Jo2012Circadian} showed that deseasoning daily and weekly cycles from the time series of mobile phone communication cannot entirely remove the burstiness in the time series, implying that the remaining burstiness might be due to various other factors affecting the human dynamics.

Note that the above studies have focused on the individual time series. However, considering the fact that humans are social beings, the effects of social interaction on the bursty human dynamics should not be ignored. Here a fundamental question arises: Are bursts in human dynamics more likely to be the consequences of intrinsic bursty characteristics of individuals or more driven by the interaction between individuals? Borrowing the terms in the network science~\cite{Borgatti2009Network, Barabasi2016Network, Menczer2020First} we refer to this question as \emph{node burstiness versus link burstiness}. Karsai et al.~\cite{Karsai2012Correlated} addressed a similar question by analyzing the time series of individuals (nodes) and their ties (links) in the mobile phone communication dataset. They found that the bursty behavior of nodes is dominated by that of links incident to the nodes, hence concluded that ``burstiness is a property of the links rather than of the nodes''. Here we remark that the mobile phone communication dataset consists of interaction events only. In particular, each mobile phone call can be described by a tuple $(i,j,t)$ in that a caller $i$ makes a call to a receiver $j$ at time $t$~\cite{Karsai2011Small, Miritello2011Dynamical, Karsai2012Correlated, Saramaki2015Seconds}. Therefore the mobile phone communication dataset has also been studied in the framework of temporal networks~\cite{Holme2012Temporal, Masuda2016Guide, Holme2019Temporal}, because the temporal network can be defined as a set of interaction events or tuples $(i,j,t)$. We note that a temporal network model has recently been suggested to reconcile node burstiness and link burstiness~\cite{Hiraoka2020Modeling} rather than to contrast one with the other. However, such a temporal-network approach may systematically preclude individual activities that are not described by tuples $(i,j,t)$, such as watching a movie alone or posting to a personal blog spontaneously. These events that do not necessarily imply interaction between individuals are referred to as standalone events and denoted by $(i,t)$; they can also be important in understanding the origin of bursty human dynamics.

As for an illustrative case study regarding the origin of bursty human dynamics, one can analyze the edit history of the self-organized online collaborative encyclopedia, Wikipedia~\cite{Yun2019Early, Yasseri2012Dynamics}, among others such as the online forum dataset~\cite{Panzarasa2015Emergence} and Twitter dataset~\cite{Ross2015Understanding}. It is well-known that the temporal patterns of edit sequences of Wikipedia articles as well as those of editors are bursty~\cite{Radicchi2009Human, Yasseri2012Dynamics, Yun2016Intellectual, Wang2016Modeling, Kwon2016Double, Zha2016Unfolding, Gandica2017Stationarity, Jeong2017Interval, Jo2020Bursttree}. Unlike mobile phone communication, in which every event is a call or a message between users, editing Wikipedia articles is not evidently a communicational act between editors and thus could be described as a collection of standalone events. Nevertheless, it has been repeatedly shown that interaction between editors does exist and in fact plays an important role in revising articles. Such interaction can be inferred, e.g., by comparing different versions (or revisions) of the article that were modified by different editors~\cite{Maniu2011Building, Yasseri2012Dynamics, Jurgens2012Temporal, Lerner2017Third, Lerner2020Free}, by considering the orders of editors who edited the same article~\cite{Iba2010Analyzing, Keegan2012Staying, Wu2015Integration, Ashford2019Understanding}, and by relating all the active editors who edited the same article to each other~\cite{Laniado2011Coauthorship}. Other works considered the interaction between editors not only in articles but also in their talk pages~\cite{Turek2010Learning, Jankowski-Lorek2016Verifying}. Using the definition of the pairwise interaction between editors, social networks of editors were constructed and analyzed~\cite{Iba2010Analyzing, Laniado2011Coauthorship, Maniu2011Building, Ashford2019Understanding}. Finally, bipartite networks between articles and editors have also been studied~\cite{Jesus2009Bipartite, Wu2011Characterizing, Keegan2012Editors, Jurgens2012Temporal}.

In this paper, by analyzing the edit sequences of Wikipedia articles we aim to understand the role of individuals in the bursty collective dynamics by looking at how spontaneous individual editors are in initiating bursty periods of editing and to what extent such editors' behaviors are driven by interaction with other editors in those periods. These questions can be formulated in terms of \emph{individual-driven burstiness versus interaction-driven burstiness}, which can be seen as an elaboration of the issue on node burstiness versus link burstiness. To infer the interaction between editors we detect bursty periods for a given timescale~\cite{Karsai2012Universal} from the edit sequence of each Wikipedia article. Here the editors who edited the same article during the same bursty period are assumed to have interacted with each other, irrespective of the contents modified by editors. Our approach is similar to that by Karsai et al.~\cite{Karsai2012Correlated} as bursty periods are detected from a sequence of calls or edits, but different in the sense that we focus on the information on who initiates the bursty periods, i.e., who makes the first edit of the bursty period.

As human activities are sometimes spontaneous and other times interaction-driven, the individual human behavior might not be understood only by one of individual-driven burstiness and interaction-driven burstiness. Therefore, we take an approach of quantifying the degree of initiative (DoI) of the editor of interest in a given article using the statistics of bursty periods containing the editor's edits. Since the bursty periods are detected for a given timescale~\cite{Karsai2012Universal}, the DoI is also a function of the timescale. By scanning the entire range of the timescale, we obtain the DoI curve, from which the area under the curve (AUC) is calculated. A large value of the AUC can be interpreted as the dominance of individual-driven burstiness over interaction-driven burstiness and vice versa. To investigate features in the edit sequences that are relevant in understanding the observed AUC values, we correlate the AUC value with several measures for temporal and editorial correlations. Finally, we confirm the empirical findings by devising and analyzing a model capturing the essential features of edit sequences.

The paper is organized as follows. In Sec.~\ref{sect:empirical} we describe our approach for analyzing the edit sequences of Wikipedia articles by means of bursty periods, the degree of initiative, and the area under the curve. In Sec.~\ref{sect:model} we devise and analyze a model capturing the essential features of the edit sequences to confirm the empirical findings in Sec.~\ref{sect:empirical}. Finally, we conclude our work with some remarks on future works in Sec.~\ref{sect:concl}.

\section{Empirical analysis}\label{sect:empirical}

\subsection{Dataset}

To scrutinize the issue on individual-driven burstiness versus interaction-driven burstiness for the origin of bursts in human dynamics, we analyze a dataset of the English Wikipedia dump on October 2, 2015~\cite{English}. The dataset originally contains 11,994,178 general articles, which we call articles in our paper. Among them, we analyze articles that have not been merged with other articles, the number of which is 4,978,964. The number of editors who contributed those articles is 40,057,921. The preprocessed data for each article contains the temporal and editorial information of edits; each edit is recorded with the editor ID and the timestamp in the resolution of seconds. See the sample of the data in Fig.~\ref{fig1}(a,b). The data analyzed in our work are available in the public repository~\cite{ChoiTime}. We remark that we do not take any other information on edits, such as the size and content of edits, into account for the analysis.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{./FIG1.pdf}
\caption{Diagram for the data analysis procedure: (a) A part of the preprocessed data of Wikipedia edit history. Each line contains the timestamp, article ID, and editor ID of the edit. (b) A part of the edit sequence of a given article (e.g., the article 5042951) containing the timestamp and the editor ID per edit. (c) Visualization of the part of the edit sequence for the article 5042951. Each edit is denoted by a vertical line in the horizontal time axis. Among editors who edited the article, we choose a particular editor or ego (e.g., the editor 13014), while all other editors are called alters. Edits by the ego and alters are denoted by red solid lines and blue dotted lines, respectively. For a given timescale $\Delta t$, bursts are detected and categorized into three groups: initiated (red area), uninitiated-but-involved (orange area), and uninvolved bursts (blue area). (d) The degree of initiative (DoI), $r$ in Eq.~\eqref{eq:doi}, for the editor 13014 in the article 5042951, with $\Delta t$ in seconds. See the text for the definitions of $\Delta t_1$ and $\Delta t_2$. (e) The DoI curve as a function of the normalized timescale $\Delta t'$ in Eq.~\eqref{eq:normalize}, from which the area under the DoI curve is calculated as $\approx0.45$.}
\label{fig1}
\end{figure}

\subsection{Degree of initiative}\label{subsec:doi}

An edit sequence of each Wikipedia article is given as a time-ordered set of edits. For the article with $n$ edits, the $i$th edit for $i=1,\cdots,n$ is associated with the editor $c_i$ and timing $t_i$. Among the editors of the article, a particular editor is chosen to be called \emph{ego}, while all other editors are called \emph{alters}. For simplicity we do not distinguish alters; each edit is made either by the ego or by the alter, i.e., $c_i\in \{\textrm{E},~\textrm{A}\}$, where E and A stand for the ego and the alter, respectively. Note that from the sequence of edit timings, i.e., $\{t_i\}_{i=1,\cdots,n}$, one obtains the sequence of inter-edit times (IETs) by the definition of $\tau_i\equiv t_i-t_{i-1}$.

To define the interaction between editors, we assume that editors who edited the same article within the same chunk of time period have interacted with each other. Precisely, we adopt the notion of bursty trains or bursts~\cite{Karsai2012Universal, Karsai2018Bursty}: For a given timescale $\Delta t$, a burst is defined as a set of consecutive edits such that the IET between any two consecutive edits in the burst is smaller than or equal to $\Delta t$, while edits in different bursts are separated by IETs larger than $\Delta t$. See Fig.~\ref{fig1}(c) for an illustrative example. Then we categorize the detected bursts, whose number is denoted by $m$, into three groups according to the editorial information: (i) A burst of which the first edit is made by the ego is called an \emph{initiated burst}. (ii) A burst that is not initiated by the ego but contains ego's edits is called an \emph{uninitiated-but-involved burst}. (iii) A burst without ego's edits is called an \emph{uninvolved burst}. The numbers of initiated bursts, uninitiated-but-involved bursts, and uninvolved bursts are respectively denoted by $m_{\textrm{init}}$, $m_{\textrm{inv}}$, and $m_{\textrm{other}}$, satisfying $m=m_{\textrm{init}}+m_{\textrm{inv}}+m_{\textrm{other}}$. Using the first two numbers we define the degree of initiative (DoI) for the ego in the given article as follows:
\begin{align}
    r \equiv\frac{m_\mathrm{init}}{m_\mathrm{init}+m_\mathrm{inv}},
    \label{eq:doi}
\end{align}
which essentially quantifies how often the ego initiates bursts while being active in editing.

For a fixed timescale $\Delta t$, a larger value of $r$ implies more initiative behavior of the ego at that timescale. However, it is not obvious which timescale is the most relevant for understanding the initiative behavior. Therefore, we study the DoI for the entire range of the timescale for systemic investigation. If $\Delta t$ is smaller than the minimum IET of the edit sequence, each edit constitutes a burst of size one on its own, implying that $m_{\textrm{inv}}=0$, hence $r=1$. As $\Delta t$ increases, bursts are merged with each other, and the total number of bursts, $m$, decreases. In particular, $m_{\textrm{init}}$ either decreases or remains the same, but never increases, which is the main driving force for the overall decreasing $r$. The increasing behavior of $r$ is found only when $m_{\textrm{inv}}$ decreases. Precisely, $m_{\textrm{inv}}$ decreases by one whenever an uninitiated-but-involved burst is merged either with the preceding initiated burst or with another uninitiated-but-involved burst. However, if the number of ego's edits is much smaller than the total number of edits, as is the case in our work, it is more common to find a merger of an uninvolved burst and its following initiated burst into an uninitiated-but-involved burst, leading to the overall increasing $m_{\textrm{inv}}$, hence the overall decreasing $r$. Finally, if $\Delta t$ is equal to or larger than the maximum IET in the edit sequence, all edits belong to a single burst. This burst is either an initiated burst or an uninitiated-but-involved burst, ending up with $r=1$ or $0$, respectively, depending on who made the first edit to the article. Such dependence on the initial condition is somewhat excessive and makes the further analysis less robust. For the edit sequences where the ego makes the first edit, to ensure $r=0$ for $\Delta t$ equal to or larger than the maximum IET in the edit sequence, we add a dummy edit by an alter before the beginning of the edit sequence; we suppose the zeroth edit made by an alter $c_0=\textrm{A}$ at time $t_0=t_1-\tau_1$ for some IET $\tau_1$ (to be discussed).

To guarantee the statistical significance of the analysis, we first choose articles with more than $10^4$ edits from the Wikipedia dataset, leaving us with $697$ articles. For each article we choose human editors who have made at least $1\%$ of edits to the article~\footnote{The ego can be either a human editor or a bot~\cite{Zheng2019Roles}. In our work we focus on the behavior of human editors.}. In total there are $3,099$ such editors, each of whom is designated as the ego. As a result we are left with $4,634$ article-ego pairs to be analyzed. As an example, we plot in Fig.~\ref{fig1}(d) the DoI curve for the editor 13014, who edited $559$ times the article 5042951 with $n=20,683$. As expected, the DoI curve overall decreases with the timescale $\Delta t$.

In each of the DoI curves obtained from the dataset, we observe two characteristic timescales: At the first timescale, denoted by $\Delta t_1$, the DoI curve starts to deviate from $r=1$ as uninitiated-but-involved bursts appear, i.e., $m_{\rm inv}>0$. At the second timescale, denoted by $\Delta t_2$, the DoI curve reaches zero as all initiated bursts disappear, i.e., $m_{\rm init}=0$. The timescale $\Delta t_1$ is equal to the minimum IET between the ego's edit and its preceding edit by the alter, provided that the alter's edit is not part of an initiated burst detected at the timescale $\Delta t_1$. This may indicate a tendency that the ego makes edits only after $\Delta t_1$ since the edits by the alters. If editing by different editors in the same burst may imply the interaction between those editors, $\Delta t_1$ could be interpreted as the minimal timescale for the ego's interaction with other editors. Similarly, the timescale $\Delta t_2$ can be a good proxy for the maximum IET between the ego's edit and its preceding edit by the alter. That is, the ego tends to make edits within $\Delta t_2$ since the edits by the alters, implying that $\Delta t_2$ could be interpreted as the maximal timescale for the ego's interaction with other editors. We find in Fig.~\ref{fig2}(a) that the population of editors can be well described by some typical values of $\Delta t_1$ and $\Delta t_2$. Here $\Delta t_1$ turns out to be shorter than the order of one day in almost all article-ego pairs, possibly due to the fact that we have chosen article-ego pairs with relatively active editors. We also find $14$ article-ego pairs (for three distinct egos) with $\Delta t_1=\Delta t_2$, because these pairs have the editorial structure of A$\cdots$AE$\cdots$EA$\cdots$A. It implies that $r$ drops from $1$ to $0$ at the timescale $\Delta t_1$ ($=\Delta t_2$), hence these pairs will be excluded for the further analysis in the next subsection.

Finally, we remark that the number of article-ego pairs to which we have added the dummy zeroth edit by the alter is quite small ($14$ among $4,634$ pairs) and that the IET between the zeroth and first edits, $\tau_1$, has been set to be $\Delta t_1$ not to affect the empirical results for $\Delta t_1$ and $\Delta t_2$.

\subsection{Area under the DoI curve}\label{subsec:auc}

To fully characterize the DoI curves, we also look at the overall decreasing behavior of each of those curves in the range of $[\Delta t_1, \Delta t_2]$. From the empirical results of the DoI curves, we find various decreasing patterns such as in convex, concave, and linear manners. To quantify such diverse behaviors by a single value, we calculate the area under the DoI curve (AUC) in a normalized range of the timescale. For $\Delta t\in [\Delta t_1, \Delta t_2]$, we define a normalized timescale as
\begin{align}
    \Delta t' \equiv \frac{\log \Delta t - \log \Delta t_1}{\log \Delta t_2 - \log \Delta t_1}.
    \label{eq:normalize}
\end{align}
Note that $0\leq \Delta t'\leq 1$. The DoI curve as a function of $\Delta t'$ is used to calculate the AUC, e.g., as shown in Fig.~\ref{fig1}(e). The value of AUC represents the overall initiative behavior of the ego for the entire range of relevant timescales. Using the AUC one can also compare different article-ego pairs more conveniently irrespective of their characteristic timescales of $\Delta t_1$ and $\Delta t_2$. In sum, each DoI curve for the article-ego pair can be characterized in terms of three quantities, namely, $\Delta t_1$, $\Delta t_2$, and the AUC.

Empirical results of $\Delta t_1$, $\Delta t_2$, and the AUC obtained for $4,620$ article-ego pairs are summarized in Fig.~\ref{fig2}(b), where we take the average of AUC values, denoted by $\langle {\rm AUC}\rangle$, for the set of article-ego pairs having the same values of $\Delta t_1$ and $\Delta t_2$, and plot $\langle {\rm AUC}\rangle$ in the space of $(\Delta t_1, \Delta t_2)$. The value of $\langle {\rm AUC}\rangle$ turns out to be overall independent of $\Delta t_1$ and $\Delta t_2$.

\begin{figure}[t]
\includegraphics[width=\columnwidth]{./FIG2.pdf}
\caption{(a) Heatmap of the number of article-ego pairs having the same characteristic timescales $\Delta t_1$ and $\Delta t_2$. (b) Heatmap of the averaged value of the area under the DoI curve, denoted by $\langle {\rm AUC}\rangle$, for the pairs having the same characteristic timescales $\Delta t_1$ and $\Delta t_2$. In panels (a,~b), log-binned $\Delta t_1$ and $\Delta t_2$ have been used for the calculation. Both horizontal and vertical dashed lines denote one day. (c--e) Heatmaps of the number of article-ego pairs in the spaces of ($B_{_{\rm E}}$, AUC) (c), ($B_{_{\rm A}}$, AUC) (d), and ($M_c$, AUC) (e). See the text for the definitions of $B_{_{\rm E}}$, $B_{_{\rm A}}$, and $M_c$. In each of panels (c--e), we also present the value of the Pearson correlation coefficient (PCC). (f) Heatmap of $\langle {\rm AUC}\rangle$ for the pairs having the same values of $B_{_{\rm E}}$ and $M_c$. In all panels the gray color means no data in the area.}
\label{fig2}
\end{figure}

We investigate features in the edit sequences that are relevant in understanding the observed AUC values. For this, we characterize the edit sequences by two types of features, i.e., temporal correlations and editorial correlations. As for the features of the first type, we consider the set of IETs followed by the ego's edits, i.e., $\mathcal{T}_{_{\rm E}}\equiv  \{\tau_i|c_i={\rm E}\}$, as well as the set of IETs followed by alters' edits, i.e., $\mathcal{T}_{_{\rm A}}\equiv \{\tau_i|c_i={\rm A}\}$. Each IET in $\mathcal{T}_{_{\rm E}}$ indicates how long the ego waits for the next edit since the latest edit, whether it was made by the ego or by the alter. Thus heterogeneity of IETs in $\mathcal{T}_{_{\rm E}}$ can reveal the temporal property of the ego's editing behavior. In general, heterogeneity of interevent times in event sequences has been extensively discussed to show the temporal correlations in various empirical time series~\cite{Karsai2018Bursty}. To quantify the heterogeneity of IETs in $\mathcal{T}_{_{\rm E}}$, we adopt the burstiness measure suggested for finite event sequences~\cite{Kim2016Measuring}:
\begin{align}
    B_{_{\rm E}} \equiv \frac{\sqrt{n_{_{\rm E}}+1} \sigma_{_{\rm E}} - \sqrt{n_{_{\rm E}}-1}\mu_{_{\rm E}}}{(\sqrt{n_{_{\rm E}}+1}-2)\sigma_{_{\rm E}}+\sqrt{n_{_{\rm E}}-1}\mu_{_{\rm E}}},
    \label{eq:burstiness}
\end{align}
where $n_{_{\rm E}}$ is the number of IETs in $\mathcal{T}_{_{\rm E}}$. $\sigma_{_{\rm E}}$ and $\mu_{_{\rm E}}$ denote the standard deviation and average of IETs in $\mathcal{T}_{_{\rm E}}$, respectively. The range of $B_{_{\rm E}}$ is $[-1,1]$. The positive $B_{_{\rm E}}$ indicates a bursty temporal pattern in the ego's editing behavior, while $B_{_{\rm E}}=0$ if the ego's editing behavior can be described by a Poisson process. The negative $B_{_{\rm E}}$ is observed for relatively regular temporal patterns. We also consider the burstiness measure for the alters obtained from $\mathcal{T}_{_{\rm A}}$, denoted by $B_{_{\rm A}}$, in a similar manner.

The second type of feature we consider is the correlation between two consecutive editors in the editor sequence, which can be quantified in terms of the Pearson correlation coefficient (PCC):
\begin{align}
    \label{eq:memory_c}
    M_c \equiv \frac{1}{n - 1}\sum_{i=1}^{n-1}\frac{ [h(c_i) - \mu_1] [h(c_{i+1}) - \mu_2]}{\sigma_1 \sigma_2},
\end{align}
where the editorial information is transformed to the numerical value by the function $h$, namely, $h(c_i={\rm E})=1$ and $h(c_i={\rm A})=0$. $\mu_1$ ($\mu_2$) and $\sigma_1$ ($\sigma_2$) denote the average and standard deviation of $h$ values except for the last (the first) editor, respectively. Positive $M_c$ implies a tendency of the ego's edit (alter's edit) to be followed by the ego's edit (alter's edit). The opposite tendency is implied by the negative $M_c$, while $M_c=0$ indicates the absence of correlations between two consecutive editors.

For each article-ego pair we calculate the values of $B_{_{\rm E}}$, $B_{_{\rm A}}$, and $M_c$ to correlate them with the AUC value for the pair. Figure~\ref{fig2}(c--e) shows heatmaps of the number of article-ego pairs in the spaces of ($B_{_{\rm E}}$, AUC), ($B_{_{\rm A}}$, AUC), and ($M_c$, AUC), respectively. Then for each case, we calculate the value of PCC $\rho$: We observe (i) a negative correlation between $B_{_{\rm E}}$ and AUC values ($\rho=-0.15$), (ii) a slightly positive correlation between $B_{_{\rm A}}$ and AUC values ($\rho=0.063$), and (iii) a strongly positive correlation between $M_c$ and AUC values ($\rho=0.68$).

The observation (i) implies the tendency that the stronger temporal correlation in the ego's editing behavior (i.e., larger $B_{_{\rm E}}$) leads to the smaller AUC values. This can be intuitively understood by the following argument: The larger burstiness measure indicates a more heterogeneous IET distribution, typically with a higher peak at small IETs and a heavier tail at large IETs. As $\Delta t$ increases, small IETs followed by the ego's edits enhance the merger of initiated bursts and their preceding bursts, resulting in the fast decaying $m_{\rm init}$ and DoI curve. On the other hand, for the range of large $\Delta t$, large IETs followed by the ego's edits slow down the merger of initiated bursts and their preceding bursts, resulting in the slow decaying $m_{\rm init}$ and DoI curve. Combining these two effects, one gets the initially fast-decaying and then slow-decaying DoI curve, hence the overall convex DoI curve. The more convex a curve is, the smaller the AUC value is calculated from it. This argument explains the observed negative correlation between $B_{_{\rm E}}$ and AUC values.

The observation (ii), i.e., the slightly positive correlation between $B_{_{\rm A}}$ and AUC values, can be explained similarly but in an opposite manner to the case with the observation (i). More importantly, we find that the effect of $B_{_{\rm A}}$ on the AUC is not as significant as that of $B_{_{\rm E}}$ on the AUC. It is probably because the majority of edits in the article are made by the alters, implying that most of IETs followed by the alters' edits have only marginal effects on the statistics of $m_{\rm init}$ and $m_{\rm inv}$.

Finally, the observation (iii), i.e., a strongly positive correlation between $M_c$ and AUC values, allows us to argue that as the ego's edits become more segregated from the alters' edits due to stronger editorial correlation (i.e., larger $M_c$), the number of uninitiated-but-involved bursts, $m_{\rm inv}$, would be smaller than otherwise. The overall smaller $m_{\rm inv}$ leads to the more elevated DoI curve [see Eq.~\eqref{eq:doi}], hence the larger AUC value.

Let us briefly discuss implications of the empirical results. For the sake of simplicity, we assume that the temporal properties of the alters' edits, characterized by $B_{_{\rm A}}$, are given or fixed. This assumption allows us to focus on the effects of the ego's temporal pattern and the editorial correlation on the AUC value. On the one hand, a large AUC value would be achieved for a small $B_{_{\rm E}}$ and large $M_c$, corresponding to the top left area in Fig.~\ref{fig2}(f). Let us consider an edit sequence where the ego's edits are well segregated from those of alters (i.e., large $M_c$), while the ego's temporal pattern is largely homogeneous ($0\lesssim B_{_{\rm E}}\ll 1$). It implies that the ego is rarely affected by interaction with other editors and more likely to be described by the individual behavior. Therefore, the ego's editorially correlated but temporarily random behavior can be interpreted as the individual-driven burstiness. See Fig.~\ref{fig3}(a) for an exemplary time series corresponding to the case dominated by the individual-driven burstiness.

\begin{figure*}[t]
\includegraphics[width=0.85\textwidth]{./FIG3.pdf}
\caption{(a) A part of time series for the article-ego pair that is characterized by small $B_{_{\rm E}}$ and large $M_c$ (hence, individual-driven burstiness). (b) A part of time series for the article-ego pair that is characterized by large $B_{_{\rm E}}$ and small $M_c$ (hence, interaction-driven burstiness). Red (blue) vertical lines denote the ego's (the alter's) edits. Time intervals have been chosen to contain $20$--$30$ edits by the ego in both cases.}
\label{fig3}
\end{figure*}

On the other hand, a small AUC value can be found for a large $B_{_{\rm E}}$ and small $M_c$, corresponding to the bottom right area in Fig.~\ref{fig2}(f). One can consider an edit sequence where the ego's edits are well mixed with the alters' edits (small $M_c$), and the IETs between the ego's edits and their preceding edits by alters are mostly small but often very large (large $B_{_{\rm E}}$). Therefore, the ego seems to interact with other editors so that the ego's such behavior can be described by the interaction-driven burstiness. See Fig.~\ref{fig3}(b) for an exemplary time series corresponding to the case dominated by the interaction-driven burstiness.

In conclusion, by correlating the observed AUC value with the burstiness measures for the IETs followed by the ego and by the alters and the Pearson correlation coefficient between two consecutive editors, we find that weaker (stronger) temporal correlations in the ego's (alters') editing behavior and/or stronger editorial correlations tend to result in the larger values of AUC. These empirical findings are explained in terms of individual-driven burstiness and interaction-driven burstiness. For more systematic investigation of the mechanisms behind such tendencies we study a model with tunable temporal and editorial correlations in the next Section.

\section{Model}\label{sect:model}

To investigate the mechanisms behind the empirical observation for the DoI curves and corresponding AUC values, we devise and analyze a model that generates edit sequences with tunable temporal and editorial correlations. For the temporal correlations, we consider the distribution of the IETs followed by the ego's edits, denoted by $P_{\rm _E}(\tau)$. We also consider the IET distribution for the alter, denoted by $P_{\rm _A}(\tau)$, in a similar manner. For the editorial correlations, we consider the conditional probability $p(c_i|c_{i-1})$, namely, $q_{\rm _{E|E}}$, $q_{\rm _{A|E}}$, $q_{\rm _{E|A}}$, and $q_{\rm _{A|A}}$. For example, $q_{\rm _{E|A}}$ denotes the probability that the ego's edit follows the alter's edit, irrespective of the timings of edits. These probabilities are not independent of each other but to satisfy the following conditions:
\begin{eqnarray}
q_{\rm _{E|E}} + q_{\rm _{A|E}} &=& 1,\label{eq:q1}\\
q_{\rm _{E|A}} + q_{\rm _{A|A}} &=& 1,\label{eq:q2}\\
q_{\rm _{A|E}}q_{\rm _E} &=& q_{\rm _{E|A}}q_{\rm _A},
\label{eq:q3}
\end{eqnarray}
where $q_{\rm _E}$ ($q_{\rm _A}=1-q_{\rm _E}$) denotes the fraction of the ego's edits (alter's edits) in the edit sequence. The condition in Eq.~\eqref{eq:q3} indicates the asymptotic balance between the frequency of the sequence ``EA'' and that of ``AE'' in the editor sequence. As a result, we are left with two independent parameters, i.e., $q_{\rm _E}$ and $q_{\rm _{E|E}}$.

To generate an edit sequence, we begin with an alter's edit at the initial time, i.e., $c_0={\rm A}$ and $t_0=0$. Then the $i$th edit for $i=1,\cdots,n$ is generated as follows:
\begin{enumerate}
\item The editor of the $i$th edit, i.e., $c_i\in \{{\rm E}, {\rm A}\}$, is determined using $p(c_{i}|c_{i-1})$, conditioned by the editor of the $(i-1)$th edit, i.e., $c_{i-1}$.
\item If $c_i={\rm E}$ (A), the preceding IET $\tau_i$ is randomly drawn from $P_{\rm _E}(\tau)$ [$P_{\rm _A}(\tau)$] to determine the timing of the $i$th edit as $t_{i}=t_{i-1}+\tau_{i}$.
\end{enumerate}
Once the edit sequence of $n+1$ edits is generated, it is analyzed to obtain the DoI curve and its corresponding AUC value.

We analyze our model by deriving the analytic form of the DoI, $r$, in terms of $q_{\rm _E}$ and $q_{\rm _{E|E}}$ as well as the arbitrary functional forms of the IET distributions $P_{\rm _E}(\tau)$ and $P_{\rm _A}(\tau)$. Let us consider an edit sequence of $n+1$ edits, i.e., $n$ IETs. The distribution of these IETs, denoted by $P_{\rm all}(\tau)$, is related to $P_{\rm _E}(\tau)$ and $P_{\rm _A}(\tau)$ as
\begin{align}
    P_{\rm all}(\tau) = q_{\rm _E} P_{\rm _E}(\tau) + q_{\rm _A} P_{\rm _A}(\tau).
\end{align}
We denote by $m$ the number of bursts detected from the edit sequence for a given timescale $\Delta t$. We define the fractions of bursts of different types discussed in the previous Section as follows:
\begin{eqnarray}
    p_{\rm init} &\equiv& \frac{m_{\rm init}}{m},\\
    p_{\rm inv} &\equiv& \frac{m_{\rm inv}}{m},\\
    p_{\rm other} &\equiv& \frac{m_{\rm other}}{m}.
\end{eqnarray}
Then the DoI in Eq.~\eqref{eq:doi} is rewritten as
\begin{align}
    r=\frac{p_{\rm init}}{1-p_{\rm other}}.
    \label{eq:anal_r}
\end{align}

Firstly, we derive the analytic form of $p_{\rm init}$. Since each IET larger than $\Delta t$ separates two consecutive bursts detected at the timescale $\Delta t$, the number of IETs larger than $\Delta t$, denoted by $n'$, is related to $m$ as
\begin{align}
    n'=m-1.
\end{align}
Further, $n'$ can also be written in terms of the cumulative IET distribution as
\begin{align}
    n'=n[1-F(\Delta t)],
\end{align}
where
\begin{align}
    F(\Delta t)\equiv \int_0^{\Delta t} P_{\rm all}(\tau)d\tau.
\end{align}
Therefore, one gets for $m\gg 1$~\cite{Jo2017Modeling, Jo2019Bursty}
\begin{align}
    m\simeq n[1-F(\Delta t)].
    \label{eq:m}
\end{align}
Each initiated burst is detected whenever an ego's edit follows any previous edit after an IET larger than $\Delta t$. Precisely,
\begin{align}
    m_{\rm init}=n \Pr[c_i={\rm E}\ \cap \tau_i> \Delta t] = n q_{\rm _E}[1-F_{\rm _E}(\Delta t)],
    \label{eq:minit}
\end{align}
where
\begin{align}
    F_{\rm _E}(\Delta t)\equiv \int_0^{\Delta t} P_{\rm _E}(\tau)d\tau.
\end{align}
We drop $\Delta t$ from now on for the sake of simplicity. From Eqs.~\eqref{eq:m} and~\eqref{eq:minit} we obtain $p_{\rm init}$ as
\begin{align}
    p_{\rm init} = \frac{q_{\rm _E}(1-F_{\rm _E})}{1-F}.
    \label{eq:anal_pinit}
\end{align}

Secondly, for the derivation of $p_{\rm other}$ one needs to know the fraction of bursts only consisting of the alter's edits among bursts whose first edit is made by the alter. For this, we define the fraction of bursts only with the alter's edits as
\begin{align}
    f\equiv \frac{m_{\rm other}}{m_{\rm inv}+m_{\rm other}}
    \label{eq:f_define}
\end{align}
to get
\begin{align}
    p_{\rm other}=(1-p_{\rm init})f.
    \label{eq:anal_pother}
\end{align}
Here $f$ can be explicitly written as
\begin{align}
    f=\frac{\sum_{b=1}^{\infty}\Pr[c_2=\cdots=c_b={\rm A} |c_1={\rm A}]} {\sum_{b=1}^{\infty} \sum_{c_2\in\{{\rm E},{\rm A}\}}\cdots \sum_{c_b\in\{{\rm E},{\rm A}\}} \Pr[c_2,\cdots,c_b|c_1={\rm A}]}.
    \label{eq:f_pr}
\end{align}
To calculate the numerator on the right hand side of Eq.~\eqref{eq:f_pr}, we consider the probability that once an alter's edit is made, the next event is also made by the alter within an IET smaller than or equal to $\Delta t$. This probability is given by $q_{\rm _{A|A}} F_{\rm _A}$, where
\begin{align}
    F_{\rm _A}(\Delta t)\equiv \int_0^{\Delta t} P_{\rm _A}(\tau)d\tau.
\end{align}
Provided that a burst is initiated by an alter ($c_1 = {\rm A}$), the probability that this burst becomes an uninvolved burst of size $b$ is given by the product of $(q_{\rm _{A|A}} F_{\rm _A})^{b-1}$ and the probability that an IET is larger than $\Delta t$. The latter probability is to guarantee that the burst size is exactly $b$. Therefore, the numerator on the right hand side of Eq.~\eqref{eq:f_pr} is calculated as
\begin{align}
    \sum_{b=1}^{\infty} (q_{\rm _{A|A}} F_{\rm _A})^{b-1} (1-F)=\frac{1-F}{1-q_{\rm _{A|A}}F_{\rm _A}}.
    \label{eq:anal_fnum}
\end{align}

\begin{figure*}[!th]
\includegraphics[width=0.9\linewidth]{./FIG4.pdf}
\caption{(a--c) Analytic results of the DoI curve of the model in Eq.~\eqref{eq:doi_final} using $P_{_{\rm E}}(\tau)$ and $P_{_{\rm A}}(\tau)$ both in a form of the power-law distribution with exponential cutoff in Eq.~\eqref{eq:powerlaw} (solid curves). In all cases we use $q_{\rm _E}=0.1$, $\tau_{\rm min}=2^2$ and $\tau_{\rm c}=2^{10}$ for $P_{_{\rm E}}(\tau)$, and $\tau_{\rm min}=1$ and $\tau_{\rm c}=2^{15}$ for $P_{_{\rm A}}(\tau)$. We plot the DoI curves for several values of $\alpha_{_{\rm E}}$ when $q_{_{\rm E|E}}=0.1$ and $\alpha_{_{\rm A}}=1.5$ (a), for several values of $\alpha_{_{\rm A}}$ when $q_{_{\rm E|E}}=0.1$ and $\alpha_{_{\rm E}}=1.5$ (b), and for several values of $q_{_{\rm E|E}}$ when $\alpha_{_{\rm E}}=\alpha_{_{\rm A}}=1.5$ (c). These analytic results are confirmed by the simulation results from $100$ generated edit sequences of up to $n=2^{23}$ edits (symbols). The error bars denote the standard deviations. (d--f) Numerical results of the AUC value calculated from the analytic result of the DoI curve in Eq.~\eqref{eq:doi_final} for various combinations of the parameter values, i.e., in the space of $(\alpha_{_{\rm A}},\alpha_{_{\rm E}})$ when $q_{_{\rm E}}=q_{_{\rm E|E}}=0.1$ (d), in the space of $(q_{_{\rm E}},q_{_{\rm E|E}})$ when $\alpha_{_{\rm E}}=\alpha_{_{\rm A}}=1.5$ (e), and in the space of $(\alpha_{_{\rm E}},q_{_{\rm E|E}})$ when $q_{_{\rm E}}=0.1$ and $\alpha_{_{\rm A}}=1.5$ (f). In all cases we use $\tau_{\rm min}=2^2$ and $\tau_{\rm c}=2^{10}$ for $P_{_{\rm E}}(\tau)$ and $\tau_{\rm min}=1$ and $\tau_{\rm c}=2^{15}$ for $P_{_{\rm A}}(\tau)$. The gray area in the panel (e) shows the nonexistent parameter space due to conditions in Eqs.~\eqref{eq:q1}--\eqref{eq:q3}.}
\label{fig4}
\end{figure*}

The denominator on the right hand side of Eq.~\eqref{eq:f_pr} can be calculated by enumerating all possible combinations of E and A following the first edit in the burst made by the alter. Let us rewrite each term in the summation over $b$ as $k_b(1-F)$, i.e.,
\begin{align}
k_b\equiv  \frac{\sum_{c_2\in\{{\rm E},{\rm A}\}}\cdots \sum_{c_b\in\{{\rm E},{\rm A}\}} \Pr[c_2,\cdots,c_b|c_1={\rm A}]}{1-F}.
\end{align}
One gets the following results:
\begin{eqnarray}
    k_1 &=& \mathrm{Pr}[\emptyset|\mathrm{A}]=1,\\
    k_2 &=& \mathrm{Pr}[\mathrm{A}|\mathrm{A}]+\mathrm{Pr}[\mathrm{E}|\mathrm{A}]=q_\mathrm{_{A|A}}F_\mathrm{_A}+q_\mathrm{_{E|A}}F_\mathrm{_E},
\end{eqnarray}
and for $b\geq 3$
\begin{eqnarray}
    k_b &=& q_\mathrm{_{A|A}}F_\mathrm{_A} k_{b-1} +
    q_\mathrm{_{E|A}}q_\mathrm{_{A|E}}F_\mathrm{_E} F_\mathrm{_A}\sum_{l=0}^{b-3} (q_\mathrm{_{E|E}} F_\mathrm{_E})^l k_{b-l-2}
    \nonumber\\
    && + q_\mathrm{_{E|A}} q_\mathrm{_{E|E}}^{b-2}{F_\mathrm{_E}}^{b-1}.
    \label{eq:kb3}
\end{eqnarray}
In Eq.~\eqref{eq:kb3} the first term on the right hand side contains all cases with $c_2={\rm A}$, while the final term accounts for the case with $c_2=\cdots =c_b={\rm E}$. The second term includes all other cases. Summing up $k_b$ over all $b$s, we get the denominator on the right hand side of Eq.~\eqref{eq:f_pr} as
\begin{align}
    &\sum_{b=1}^{\infty}k_b (1-F)\nonumber\\ &=\frac{(1-q_\mathrm{_{E|E}} F_\mathrm{_E}+q_\mathrm{_{E|A}}F_\mathrm{_E})(1-F)}{(1-q_\mathrm{_{A|A}} F_\mathrm{_A})(1-q_\mathrm{_{E|E}} F_\mathrm{_E})-q_\mathrm{_{E|A}}q_\mathrm{_{A|E}}F_\mathrm{_E} F_\mathrm{_A}}.
\label{eq:anal_fden}
\end{align}
Then by using Eqs.~\eqref{eq:anal_r},~\eqref{eq:anal_pinit},~\eqref{eq:anal_pother},~\eqref{eq:anal_fnum}, and~\eqref{eq:anal_fden}, for a given timescale $\Delta t$ we finally obtain the analytic result of the DoI for the entire range of $q_{\rm _E}$ and $q_{\rm _{E|E}}$ as well as for arbitrary functional forms of $P_{\rm _E}(\tau)$ and $P_{\rm _A}(\tau)$:
\begin{widetext}
\begin{align}
    r = \frac{ q_{\rm _E}(1-F_{\rm _E})(1-q_{\rm _{A|A}}F_{\rm _A})(1-q_{\rm _{E|E}}F_{\rm _E}+q_{\rm _{E|A}}F_{\rm _E})}{(1-F) q_{\rm _{E|A}}F_{\rm _E}(1-q_{\rm _{A|A}}F_{\rm _A}+q_{\rm _{A|E}}F_{\rm _A}) +q_{\rm _E}(1-F_{\rm _E}) [(1-q_\mathrm{_{A|A}} F_\mathrm{_A})(1-q_\mathrm{_{E|E}} F_\mathrm{_E})-q_\mathrm{_{E|A}}q_\mathrm{_{A|E}}F_\mathrm{_E} F_\mathrm{_A}]}.
    \label{eq:doi_final}
\end{align}
\end{widetext}

\begin{figure*}[!th]
\includegraphics[width=0.88\linewidth]{./FIG5.pdf}
\caption{(a) Non-monotonic dependence of the burstiness parameter $B$ in Eq.~\eqref{eq:originalB} on the power-law exponent $\alpha$ for the IET distribution given by Eq.~\eqref{eq:powerlaw} with $\tau_{\rm min}=1$ and various values of $\tau_c$. (b) The aggregate distribution of IETs followed by the ego's edits over all article-ego pairs (red circles), with the estimated values of power-law exponent $\alpha_{\rm{_E}}=1.15(1)$ (depicted by a dashed line) and burstiness measure $B_{\rm{_E}}\approx 0.69$. (c) The aggregate distribution of IETs followed by the alters' edits over all article-ego pairs (red circles), with the estimated values of power-law exponent $\alpha_{\rm{_A}}=1.11(2)$ (depicted by a dashed line) and burstiness measure $B_{\rm{_A}}\approx 0.74$.}
\label{fig5}
\end{figure*}

We demonstrate our analytic result in Eq.~\eqref{eq:doi_final} by adopting the power-law IET distributions with the power-law exponent $\alpha$, lower bound $\tau_{\rm min}$, and exponential cutoff $\tau_{\rm c}$:
\begin{align}
    P(\tau) = \frac{\tau_{\rm c}^{\alpha-1}}{\Gamma(1-\alpha,\tau_\mathrm{min}/\tau_{\rm c})} \tau^{-\alpha} e^{-\tau/\tau_{\rm c}}\theta(\tau-\tau_{\rm min}),
    \label{eq:powerlaw}
\end{align}
where $\Gamma(\cdot,\cdot)$ is the upper incomplete Gamma function and $\theta(\cdot)$ is the Heaviside step function. This choice is based on empirical results in the literature~\cite{Karsai2018Bursty}. For denoting the parameters for the IET distribution of the ego (the alter) we add the subscript E (A) to those parameters, such as $\alpha_{\rm _E}$ ($\alpha_{\rm _A}$). The analytic results of the DoI curve for various combinations of parameter values of $\alpha_{\rm _E}$, $\alpha_{\rm _A}$, and $q_{\rm _{E|E}}$ are shown as solid curves in Fig.~\ref{fig4}(a--c), where in all cases we use $q_{\rm _E}=0.1$, $\tau_{\rm min}=2^2$ and $\tau_{\rm c}=2^{10}$ for the ego, and $\tau_{\rm min}=1$ and $\tau_{\rm c}=2^{15}$ for the alter. These analytic results are successfully confirmed by the simulation results with $100$ generated edit sequences of up to $n=2^{23}$ edits, as depicted by symbols in Fig.~\ref{fig4}(a--c).

Next, we numerically calculate the AUC value from the analytic result of the DoI curve in Eq.~\eqref{eq:doi_final}, for which we identify the values of $\Delta t_1$ and $\Delta t_2$. Since $\Delta t_1$ is the maximum value of $\Delta t$ satisfying $m_{\rm inv}=0$, i.e., $f=1$ in Eq.~\eqref{eq:f_define}, $\Delta t_1$ turns out to be the same as $\tau_{\rm min}$ for the ego's IET distribution. $\Delta t_2$ is determined as the minimum value of $\Delta t$ satisfying $m_{\rm init}=0$ or $p_{\rm init}=0$, i.e., $F_{\rm _E}=1$ by Eq.~\eqref{eq:anal_pinit}. However, since $F_{\rm _E}(\Delta t)=1$ can be achieved only when $\Delta t\to \infty$, we instead obtain $\Delta t_2$ satisfying the condition $F_{\rm _E}(\Delta t_2)=1-10^{-6}$, which is essentially of the order of $\tau_c$ for $P_{\rm _E}(\tau)$.

Figure~\ref{fig4}(d--f) shows how the AUC value depends on the parameter values for temporal and editorial correlations. We observe in Fig.~\ref{fig4}(d) that the smaller $\alpha_{\rm _E}$ and/or larger $\alpha_{\rm _A}$ lead to the larger AUC values. To compare this finding with the empirical observation in Subsec.~\ref{subsec:auc}, one needs to understand the relation between the shape of the IET distribution and the burstiness parameter derived from it. Once the functional form of the IET distribution $P(\tau)$ is given as in Eq.~\eqref{eq:powerlaw}, the standard deviation $\sigma$ and the mean $\mu$ of the IET are calculated to obtain the value of the burstiness parameter~\cite{Goh2008Burstiness}:
\begin{align}
    B\equiv\frac{\sigma-\mu}{\sigma+\mu},
    \label{eq:originalB}
\end{align}
which can also be derived from Eq.~\eqref{eq:burstiness} in the case with $n\to\infty$. As a result, we observe the non-monotonic dependence of $B$ on the power-law exponent $\alpha$ for fixed values of $\tau_{\rm min}$ and $\tau_c$, as depicted in Fig.~\ref{fig5}(a). This implies that a heavier tail of the IET distribution with a smaller value of $\alpha$ does not necessarily lead to the larger value of $B$, in particular, in the presence of the exponential cutoff to the power-law tail. For the IET distributions with $\tau_{\rm min}=1$ and various values of $\tau_c$, we find an increasing behavior of $B$ as a function of $\alpha$ for the range of $\alpha\lesssim 1.5$. Assuming that the empirical IET distributions of the ego and the alters follow the functional form in Eq.~\eqref{eq:powerlaw}, we roughly estimate values of the power-law exponent of the empirical IET distribution for the ego and for the alters to find them within the range of $\alpha\lesssim 1.5$, as shown in Fig.~\ref{fig5}(b,~c). Conclusively, the analytic result that the smaller $\alpha_{\rm _E}$ and larger $\alpha_{\rm _A}$ lead to the larger AUC value is consistent with the empirical finding of the negative correlation between $B_{\rm _E}$ and AUC values and of the positive correlation between $B_{\rm _A}$ and AUC values.

In Fig.~\ref{fig4}(e) we observe that the larger $q_{\rm _{E|E}}$ for a fixed $q_{\rm _E}$ leads to the larger AUC values, which indeed confirms the empirical tendency that the larger AUC values are observed for the stronger editorial correlations as shown in Fig.~\ref{fig2}(e). Finally, Fig.~\ref{fig4}(f) shows how the AUC value varies when the effects of $\alpha_{\rm _E}$ and $q_{\rm _{E|E}}$ interplay with each other.

By devising and analyzing the model we could understand the underlying mechanisms behind the empirical results more rigorously in terms of the effects of temporal and editorial correlations on the DoI curves and their corresponding AUC values. Thus our modeling approach helps us to better understand the issue on individual-driven burstiness versus interaction-driven burstiness in human dynamics.

\section{Conclusion}\label{sect:concl}

Although the origin of bursty temporal patterns observed in various human behaviors has been extensively investigated, its understanding still remains incomplete~\cite{Karsai2018Bursty}. Considering the fact that humans are social beings, one can ask the fundamental question of whether the bursty human dynamics is dominated by the characteristics of individuals or by the interaction between them. In this paper we have addressed this question by analyzing the Wikipedia edit history to see how spontaneous individual editors are in initiating bursty periods of editing and to what extent such editors' behaviors are driven by interaction with other editors in those periods. This question is referred to as individual-driven burstiness versus interaction-driven burstiness.

After detecting the bursty periods or bursts~\cite{Karsai2012Universal} from the edit sequence of each Wikipedia article, we quantify the degree of initiative (DoI) of an individual editor or ego using the statistics of the bursts containing the ego's edits. All other editors are called alters. Since the bursts are detected for a given timescale, the DoI is also a function of the timescale. Then scanning the entire range of timescale in the article, we obtain the DoI curve, from which the area under the curve (AUC) is calculated. The large value of AUC for an article-ego pair implies the dominance of individual-driven burstiness over interaction-driven burstiness and vice versa. By correlating the AUC value with several measures for temporal and editorial correlations, we find the tendency of the AUC values to be larger for weaker (stronger) temporal correlations of the ego (the alters) and/or stronger editorial correlations in the edit sequences. We also successfully confirm these empirical findings by devising and analyzing a model capturing the essential features of edit sequences. Our approach enables us to better understand the origin and underlying mechanisms of bursts in human social dynamics. We also remark that our approach can be applied to any other time series that can be characterized by a sequence of events with both temporal and contextual information. Here the context of the event indicates a situation in which the event occurs, or any other attributes of the event~\cite{Jo2012Spatiotemporal, Jo2013Contextual}, such as editors of editing events in our case and discussion topics of posted messages in online forum~\cite{Panzarasa2015Emergence}.

Finally, we discuss possible future works. For the comprehensive description of the temporal and editorial correlations, one can consider higher-order temporal and editorial correlations than those used in our work, e.g., by means of a burst-tree decomposition method for characterizing arbitrary-order temporal correlations in event sequences~\cite{Jo2020Bursttree} and $\epsilon$-machine for detecting patterns in the sequence of symbols (e.g., editors in our case)~\cite{Shalizi2001Computational, Shalizi2002Algorithm}. In our work we have focused only on whether the editor of interest initiated bursts or got involved in bursts initiated by other editors. In addition to such initiative behaviors the organization of edits in each burst may also be important to understand the interaction dynamics within bursts, e.g., as studied in Refs.~\cite{Wu2010Evidence, Karsai2012Correlated}. In addition, considering that an editor can edit multiple articles, the same editor may show different initiative behaviors when editing different articles, which can be studied in the future. Further, beyond the article-editor pairs one can extend our approach to encompass the entire web of articles and editors in terms of bipartite temporal networks~\cite{Jurgens2012Temporal, Keegan2012Editors, Zeng2013Trend, Bravo-Hermsdorff2019Gender}. The bipartite temporal network for Wikipedia edit history can consist of editing events, denoted by a tuple $(i,u,t)$, indicating that an editor $i$ makes an edit to an article $u$ at time $t$. Finally, one can also consider information on the interaction between editors from other types of Wikipedia pages, e.g., talk pages~\cite{Viegas2007Talk, Turek2010Learning, Jankowski-Lorek2016Verifying}, or by means of development tools such as Editor Interaction Analyser~\cite{Editor}. One can gain a much deeper insight into the issue on individual-driven burstiness versus interaction-driven burstiness in the context of collective dynamics, thus one can better understand the origin and underlying mechanisms of bursts in human social dynamics.

\begin{acknowledgments}
The authors thank Woo-Sik Son for providing us with the preprocessed dataset of the English Wikipedia, and Heetae Kim and Jinhyuk Yun for useful discussion. H.-H.J. was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF-2018R1D1A1A09081919).
\end{acknowledgments}


\section{INTRODUCTION} \label{sec_intro}
From critical infrastructures and industrial process control to autonomous driving and various biomedical applications, dynamical control systems are increasingly able to be instrumented with new sensing and actuation capabilities. These cyber-physical systems (CPS) comprise growing webs of interconnected feedback loops and must operate efficiently and resiliently in dynamic and uncertain environments. As these systems become large, devising both model-based \cite{murguia2019model,umsonst2018anomaly} and data-driven \cite{renganathan2020distributionally,khojasteh2020learning,martinez_tuning} methods for detecting anomalies (such as component failures or malicious attacks) are critical for their robust and efficient operation. Such critically important cyber-physical networks have become an attractive target to attackers. These systems are large and complex and are often not monitored well enough, enabling attackers to manipulate the system without being detected and cause damage \cite{zetter2014countdown,cherepanov2017industroyer}.

To simplify the analysis and design, often such complex cyber-networks are modeled as a discrete-time linear time invariant system with Gaussian noises. However, this can lead to a significant miscalculation of probabilities and risk if the underlying processes behave differently, for example due to various nonlinearities or malicious attacks. In the context of attacks, it is possible for an attacker to modify the sensor outputs and effectively generate aggressive and strategic noise profiles to sabotage the operation of the system. With stochastic optimization techniques, particularly using the emerging area of distributionally robust optimization (DRO) approaches \cite{dr_wiesemann}, these limitations can be recognized and addressed. DRO enables modelers to explicitly incorporate inherent ambiguity in probability distributions into optimization problems. DRO approaches can be categorized based on the form of the ambiguity set. There are several different parameterizations, including those based on moments, support, directional derivatives \cite{dr_goh}, and Wasserstein balls \cite{dr_peyman}. In practice, we have access to only a finite amount of historical data. However, it is possible to use finite historical data and guarantee resiliency in such critical cyber-physical networks. Here, we propose to use moment-based DRO methods to improve modeling and reduce false alarm rates in cyber-physical networks.

In the context of attacks, the detector tuning has a direct implication on the effect an attacker can have while still remaining stealthy \cite{navid,umsonst2018anomaly}. A model-based approach to attack detection uses a detector that raises alarms when there is a large enough discrepancy between the actual and predicted measurements, a statistic termed the residual. The detector's sensitivity can be increased by decreasing the threshold of detection, but there is an inherent trade-off between sensitivity and the rate at which false alarms are generated. Keeping false alarms to a manageable level requires adjusting sensitivity and the tuning of the detector threshold is typically informed by the distribution of the residual. Authors in \cite{navid_GMM} used Gaussian Mixture Model to approximate the arbitrary noise distributions and obtained a detector threshold corresponding to a desired false alarm rate.

However, when noise distributions are only known to an ambiguity set, traditional tools and approximations no longer suffice to select the threshold and so we turn to a distributionally robust approach. Interest in a DRO-informed perspective on detector tuning is supported by recent work on using a Wasserstein metric \cite{martinez_tuning}. Our moment-based ambiguity set formulation includes all distributions with fixed moments up to some order. The problem of designing anomaly detector thresholds subject to moment constraints of system uncertainties can be addressed using Generalized Moment Problems described in \cite{lasserre2010moments, lasserre2001global}. Further, the conditions for a truncated (finite) moment sequence to represent a probability measure were studied in \cite{curto2008extremal, helton2012semidefinite}. Authors in \cite{bertsimas_popescu} proposed semidefinite programs to compute a probability bound for a random variable lying in a set with known moments up to some order. These techniques can be utilized to design detector thresholds for residual distributions consistent with finite fixed moments up to some order.

\textit{Contributions:} This paper is a significant extension of our previous work \cite{renganathan2020distributionally} where we used a moment-based ambiguity set formulation with fixed first two moments (Proposition \ref{prop_dr_alfa}) to obtain a detector threshold via generalized Chebyshev inequality. The main contributions of the present paper are:
\begin{itemize}
    \item We propose an approach to construct moment-based ambiguity sets with fixed moments up to $k\textsuperscript{th}$ order for the anomaly detection measure and design an anomaly detector threshold for CPSs that exhibit non-Gaussian uncertainties. The approach can utilize either residual moments obtained through a dynamic model or estimated directly from residual data.
    \item We use a semidefinite program (SDP) defined using higher-order moments of the detection measure to find a sharper probability bound for classifying the residual with an improved detector threshold (Theorem \ref{thm_alfa_higher_moments}).
    \item We illustrate empirically that inclusion of higher-order moments results in tightened threshold (Lemma \ref{lemma_1}), thereby restricting an attacker's potential impact, and also prove that the volume of the attack-reachable set shrinks with a tightened threshold (Corollary  \ref{cor_reach_shrinks}).
\end{itemize}
While anomaly detection is widely studied in CPS literature, our distributionally robust approach using higher-order moments of the detection measure data marks the novel contribution of this paper. The rest of the paper is organized as follows. Section \ref{sec_momemt_problem} formulates the problem using moment-based ambiguity sets. In Section \ref{sec_dr_threshold_design}, the design of anomaly detector threshold is discussed. Section \ref{sec_reach_sets} describes the procedure to find the boundary of the reachable sets obtained using distributionally robust tuned detector. Section \ref{sec_numerical_simulation} presents the numerical results with inferences. Finally, Section \ref{sec_conclusion} concludes and summarizes the future research directions.
\section{Problem Formulation Using Higher-Order Moment-Based Ambiguity Sets} \label{sec_momemt_problem}
In this section, we propose a framework for designing an anomaly detector threshold for cyberphysical systems. The approach can utilize either model-based propagation of residual moments or data-driven estimation of detection measure moments directly from data.
\subsection{Model-based Problem Formulation}
Here, we model an uncertain cyber-physical system as a stochastic discrete-time linear system
\begin{align} \label{eqn_uncertain_cps}
    x_{t+1} &= A x_t + B u_t + w_t, \quad t \in \mathbb{N} \\
    y_t &= C x_t + v_t,
\end{align}
where $x_t \in \mathbb{R}^n, u_t \in \mathbb{R}^m$ are the system state and input respectively at time $t$. The matrices $A$ and $B$ denote the system matrix and control input matrix, respectively. The output $y_t \in \mathbb{R}^p$ aggregates a linear combination of the states with the observation matrix $C \in \mathbb{R}^{p \times n}$. We assume that the pair $(A,C)$ is detectable and $(A,B)$ is stabilizable. The process noise $w_t \in \mathbb{R}^{n}$ and the sensor noise $v_t \in \mathbb{R}^{p}$ are modeled as zero-mean random vectors independent and identically distributed across time with covariance matrix $\Sigma_{w}, \Sigma_{v}$ respectively. Let $\kappa = (k_1,\dots,k_n)^{\top}$ with $k_j \in \mathbb{Z}_{+}$ non-negative integers and $J_k = \{ \kappa | \sum^{n}_{i = 1} k_i \leq k \}$. Further, assume that the feasible first $k$ moment sequence of the distributions of the $w_t, v_t$ namely $M^{\kappa}_{w}, M^{\kappa}_{v}, \kappa \in J_k$ respectively are known. The distributions $P_w$ of $w_t$ and $P_v$ of $v_t$ are unknown (and not necessarily Gaussian and possibly heavy-tailed\footnote{For the purposes of this paper, we consider heavy-tailed distributions as those whose moments above a certain order may be infinite, in which case their tails are heavier than a Gaussian. We assume the moments up to order $k$ are finite.}) and will be assumed to belong to the $k$-moments-based ambiguity sets of distributions $\mathcal{P}^{w}_{k}$ and $\mathcal{P}^{v}_{k}$ respectively defined as follows
\begin{align} \label{eqn_ambig_set_w_high}
    \mathcal{P}^{w}_{k} = \{ P_{w} \mid \mathbb{E}[w^{k_1}_{t_1} w^{k_2}_{t_2} \dots w^{k_n}_{t_n} ] = M^{\kappa}_{w}, \ \kappa \in J_k \}, \\
   \mathcal{P}^{v}_{k} = \{ P_{v} \mid \mathbb{E}[v^{k_1}_{t_1} v^{k_2}_{t_2} \dots v^{k_p}_{t_p} ] = M^{\kappa}_{v}, \ \kappa \in J_k \}  \label{eqn_ambig_set_v_high}.
\end{align}
When the actual measurement $y_t$ is corrupted by an additive attack, $\delta_t \in \mathbb{R}^p$, the true output of the system fed to the controller becomes
\begin{align}
    \Bar{y}_t = y_t + \delta_t = C x_t + v_t + \delta_t.
\end{align}
We utilize a steady-state Kalman filter to construct a state estimate $\hat{x}_{t}$ to minimize the squared norm of the estimation error $e_t = x_t - \hat{x}_t$ where,
\begin{align}
    \hat{x}_{t+1} &= A \hat{x}_t + B u_t + L (\bar{y}_t - C \hat{x}_t),
\end{align}
and the estimation error evolves as
\begin{align}\label{eqn_error_dynamics}
    e_{t+1} &= (A - LC) e_t + w_t - L v_t - L \delta_t.
\end{align}
In the absence of attacks (that is, $\delta_t = 0$) and when the covariance matrices of the noises are fixed and known, the Kalman gain $L = PC^{\top}(CPC^{\top} + \Sigma_{v})^{-1}$ minimizes the steady state covariance matrix
\begin{align}
P &:= \lim_{t \rightarrow \infty} P_t := \mathbf{E}[e_t e^{\top}_t].
\end{align}
Since $(A,C)$ is assumed to be detectable, the existence of $P$ is guaranteed and it can be found through the solution of an algebraic Ricatti equation. We define a residual sequence $r_t$ as the difference between the actual received output $\Bar{y}_t$ and the predicted output $C \hat{x}_t$ as,
\begin{align} \label{eqn_r}
    r_{t} &= \bar{y}_t - C \hat{x}_t = C e_t + v_t + \delta_t,
\end{align}
and in the attack free setting, $r_t$ falls according to a zero mean distribution with covariance
\begin{align} \label{eqn_residual}
    \Sigma_{r} = \mathbf{E}[r_t r^{\top}_t] = C P C^{\top} + \Sigma_{v}.
\end{align}
Since \eqref{eqn_r} is linear, it is possible to obtain the fixed and first $k$ moments of the random variable $r_{t}$ by propagating the corresponding moments of the primitive random variables $w_t, v_t$. Thus, the distribution $P_r$ of $r_t$ (not necessarily Gaussian) belongs to an ambiguity set $\mathcal{P}^{r}_{k}$ given by
\begin{align} \label{eqn_ambig_set_r}
    \mathcal{P}^{r}_{k} = \{ P_{r} \mid \mathbb{E}[r^{k_1}_{t_1} r^{k_2}_{t_2} \dots r^{k_p}_{t_p} ] = M^{\kappa}_{r}, \ \kappa \in J_k \}.
\end{align}
We define detection measure $q_t$ as a quadratic function of $r_t$
\begin{align} \label{eqn_q}
    q_t = r^{\top}_t \Sigma^{-1}_{r} r_t,
\end{align}
which will be compared to a threshold for anomaly detection.

\subsection{Data-Driven Moment Estimation from Residual Data}
An alternative to obtaining residual moments by propagating the moments of primitive random variables through the system model is to instead collect residual data $r_t$ (from attack-free operation) and estimate residual moments or the moments of $q_t$ directly from the data. Such a data-driven approach allows our proposed tuning approaches to be used in much broader settings where it is difficult to propagate moments through a model (or even to obtain a model), but where residual data is easily generated from sensors and a state estimator. Higher-order moments require increasingly more data to obtain accurate estimates. Determining the required amount of data is possible using finite-sample measure concentration results given as in  \cite{shawe2003estimating}, but we leave such an analysis for future work. The proposed methods for setting thresholds for anomaly detection can be used together with ambiguity sets built upon data-driven residual estimates, although the false alarm rates will also be affected by sampling errors. Moment estimation uncertainty could be accommodated using the same generalized moment problem computations we propose here, but with assumed \emph{bounds} on moment estimates rather than having them fixed to exact known values. Such problems can also be reformulated in a computationally tractable manner, as described in, e.g., Section 3 of \cite{delage2010distributionally}, which incorporate uncertainty sets for estimated moments. These formulations would simply impose further constraints on our primal problem defined in subsection III.B of our revised manuscript. It is possible and would be interesting to explore how to use statistical confidence intervals for data-driven moment estimates to inform bounds used in the SDP and obtain end-to-end guarantees on false alarm rates associated with certain thresholds. This will be pursued in future work.

In either setting, the \textit{feasible} first $k$ moment sequence of the distribution of $q_t$ denoted by $M^{k}_{q}$ is assumed to be known either from its primitive variables through the model-based approach or estimated from data. Then, the $k$-moments-based ambiguity set of the scalar random variable $q_t$ is defined as
\begin{align} \label{eqn_ambig_set_q_high}
    \mathcal{P}^{q}_{k} = \{ P_{q} \mid  \mathbb{E}[q^{k}_{t}] = M^{k}_{q} \}.
\end{align}
\begin{proposition} \label{def_rs_reach}
\cite{bertsimas_popescu} Consider a univariate random variable $X$ defined on $\Omega=\mathbb{R}_{+}$, endowed with its Borel sigma algebra of events. A sequence $\bar{\sigma} = (M_1, M_2, \dots, M_k)^{\top}$ is a \textbf{feasible} $(1,k,\Omega)-$moment vector of the random variable $X$, if and only if the matrices $R_k \succeq 0$ and $R_{k-1} \succeq 0$, where for any integer $l \geq 0$, the matrices are defined as
\begin{align} \label{eqn_mom_feas_cdtn}
\text {\scriptsize
    $R_{2l} = \begin{bmatrix} 1 & M_1 & \dots & M_l \\ M_1 & M_2 & \dots & M_{l+1} \\ \vdots & \vdots & \ddots & \vdots \\ M_l & M_{l+1} & \dots & M_{2l} \end{bmatrix}, R_{2l+1} = \begin{bmatrix} M_1 & \dots & M_{l+1} \\ M_2 & \dots & M_{l+2} \\ \vdots & \ddots & \vdots \\ M_{l+1} & \dots & M_{2l+1}
    \end{bmatrix}.$}
\end{align}
\end{proposition}

\section{Design of Anomaly Detector Thresholds} \label{sec_dr_threshold_design}
Given a threshold \footnote{The first and second subscripts in the threshold separated by a comma denote the random variable and number of moments respectively.}$\alpha_{q,k} > 0$ and the distance measure $q_t$, alarm time(s) $t^{\star}$ are produced according to the following rules,

\begin{align} \label{eqn_detectorthresh}
    \begin{cases} q_t \leq \alpha_{q,k}, &\text{no alarm} \\
    q_t > \alpha_{q,k}, &\text{alarm: } t^{\star} = t.
    \end{cases}
\end{align}
Even in the absence of attacks, the detector is expected to generate false alarms due to the infinite support of $v_t$, because some values drawn from $P_{q}$ will exceed the threshold $\alpha_{q,k}$. If $P_{q}$ is known, then it is possible to extract an optimum threshold value $\alpha_{q}^*$ from the corresponding cumulative distribution function $F_{q}$ for a desired false alarm rate, $\mathcal{A}$. For example, if $r_t$ is Gaussian, $q_t$ would be a chi-squared random variable, and the optimum threshold $\alpha^{\star}_{q}$ corresponding to the desired false alarm rate $\mathcal{A} = \mathcal{A}^{\star}$ is then
\begin{equation} \label{eqn_chi_squared_threshold}
\alpha^{\star}_{q} = \alpha_{\chi^{2}} := 2 \mathbb{P}^{-1}\left(1-\mathcal{A}^{\star}, \frac{p}{2}\right),
\end{equation}
where $\mathbb{P}^{-1}(\cdot, \cdot)$ denotes the inverse regularized lower incomplete gamma function \cite{murguia2019model}. When the complete distribution is not available, tuning methods using \eqref{eqn_chi_squared_threshold} may design thresholds that generate actual false alarm rates significantly higher than what is desired. With the distributionally robust approach, we aim to achieve a false alarm rate less than $\mathcal{A}$, and the detector threshold $\alpha^{\star}_{r,k}$ is selected such that
\begin{align}  \label{eqn_r_ineq}
\sup_{P_{r} \in \mathcal{P}^{r}_{k}} P_{r}(r^{\top}_{t} \Sigma^{-1}_{r} r_{t} \leq \alpha^{\star}_{r,k}) &= 1 - \mathcal{A}.
\end{align}
When the first two moments $(k = 2)$ of $r_t$ are known, the following proposition using the generalized Chebyshev inequality explained in \cite{boyd_sdp,chen_chebyshev} can be used to obtain the worst case detector threshold $\alpha^{\star}_{r,2}$ satisfying \eqref{eqn_r_ineq}.
\begin{proposition} \label{prop_dr_alfa}
\cite{renganathan2020distributionally} Given a desired false alarm rate $\mathcal{A}$ and $r_t \sim P_{r} \in \mathcal{P}^{r}_{k}$ with $k = 2$, the optimal distributionally robust threshold $\alpha^{\star}_{r,2}$ satisfying \eqref{eqn_r_ineq} is
\begin{align} \label{eqn_DR_alpha}
    \alpha^{\star}_{r,2} = \frac{p}{\mathcal{A}}.
\end{align}
\end{proposition}
\subsection{Improved Detector Threshold With Higher-Order Moments} \label{subsec_k_moments}
It is possible to obtain a sharpened detector threshold than $\alpha^{\star}_{r,2}$, if higher-order moments are taken into account. For example, the skewness and kurtosis parameters convey asymmetry and heaviness of tails of the distribution, respectively. We can leverage such information about the true but unknown distribution revealed by the higher-order moments to tighten the required probability bound and thereby obtain an improved detector threshold. It is possible to use $r_{t}$ with higher order moments (first $k$ moments) to obtain a sharpened detector threshold $\alpha^{\star}_{r,k}$. However, it was shown in \cite{bertsimas_popescu} that for $r_t \in \mathbb{R}^{p}$ with support $\Omega = \mathbb{R}^{p}$ and $k \geq 4$, it is NP-hard to find tight bounds for the corresponding moment bound problem with rational problem data. On the other hand, they provide a semidefinite optimization problem in $k+1$ dimension for the case of a univariate random variable with $k$ moments. Hence, rather looking for higher-order moments of $r_{t}$, instead we look for higher-order moments of scalar random variable $q_{t}$. Subsequently, we use a bisection algorithm to obtain a sharpened detection threshold $\alpha^{\star}_{q,k}$ for a given $\mathcal{A}$.
\subsection{Estimating Probability Using $(1,k,\Omega)-$Moment Bound}
\begin{figure}[ht]
    \centering
    \includegraphics[scale=0.35]{poly.PNG}
    \caption{A moment based polynomial bounding the indicator function $\mathbf{1}_{\mathcal{S}_{k}}$ representing the set $\mathcal{S}_{k} = \mathbb{R}_{> \alpha_{q,k}}$ is shown here.}
    \label{fig_poly}
\end{figure}
Given the first $k$ moments of random variable $q_t$ with support $\Omega = \mathbb{R}_{\geq 0}$ and the set $\mathcal{S}_{k} = \mathbb{R}_{> \alpha_{q,k}}$ representing an alarm event as shown in Figure \ref{fig_poly}, the infinite dimensional $(1,k,\Omega)-$moment bound primal problem is given by
\begin{equation} \label{eqn_1komega_primal}
    \begin{aligned}
        &\sup_{P_{q} \in \mathcal{P}^{q}_{k}} & & P_{q}(q_{t} \in \mathcal{S}_{k})\\
        &\text{s.t.} & &\int_{\Omega} q^{r} d\mu = M^{r}_{q}, \quad r = 0,1,\dots,k.
    \end{aligned}
\end{equation}
Since \eqref{eqn_1komega_primal} is infinite dimensional, it is difficult to solve efficiently. However, we can use the linear programming duality theory to associate a dual variable $y_{r}, r = 0,1,\dots,k$ with each equality constraint of \eqref{eqn_1komega_primal} to get the corresponding dual problem
\begin{equation} \label{eqn_1komega_dual}
    \begin{aligned}
        &\text{min} & & \sum^{k}_{r = 0}y_r M^{r}_q\\
        &\text{s.t.} & &p(q) = \sum^{k}_{r = 0}y_r q^r \geq 1, \quad \forall q \in \mathcal{S}_{k},\\
        & & &p(q) = \sum^{k}_{r = 0}y_r q^r \geq 0, \quad \forall q \in \Omega.
    \end{aligned}
\end{equation}
The probability obtained as the solution to \eqref{eqn_1komega_dual} is an upper bound to the probability associated with \eqref{eqn_1komega_primal}. In general, if the moment vector $\bar{\sigma}_{k} = (M^{0}_q, M^{1}_q, \dots M^{k}_q)$ is an interior point of the set $\mathcal{M}_{k}$ of all feasible moment vectors, then strong duality exists between \eqref{eqn_1komega_primal} and \eqref{eqn_1komega_dual} enabling us to obtain a \emph{tight} bound on $P_{q}(q_{t} \in \mathcal{S}_{k})$. To achieve a desired false alarm rate $\mathcal{A}$, we can tune the threshold $\alpha_{q,k}$ defining the set $\mathcal{S}_{k}$ such that solution to \eqref{eqn_1komega_dual} is $\mathcal{A}$. The following lemma establishes the general trend observed between the values of the tuned threshold $\alpha_{q,k}$ for increasing values of $k$.

\begin{lemma} \label{lemma_1}
Let $\mathcal{A}$ be the desired false alarm rate and $q_{t} \sim P_{q} \in \mathcal{P}^{q}_{k}$ as in \eqref{eqn_ambig_set_q_high}. Then, $\forall j,k \in \mathbb{Z}_{\geq 1}$ with $j < k$, we have
\begin{align}
    \alpha^{\star}_{q,k} \leq \alpha^{\star}_{q,j}
\end{align}
when the associated moments agree up to order $j$.
\end{lemma}
\begin{proof}
When the associated moments up to order $j<k$ agree, clearly we have $\mathcal{P}^{q}_{k} \subseteq \mathcal{P}^{q}_{j}$, since additional moment constraints restrict the set of distributions.
For a fixed threshold $\alpha^{\star}_{q,j}$ tuned for the false alarm rate $\mathcal{A}$, it follows that
\begin{align}
    \sup_{P_{q} \in \mathcal{P}^{q}_{k}} P_{q}(q_{t} > \alpha^{\star}_{q,j}) \leq \sup_{P_{q} \in \mathcal{P}^{q}_{j}} P_{q}(q_{t} > \alpha^{\star}_{q,j}) = \mathcal{A}.
\end{align}
Thus, to achieve a desired false alarm rate $\mathcal{A}$ under the constraint that $P_{q} \in \mathcal{P}^{q}_{k}$, the threshold $\alpha^{\star}_{q,k}$ must be non-increasing in $k$, satisfying $\alpha^{\star}_{q,k} \leq \alpha^{\star}_{q,j}$ for $j < k$.
\end{proof}

\begin{theorem} \label{thm_alfa_higher_moments}
Let $\epsilon > 0$, $k \in \mathbb{Z}_{\geq 2}$ and consider $q_t$ with its first $k$ moments of distribution $(M^1_q,M^2_q,\dots,M^k_q)$ (we let $M^0_q = 1$) defined on $\mathbb{R}_{+}$ being known (or estimated from the given data) and the associated anomaly detector threshold $\alpha_{q,k}$ introduced in \eqref{eqn_detectorthresh}, which is intended to achieve a desired false alarm rate $\mathcal{A}$. Suppose that $\alpha_{q,k}$ is obtained by solving the following bisection algorithm (knowing that the desired $\alpha^{\star}_{q,k} \in [\alpha_{l}, \alpha_{u}]$ and $\alpha_{u} = \alpha^{\star}_{q,k-1}$ by Lemma \ref{lemma_1})
\begin{algorithm}
\caption{\texttt{Bisection} Subroutine}
\While{$\alpha_u - \alpha_l > \epsilon$}
{
$\alpha^{\star}_{q,k} \gets (\alpha_u + \alpha_l)/2$ \\
$p_{min} \gets $\text{ optimal value of} \eqref{eqn_sdp_higher_moments} with $\alpha_{q,k} = \alpha^{\star}_{q,k}$ \\
\eIf{$p_{min} > \mathcal{A}$}
{$\alpha_l \gets \alpha^{\star}_{q,k}$}
{$\alpha_u \gets \alpha^{\star}_{q,k}$}
}
\label{alg_bisect}
\end{algorithm}
where the SDP that gives the tight upper bound on $P_{q}(q_t \geq \alpha_{q,k})$ in the third line is
\begin{equation} \label{eqn_sdp_higher_moments}
    \begin{aligned}
        &\text{min} & & \sum^{k}_{r = 0}y_r M^{r}_{q}\\
        &\text{s.t.} & & \sum_{i,j:i+j=2l-1} x_{ij} = 0, &l = 1,\dots,k, \\
        & & & (y_{0}-1) + \sum^{k}_{r = 1}y_r \alpha^{r}_{q,k} = x_{00}, \\
        & & & \sum^{k}_{r = l}y_r \binom{r}{l} \alpha^{r-l}_{q,k} = \sum_{i,j:i+j=2l} x_{ij}, &l = 1,\dots,k, \\
        & & & X \succeq 0, \\
        & & & \sum_{i,j:i+j=2l-1} z_{ij} = 0, &l = 1,\dots,k, \\
        & & & \sum^{l}_{r = 0}y_r \binom{k-r}{l-r} \alpha^{r-l}_{q,k} = \sum_{i,j:i+j=2l} z_{ij}, &l = 0,\dots,k, \\
        & & & Z \succeq 0.
    \end{aligned}
\end{equation}
with variables $X, Z\in \mathbb{R}^{(k+1) \times (k+1)}, y_r, r = 0,1,\dots,k$. Then with this optimal detector threshold $\alpha^{\star}_{q,k}$, the false alarm rate under the worst-case distribution of the $q_t$ specified by the $k$-moments-based ambiguity set \eqref{eqn_ambig_set_q_high} is at most $\mathcal{A}$.
\end{theorem}
\begin{proof}
To get an upper bound on required probability in \eqref{eqn_1komega_primal}, it suffices to check for the polynomials defined in \eqref{eqn_1komega_dual} to be non-negative in their respective sets. That is, the polynomial $p(q)$ satisfies $p(q) - 1 \geq 0, \forall q \in \mathcal{S}_{k}$ if and only if there exists a positive semidefinite matrix $X \in \mathbb{R}^{(k+1) \times (k+1)}$ that satisfies the first four constraints defined in \eqref{eqn_sdp_higher_moments}. Similarly, the same polynomial $p(q)$ satisfies $p(q) \geq 0, \forall q \in [0, \alpha_{q,k}]$ if and only if there exists a positive semidefinite matrix $Z \in \mathbb{R}^{(k+1) \times (k+1)}$ that satisfies the last three constraints defined in \eqref{eqn_sdp_higher_moments}. Then using the results given in Theorem 3.2 of \cite{bertsimas_popescu}, for a fixed threshold $\alpha_{q,k}$, the optimal value of the semidefinite program in \eqref{eqn_sdp_higher_moments} yields the worst-case false alarm rate generated by the worst-case $P_{q} \in \mathcal{P}^{q}_{k}$ in \eqref{eqn_ambig_set_q_high}. To ensure that the optimal value of \eqref{eqn_sdp_higher_moments} is equal to the desired false alarm rate $\mathcal{A}$, it suffices to tune the threshold $\alpha_{q,k}$ using the bisection algorithm \ref{alg_bisect}, to get $\alpha^{\star}_{q,k}$.
\end{proof}

\subsection{Discussion}
For $k = 1,2,3$, the solution to  \eqref{eqn_sdp_higher_moments} can be obtained in closed-form using Theorem 3.3 of \cite{bertsimas_popescu} with the corresponding moments data without explicitly solving the SDP in \eqref{eqn_sdp_higher_moments}. For $k = 1,2$, we recover the Markov bound and a strictly improved Chebyshev bound respectively as the solution of \eqref{eqn_sdp_higher_moments}. Assuming that $\hat{\delta}_{k} > 0, \alpha^{\star}_{q,k} = (1+\hat{\delta}_{k}) M^{1}_{q}$, the resulting closed form solutions for the optimal detector threshold $\alpha^{\star}_{q,k}$ for a desired false alarm rate $\mathcal{A}$ with $k = 1,2$ are summarized in Table \ref{tab_alfas} with $C^{2}_{M} = \frac{M^{2}_{q} - (M^{1}_{q})^{2}}{(M^{1}_{q})^{2}}$ and these values can be used as $\alpha_{u}$ towards computing $\alpha^{\star}_{q,k}$ for any $k >2$. We omit the expression for the threshold $\alpha^{\star}_{q,3}$ for the sake of brevity.
\begin{table}[ht]
    \centering
    \begin{tabular}{|| c | c ||}
 \hline
 $k$ & $\alpha^{\star}_{q,k}$ value\\ [0.5ex]
 \hline\hline
 1 & $\frac{M^{1}_{q}}{\mathcal{A}}$ \\ \hline
 2 & $\left(1 + \sqrt{\frac{1 - \mathcal{A}}{\mathcal{A}}} \, C_{M} \right)M^{1}_{q}$ \\ [1ex]
 \hline
 \end{tabular}
    \caption{Closed-form solutions for $\alpha^{\star}_{q,k}$ with $k = 1,2.$}
    \label{tab_alfas}
\end{table}
Further, the optimal threshold $\alpha^{\star}_{q,k}$ is tight for a given $\mathcal{P}^{q}_{k}$ up to a tolerance $\epsilon > 0$ specified by the bisection in Algorithm \ref{alg_bisect}.

The exact detector threshold (which we call $\alpha^{\star}_{q}$) corresponding to $P_{q}$ is obtained only asymptotically (or equivalently when true $P_{q}$ is known exactly). That is,
    \begin{align}
        \alpha^{\star}_{q} = \begin{cases} \lim_{k \rightarrow \infty} \alpha^{\star}_{q,k}, &\text{if all moments are known}, \\
        F^{-1}_{q}(\mathcal{A}), &\text{if true $P_{q}, F_{q}$ are known}.
        \end{cases}
    \end{align}
However, while determining how many moments are needed to get a close approximation of the exact threshold is difficult in general, we find that significant improvements can be obtained from a small number of moments. Given $k$ moments of $q_t$, the complexity of this higher-order moment based approach involves solving the SDP given by \eqref{eqn_sdp_higher_moments} with variables $X, Z\in \mathbb{R}^{(k+1) \times (k+1)}, y_r, r = 0,1,\dots,k$. For large $k$, it is advisable to use, e.g., the Legendre polynomial basis instead of the standard polynomial basis for obtaining the moment-based polynomial in \eqref{eqn_sdp_higher_moments} as the former has nice orthogonal properties that improve numerical stability.

\section{Attack-Reachable Set Bounds}\label{sec_reach_sets}
With the attacker assumed to have perfect knowledge of the system dynamics, the Kalman filter, control inputs,  measurements along with read and write access to all the sensors at each time step, we show that the volume of the attack-reachable set shrinks with the tightened threshold. We define a zero-alarm attack, which generates attack sequences so that no alarms are raised during attack. With $\Sigma^{\frac{1}{2}}_{r}$ being the symmetric square root of $\Sigma_{r}$ and an attack input $\bar{\delta}_t$ such that $\bar{\delta}^{\top}_t \bar{\delta}_t \leq \alpha$ (here $\alpha = \alpha^{\star}_{q,k}$), the attack sequence $\delta_t = - Ce_t - v_t + \Sigma^{\frac{1}{2}}_{r} \bar{\delta}_t$ guarantees that no alarm is raised. Then, using a static estimator feedback $u_t = K \hat{x}_t$ with $\bar{\delta}_t$, the evolution of the system dynamics with the joint state $\xi_t=[x_t, e_t]^{\top}$ with input $\zeta_t=[w_t, \bar{\delta}_t]^{\top}$ is studied. We define a reachable set of interest, driven by the ellipsoidally bounded inputs $w_t$ and $\bar{\delta}_t$, as
\begin{equation}
\mathcal{R}_x = \left\{  x_t=[I_n,\, 0_{n\times n}]\xi_t \ \left|\
    \begin{aligned}
        &\xi_{t+1}=\hat{A}\xi_{t}+ \hat{B} \zeta_t, \\
        &\xi_1=\mathbf{0},\ \bar{\delta}^{\top}_t \bar{\delta}_t \leq \alpha,\\
        &w_t^\top\Sigma_w^{-1}w_t \leq \bar{w},\ \forall t\in\mathbb{N}
    \end{aligned}
    \right.
    \right\},
\end{equation}
where the noise threshold $\bar{w}$ obtained using \eqref{eqn_DR_alpha} satisfies,
\begin{equation}\label{eqn_DRnoise_threshold}
    \sup_{P_{w} \in \mathcal{P}^{w}_{k}} P_{w}(w^{\top}_{t} \Sigma^{-1}_{w} w_{t} \leq \bar{w}) = 1 - \mathcal{A},
\end{equation}
and $\hat{A} = \begin{bmatrix}A+BK & -BK \\ 0 & A\end{bmatrix}, \hat{B} = \begin{bmatrix}I & 0\\I & -L\Sigma_r^{1/2}\end{bmatrix}$.
Using the geometric approach presented in \cite{hashemi2020codesign} with $H_i = A_{cl}^{i} - A^{i}$, the reachable set of states is the Minkowski sum of the following ellipsoidal bound
\begin{align}\label{eqn_reach_minkowski}
\text {\footnotesize
$\mathcal{R}_{x,t}(\bar{w}, \alpha)  =\bigoplus_{i=0}^{t-2} \mathcal{E}\left(\bar{w} A^i \Sigma_w {A^{i}}^{\top}\right) \oplus \mathcal{E}\left(\alpha H_i L\Sigma_r L^{\top} H_{i}^{\top}\right).$}
\end{align}
Specifically, theorem 1 of \cite{hashemi2020codesign} provides us the exact boundary of Minkowski sum in \eqref{eqn_reach_minkowski} using an analytical formula. The following corollary highlights the effects of the tightened threshold on the size of the reachable set.
\begin{corollary} \label{cor_reach_shrinks}
The volume of the reachable set $\mathcal{R}_{x,t}(\bar{w}, \alpha)$ shrinks with the tightened threshold $\alpha^{\star}_{q,k}$ for all $k \geq 1$. Further, $\forall k_1, k_2 \geq 1$ and $k_1 < k_2$, the following inclusion holds $\mathcal{R}_{x,t}(\bar{w}, \alpha^{\star}_{q,k_{1}}) \supseteq \mathcal{R}_{x,t}(\bar{w}, \alpha^{\star}_{q,k_{2}})$.
\end{corollary}
\begin{proof}
Apply \eqref{eqn_reach_minkowski} with $k_1, k_2 \geq 1$ \& $k_1 < k_2$ along with $\alpha^{\star}_{q,k_{1}} > \alpha^{\star}_{q,k_{2}}$ (by Lemma \ref{lemma_1}), to see that volume$(\mathcal{R}_{x,t}(\bar{w}, \alpha^{\star}_{q,k_{2}})) \leq \text{volume}(\mathcal{R}_{x,t}(\bar{w}, \alpha^{\star}_{q,k_{1}}))$. Inclusion follows from Minkowski sum with $\alpha^{\star}_{q,k_{1}} > \alpha^{\star}_{q,k_{2}}$.
\end{proof}

\section{Numerical Simulation} \label{sec_numerical_simulation}
We consider an empirical system under study with the detector tuned to a false alarm rate $\mathcal{A} = 0.05$ ($5$\%). We demonstrate here simulation results when the uncertainties are zero-mean Gaussian. We compare the size of the reachable set boundary computed using \eqref{eqn_reach_minkowski} for thresholds $\alpha^{\star}_{q,1}, \alpha^{\star}_{q,2}, \alpha^{\star}_{q,4}$. We assume that the modeler is unaware of the functional form of the uncertainties and has to arrive at a detector threshold satisfying the desired false alarm rate. The SDP in \eqref{eqn_sdp_higher_moments} was solved with $\epsilon = 10^{-4}$ using SOSToolbox in Matlab with SeDuMi solver.
\begin{equation*}
    \begin{aligned}
        A &= \begin{bmatrix} 0.84 & 0.23 \\ -0.47 & 0.12 \end{bmatrix}, B = \begin{bmatrix} 0.07 & -0.32 \\ 0.23 & 0.58\end{bmatrix}, C = \begin{bmatrix} 1 & 0 \\ 2 & 1 \end{bmatrix} \\
        K &= \begin{bmatrix} 1.404 & -1.402 \\ 1.842 & 1.008 \end{bmatrix}, L = \begin{bmatrix} 0.0276 & 0.0448 \\ -0.01998 & -0.0290 \end{bmatrix}, \\
        \Sigma_v &= 6 I_p, \Sigma_w = \begin{bmatrix} 0.0225 & -0.0055 \\ -0.0055 & 0.0100 \end{bmatrix}.
    \end{aligned}
\end{equation*}
\begin{figure}
    \centering
    \includegraphics[scale=0.18]{sim_normal_reachable.png}
    \caption{The reachable sets associated with the detector thresholds $\alpha^{\star}_{q,1}, \alpha^{\star}_{q,2}, \alpha^{\star}_{q,4}, \alpha_{\chi^{2}}$ are shown in orange, green, blue and magenta colors respectively. It is evident from the size of the reachable set that given a desired false alarm rate $\mathcal{A}$, the knowledge of higher-order moments results in a tightened detector threshold and thereby restricting the attacker's ability to launch a larger attack.}
    \label{fig_sim_normal_reach}
\end{figure}
\begin{figure}
    \centering
    \includegraphics[scale=0.18]{sim_normal_indicate.png}
    \caption{The moment based polynomials in orange, green and blue bounding their respective indicator functions  $\mathbf{1}_{(\alpha^{\star}_{q,1},\infty)}, \mathbf{1}_{(\alpha^{\star}_{q,2},\infty)}, \mathbf{1}_{(\alpha^{\star}_{q,4},\infty)}$ in shaded orange, green and blue colors are shown here. Clearly, $\alpha^{\star}_{q,1}$ is very conservative and with $k = 4$, the threshold $\alpha^{\star}_{q,4}$ starts getting closer to the true threshold $\alpha_{\chi^{2}}$ given by \eqref{eqn_chi_squared_threshold}.}
    \label{fig_sim_normal_indicate}
\end{figure}

When the noises $w_t$ and $v_t$ are truly Gaussian, it is evident from Fig. \ref{fig_sim_normal_reach} that the reachable set corresponding to the thresholds $\alpha^{\star}_{q,1} = \alpha^{\star}_{r,2} = 40$ is conservative and ensures that the false alarm does not exceed 5\% but this also provides the attacker with the ability to launch a larger attack. However, with the knowledge of additional moments, the detector thresholds $\alpha^{\star}_{q,2} = 10.7684$, $\alpha^{\star}_{q,4} = 9.1315$ get tightened with $\alpha^{\star}_{q,4} \leq \alpha^{\star}_{q,2} \leq \alpha^{\star}_{q,1}$ as shown in Fig. \ref{fig_sim_normal_indicate}. Further, this threshold tightening limits the attacker's ability to launch a larger attack which is depicted through the reachable sets corresponding to the thresholds $\alpha^{\star}_{q,1}, \alpha^{\star}_{q,2}, \alpha^{\star}_{q,4}$ as shown in Fig. \ref{fig_sim_normal_reach}. Subsequently, the false alarm rate corresponding to the threshold $\alpha_{\chi^{2}}$ was 5\% as expected and with the thresholds $\alpha^{\star}_{q,4}, \alpha^{\star}_{q,2}, \alpha^{\star}_{q,1}$, it dropped to $1\%, 0.45\%, 0\%$ respectively. When noises are truly multi-variate Laplacian (which has heavier tails than normal distribution with same mean and covariance), it resulted in thresholds $\alpha^{\star}_{q,1} = 39.83$, $\alpha^{\star}_{q,2} = 17.23, \alpha^{\star}_{q,4} = 16.54$ with false alarm rates $0\%, 0.9\%, 1\%$ respectively. Thus, no matter what distributions satisfying \eqref{eqn_ambig_set_w_high}, \eqref{eqn_ambig_set_v_high} govern the noises $w_t$, $v_t$ respectively, the inclusion of higher-order moments restricts the attacker's potential impact through a tightened detector threshold.

\section{Conclusion \& Future Outlook}\label{sec_conclusion}
We have proposed a distributionally robust approach to form the $k$ moments based ambiguity set for the detection measure data and used it to tune the anomaly detectors for a desired false alarm rate. We found a detector threshold which guaranteed that the false alarm rate did not exceed a desired value using a semidefinite program. We have demonstrated the effectiveness of our proposed approach with a numerical example. Further, our approach using higher-order moments restricted the attacker's potential impact. Future works include addressing the problems associated with the data-driven formulation with moment estimation uncertainty and securing nonlinear cyberphysical systems with distributionally robust unscented Kalman filter based state estimation.

\bibliographystyle{IEEEtran}