Long-tailed distributions of inter-event times as mixtures of exponential distributions

TL;DR

This paper investigates the long-tailed distributions of human inter-event times, proposing that they can be modeled as mixtures of exponential distributions, and introduces a novel model selection method based on the minimum description length principle.

Contribution

It introduces a new approach using the minimum description length principle for selecting mixture models of exponential distributions, addressing challenges due to non-identifiability.

Findings

Mixtures with few components are preferred over complex ones.

Fitting accuracy is comparable to power-law distribution models.

Results support Poissonian explanations of human behavior.

Abstract

Inter-event times of various human behavior are apparently non-Poissonian and obey long-tailed distributions as opposed to exponential distributions, which correspond to Poisson processes. It has been suggested that human individuals may switch between different states in each of which they are regarded to generate events obeying a Poisson process. If this is the case, inter-event times should approximately obey a mixture of exponential distributions with different parameter values. In the present study, we introduce the minimum description length principle to compare mixtures of exponential distributions with different numbers of components (i.e., constituent exponential distributions). Because these distributions violate the identifiability property, one is mathematically not allowed to apply the Akaike or Bayes information criteria to their maximum likelihood estimator to carry out model selection. We overcome this theoretical barrier by applying a minimum description principle to joint likelihoods of the data and latent variables. We show that mixtures of exponential distributions with a few components are selected as opposed to more complex mixtures in various data sets and that the fitting accuracy is comparable to that of state-of-the-art algorithms to fit power-law distributions to data. Our results lend support to Poissonian explanations of apparently non-Poissonian human behavior.

Tables

Table 1. Table 1: Properties of the data sets. N 𝑁 N represents the number of individuals with at least one inter-event time (therefore, at least two events). The edges are directed in the Bitcoin, Email, College, and Sexual data sets. For these data sets, we only considered the individuals as sender but not recipient of the edges. N ≥ 100 subscript 𝑁 absent 100 N_{\geq 100} represents the number of individuals with at least 100 inter-event times. The mean and standard deviation (abbreviated as std) are for inter-event times and calculated on the basis of the individuals having at least 100 inter-event times. For the Sexual data set, N ≥ 100 subscript 𝑁 absent 100 N_{\geq 100} , the mean, and standard deviation are based on the individuals having at least 50 inter-event times.

Data	$N$	$N_{\ge 100}$	Time resolution	mean $\pm$ std of the inter-event time
Office	92	30	20 sec	$1.17\times {10}^3\pm 2.56\times {10}^3$ sec
Hypertext	112	72	20 sec	$3.91\times {10}^2\pm 1.11\times {10}^3$ sec
Reality	104	90	1.5 hour	$0.50\pm 1.85$ day
Bitcoin	2,969	35	128 sec	$5.05\times {10}^5\pm 1.58\times {10}^7$ sec
Email	777	415	1 sec	$1.26\times {10}^6\pm 9.49\times {10}^6$ sec
College	1,176	159	1 sec	$4.18\times {10}^4\pm 2.55\times {10}^5$ sec
Sexual	5,660	56	1 day	$16.9\pm 24.8$ day

Table 2. Table 2: Number of components of the EMM selected by each model selection criterion. Each entry of the table shows the selected ( k , k ∗ ) 𝑘 superscript 𝑘 (k,k^{*}) pair. For each of the seven data sets, the individual with the largest number of inter-event times is used. The number of inter-event times for that individual is shown in the second column of the table. We compared the EMMs with k = 1 𝑘 1 k=1 , 2 2 2 , …, 9 9 9 , 10 10 10 , 20 20 20 , 50 50 50 , and 100 100 100 under each criterion. For the Sexual data set, k = 1 𝑘 1 k=1 , 2, or 3 were ties when AIC LVC subscript AIC LVC \text{AIC}_{\text{LVC}} , BIC LVC subscript BIC LVC \text{BIC}_{\text{LVC}} , NML LVC subscript NML LVC \text{NML}_{\text{LVC}} , or DNML was used. However, these four k 𝑘 k values effectively produced the same exponential distribution with k ∗ = 1 superscript 𝑘 1 k^{*}=1 . Therefore, in the table, we wrote k = 1 𝑘 1 k=1 for these four criteria. Note that, for any criterion, the same k ∗ superscript 𝑘 k^{*} value induced by different initial k 𝑘 k values may yield different criterion values because the estimated EMM may be different between the different k 𝑘 k values even if the k ∗ superscript 𝑘 k^{*} value is the same. On the other hand, if different criteria are maximized at the same k 𝑘 k value, they are maximized by the same EMM in addition that they share the k ∗ superscript 𝑘 k^{*} value.

Data	$n$	AIC	BIC	${\text{AIC}}_{\text{LVC}}$	${\text{BIC}}_{\text{LVC}}$	${\text{NML}}_{\text{LVC}}$	DNML
Office	403	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$
Hypertext	659	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$
Reality	1,198	$(3,3)$	$(3,3)$	$(2,2)$	$(2,2)$	$(2,2)$	$(2,2)$
Bitcoin	604	$(5,4)$	$(3,3)$	$(4,3)$	$(4,3)$	$(4,3)$	$(4,3)$
Email	3,948	$(8,6)$	$(8,6)$	$(4,4)$	$(4,4)$	$(4,4)$	$(4,4)$
College	1,090	$(8,6)$	$(4,4)$	$(3,3)$	$(3,3)$	$(3,3)$	$(3,3)$
Sexual	118	$(1,1)$	$(1,1)$	$(1,1)$	$(1,1)$	$(1,1)$	$(1,1)$

Table 3. Table 3: Likelihood of the entire and truncated data sets. The results for the individual with the largest number of inter-event times in each data set are shown. The number of inter-event times larger than min i = 1 n ⁡ τ i superscript subscript 𝑖 1 𝑛 subscript 𝜏 𝑖 \min_{i=1}^{n}\tau_{i} and b ^ ^ 𝑏 \hat{b} as estimated by the PLFit is denoted by n ′′ superscript 𝑛 ′′ n^{\prime\prime} and n ′ superscript 𝑛 ′ n^{\prime} , respectively. The columns labeled all, τ > min i ⁡ τ i 𝜏 subscript 𝑖 subscript 𝜏 𝑖 \tau>\min_{i}\tau_{i} , and τ ≥ b ^ 𝜏 ^ 𝑏 \tau\geq\hat{b} are the likelihood values for the n 𝑛 n , n ′′ superscript 𝑛 ′′ n^{\prime\prime} , and n ′ superscript 𝑛 ′ n^{\prime} inter-event times, respectively.

all $\tau >\min {\tau }_i$ Data $n$ EMM Pareto $n^{\prime \prime }$ EMM Pareto Office 403 $-2780.8$ $-2804.4$ 387(96%) $-2696.3$ $-2744.0$ Hypertext 659 $-3826.3$ $-3534.8$ 464(70%) $-2924.6$ $-2916.4$ Reality 1,198 $843.4$ $2063.4$ 447(37%) $-341.0$ $-420.8$ Bitcoin 604 $-7968.0$ $-8519.7$ 570(94%) $-7680.8$ $-8291.7$ Email 3,948 $-36058.3$ $-38853.4$ 3,915(99%) $-35893.1$ $-38789.6$ College 1,090 $-8621.8$ $-9073.1$ 1,089($\approx$100%) $-8617.5$ $-9071.4$ Sexual 118 $-396.9$ $-423.4$ 108(92%) $-372.3$ $-416.8$

$\tau \ge \hat{b}$ Data $n^{\prime }$ EMM Pareto PLFit Office 296(73%) $-2196.1$ $-2278.2$ $-2056.7$ Hypertext 337(51%) $-2301.1$ $-2351.7$ $-2034.6$ Reality 116(10%) $-346.6$ $-448.3$ $-9.4$ Bitcoin 85(14%) $-1387.0$ $-1438.8$ $-1202.0$ Email 253(6%) $-3799.7$ $-3809.0$ $-3075.2$ College 983(90%) $-8140.2$ $-8616.6$ $-8021.2$ Sexual 75(64%) $-283.4$ $-334.3$ $-237.2$

Table 4. Table 4: Relationship between the best k ∗ superscript 𝑘 k^{*} value and the number of inter-event times (i.e., n 𝑛 n ) across the individuals within each data set. We calculated the Pearson correlation coefficient, r 𝑟 r , between k ∗ superscript 𝑘 k^{*} and n 𝑛 n , and its p 𝑝 p value. The calculation was based on the individuals having at least 100 inter-event times unless we state the threshold value in the first column. The cells without results are the cases in which all the individuals yielded the same k ∗ superscript 𝑘 k^{*} value and hence one cannot calculate the correlation coefficient.

Data		AIC	BIC	${\text{AIC}}_{\text{LVC}}$	${\text{BIC}}_{\text{LVC}}$	${\text{NML}}_{\text{LVC}}$	DNML
Office	$r$	$-0.077$	$-0.078$	—	—	—	—
	$p$	$6.9\times {10}^{-1}$	$6.8\times {10}^{-1}$	—	—	—	—
Hypertext	$r$	$0.005$	$-0.034$	$-0.004$	$-0.004$	—	—
	$p$	$9.6\times {10}^{-1}$	$7.8\times {10}^{-1}$	$9.7\times {10}^{-1}$	$9.7\times {10}^{-1}$	—	—
Reality	$r$	$-0.049$	$-0.097$	$-0.015$	$-0.018$	$-0.036$	$-0.036$
	$p$	$6.5\times {10}^{-1}$	$3.6\times {10}^{-1}$	$8.9\times {10}^{-1}$	$8.7\times {10}^{-1}$	$7.4\times {10}^{-1}$	$7.4\times {10}^{-1}$
Bitcoin	$r$	$0.230$	$0.011$	$0.041$	$0.006$	$0.152$	$0.152$
	$p$	$1.8\times {10}^{-1}$	$9.5\times {10}^{-1}$	$8.2\times {10}^{-1}$	$9.7\times {10}^{-1}$	$3.8\times {10}^{-1}$	$3.8\times {10}^{-1}$
Email	$r$	$0.490$	$0.470$	$0.239$	$0.226$	$0.229$	$0.243$
	$p$	$1.8\times {10}^{-26}$	$3.5\times {10}^{-24}$	$8.0\times {10}^{-7}$	$3.5\times {10}^{-6}$	$2.4\times {10}^{-6}$	$5.2\times {10}^{-7}$
College	$r$	$0.392$	$0.345$	$0.106$	$0.101$	$0.110$	$0.110$
	$p$	$3.2\times {10}^{-7}$	$8.4\times {10}^{-6}$	$1.8\times {10}^{-1}$	$2.0\times {10}^{-1}$	$1.7\times {10}^{-1}$	$1.7\times {10}^{-1}$
Sexual	$r$	$0.035$	$-0.049$	$-0.012$	$-0.012$	$0.162$	$0.099$
($\ge 50$)	$p$	$8.0\times {10}^{-1}$	$7.2\times {10}^{-1}$	$9.3\times {10}^{-1}$	$9.3\times {10}^{-1}$	$2.3\times {10}^{-1}$	$4.7\times {10}^{-1}$
Email	$r$	$0.454$	$0.331$	$0.073$	$0.061$	$0.022$	$0.058$
($\ge 500$)	$p$	$2.7\times {10}^{-7}$	$2.7\times {10}^{-4}$	$4.3\times {10}^{-1}$	$5.1\times {10}^{-1}$	$8.1\times {10}^{-1}$	$5.3\times {10}^{-1}$
College	$r$	$0.313$	$0.261$	$0.004$	$0.002$	$0.011$	$0.011$
($\ge 200$)	$p$	$1.4\times {10}^{-2}$	$4.2\times {10}^{-2}$	$9.8\times {10}^{-1}$	$9.9\times {10}^{-1}$	$9.3\times {10}^{-1}$	$9.3\times {10}^{-1}$

Equations

$$p(\tau ;\mathbf{\pi },\mathbf{\mu })=\sum _{j=1}^k\frac{{\pi }_j}{{\mu }_j}\exp \left(-\frac{\tau }{{\mu }_j}\right),$$

$${\gamma }_{ij}^{(t)}=\frac{{\pi }_j^{(t)}p({\tau }_i;{\mu }_j^{(t)})}{{\sum }_{j^{\prime }=1}^k{\pi }_{j^{\prime }}^{(t)}p({\tau }_i;{\mu }_{j^{\prime }}^{(t)})},$$

$${\pi }_j^{(t+1)}$$

$${\mu }_j^{(t+1)}$$

$${\hat{z}}_i$$

$${\hat{\pi }}_j(\hat{}\mathbf{z})=\frac{n_j}{n},$$

$${\hat{\mu }}_j(\mathbf{\tau },\hat{}\mathbf{z})=\frac{1}{n_j}\sum _{i=1:{\hat{z}}_i=j}^n{\tau }_i.$$

$$\text{AIC}(\mathbf{\tau })=-\log p(\mathbf{\tau };{\hat{\theta }}_{\mathrm{EM}}(\mathbf{\tau }))+2k-1$$

$$\text{BIC}(\mathbf{\tau })=-\log p(\mathbf{\tau };{\hat{\theta }}_{\mathrm{EM}}(\mathbf{\tau }))+\frac{2k-1}{2}\log n,$$

$${\text{AIC}}_{\text{LVC}}(\mathbf{\tau },\hat{}\mathbf{z})=-\log p(\mathbf{\tau },\hat{}\mathbf{z};\hat{\theta }(\mathbf{\tau },\hat{}\mathbf{z}))+2k^{\ast }-1$$

$${\text{BIC}}_{\text{LVC}}(\mathbf{\tau },\hat{}\mathbf{z})=-\log p(\mathbf{\tau },\hat{}\mathbf{z};\hat{\theta }(\mathbf{\tau },\hat{}\mathbf{z}))+\frac{k^{\ast }-1}{2}\log n+\frac{1}{2}\sum _{j=1}^{k^{\ast }}\log n_j.$$

$$p(\mathbf{\tau },\hat{}\mathbf{z};\hat{\theta }(\mathbf{\tau },\hat{}\mathbf{z}))=\prod _{j=1}^{k^{\ast }}{\left({\hat{\pi }}_j\right)}^{n_j}p(\{{\tau }_i{\}}_{i=1,{\hat{z}}_i=j}^n;{\hat{\mu }}_j)=\prod _{j=1}^{k^{\ast }}{\left(\frac{{\hat{\pi }}_j}{e{\hat{\mu }}_j}\right)}^{n_j}.$$

$$\sum _{\mathbf{x}}p(\mathbf{x};\hat{\theta }(\mathbf{x}))>1$$

$$p_{\mathrm{NML}}(\mathbf{x})=\frac{p(\mathbf{x};\hat{\theta }(\mathbf{x}))}{{\sum }_{{\mathbf{x}}^{\prime }}p({\mathbf{x}}^{\prime };\hat{\theta }({\mathbf{x}}^{\prime }))}$$

$$L_{\mathrm{NML}}(\mathbf{x})$$

$$C(n)\equiv \sum _{{\mathbf{x}}^{\prime }}p({\mathbf{x}}^{\prime };\hat{\theta }({\mathbf{x}}^{\prime }))$$

$$C(n)=\int p({\mathbf{x}}^{\prime };\hat{\theta }({\mathbf{x}}^{\prime }))\mathrm{d}{x}_1^{\prime }\cdots \mathrm{d}{x}_n^{\prime }.$$

$$p(\tau ;\mu )=\frac{1}{\mu }\exp \left(-\frac{\tau }{\mu }\right).$$

$$L_{\text{NML}}(\mathbf{\tau })$$

$$C_{\exp }(n)$$

$$\hat{\mu }=\frac{1}{n}\sum _{i=1}^n{\tau }_i^{\prime }$$

$$\log p(\mathbf{\tau };\hat{\mu }(\mathbf{\tau }))$$

$$\log C_{\exp }(n)$$

$$L_{\mathrm{NML}}(\mathbf{\tau })=n\log \hat{\mu }+n\log n-\log \Gamma (n)+\log (m_{\max }-m_{\min }).$$

$${\tilde{L}}_{\mathrm{NML}}(\mathbf{\tau })=L_{\mathrm{NML}}(\mathbf{\tau })+\ell (m_{\min })+\ell (m_{\max }).$$

$$L_{\text{LVC}}(\mathbf{\tau },\hat{}\mathbf{z})=-\log p(\mathbf{\tau },\hat{}\mathbf{z};\hat{\theta }(\mathbf{\tau },\hat{}\mathbf{z}))+\log C_{\mathrm{EMM}}(n,k^{\ast }),$$

$$C_{\mathrm{EMM}}(n,k^{\ast })=\sum _{{\hat{z}}_1^{\prime }=1}^{k^{\ast }}\cdots \sum _{{\hat{z}}_n^{\prime }=1}^{k^{\ast }}{\int }_{{\mathcal{T}}^{\prime }({\mu }_{\min }^{\prime },{\mu }_{\max }^{\prime })}p({\mathbf{\tau }}^{\prime },\hat{}{\mathbf{z}}^{\prime };\hat{\theta }({\mathbf{\tau }}^{\prime },\hat{}{\mathbf{z}}^{\prime }))\mathrm{d}{\tau }_1^{\prime }\cdots \mathrm{d}{\tau }_n^{\prime }.$$

$${\mu }_{\min }^{\prime }=\exp (m_{\min }^{\prime }),$$

$$m_{\min}^{\prime}=\Bigl{\lfloor}\log\left(\min_{j=1,\ldots,k^{*}}\hat{\mu}_{j}\right)\Bigr{\rfloor},$$

$${\mu }_{\max }^{\prime }=\exp (m_{\max }^{\prime }),$$

$$m_{\max}^{\prime}=\Bigl{\lceil}\log\left(\max_{j=1,\ldots,k^{*}}\hat{\mu}_{j}\right)\Bigr{\rceil},$$

Full text

Long-tailed distributions of inter-event times as mixtures of exponential distributions

Makoto Okada1, Kenji Yamanishi1, and Naoki Masuda2,3,4,∗  

 

1Graduate School of Information Science and Technology,

The University of Tokyo, 7-3-1 Hongo Bunkyo, Tokyo 113-8656, Japan

2Department of Engineering Mathematics,

University of Bristol, Woodland Road, Clifton, Bristol BS8 1UB, United Kingdom

3Department of Mathematics

University at Buffalo, State University of New York, Buffalo, NY 14260-2900, USA

4 Computational and Data-Enabled Science and Engineering Program,

University at Buffalo, State University of New York, Buffalo, NY 14260-5030, USA

 

∗ Corresponding author (naokimas@buffalo.edu)

Abstract

Inter-event times of various human behavior are apparently non-Poissonian and obey long-tailed distributions as opposed to exponential distributions, which correspond to Poisson processes. It has been suggested that human individuals may switch between different states in each of which they are regarded to generate events obeying a Poisson process. If this is the case, inter-event times should approximately obey a mixture of exponential distributions with different parameter values. In the present study, we introduce the minimum description length principle to compare mixtures of exponential distributions with different numbers of components (i.e., constituent exponential distributions). Because these distributions violate the identifiability property, one is mathematically not allowed to apply the Akaike or Bayes information criteria to their maximum likelihood estimator to carry out model selection. We overcome this theoretical barrier by applying a minimum description principle to joint likelihoods of the data and latent variables. We show that mixtures of exponential distributions with a few components are selected as opposed to more complex mixtures in various data sets and that the fitting accuracy is comparable to that of state-of-the-art algorithms to fit power-law distributions to data. Our results lend support to Poissonian explanations of apparently non-Poissonian human behavior.

1 Introduction

Many social and economic processes are a consequence of human behavior. Technological advances are increasingly enabling us to record various human behaviors in quantitative manners. Such quantitative understanding often challenges traditional assumptions underlying mathematical modeling of human behavior, a major instance of which is a lack of Poissonian properties in a range of human behavioral data. In other words, when a sequence of time-stamped events, such as email correspondences, is observed from a human individual, the event sequence often deviates from Poisson processes, in which inter-event times independently obey an exponential distribution. Rather, empirical inter-event times from human behavior and other data often obey long-tailed distributions [1, 2]. Such non-Poissonian processes on nodes and edges are building blocks of temporal networks on which dynamical processes such as epidemic processes often behave differently from the same processes on static networks [3, 4].

Two classes of models that generate long-tailed distributions of inter-event times are priority queue models and modulated Markov processes [3, 4]. In priority queue models, one assumes that a human individual is a queue that receives tasks across time and prioritizes some particular tasks for execution [1]. In contrast, in modulated Markov processes, one assumes that a human individual generates events according to a Poisson process (i.e., one determines whether to have an event right now at a constant rate without memory) but modulates the event rate. The event rate may be assumed to take one of a few values [5, 6, 7, 8, 9, 10] or continuously many values [11, 12]. When we assume a few states, where each state corresponds to a Poisson process, an underlying assumption about human behavior is that an individual transits among a few distinguishable states, such as an active state and an inactive state. If this is the case, we should be able to fit a mixture of exponential distributions rather than power-law distributions to empirical distributions of inter-event times with a reasonable accuracy. This is because a subset of inter-event times, i.e., those produced in each state, is expected to obey an exponential distribution.

In fact, a mixture of a small number of exponential distributions and a power-law distribution with an exponential cutoff can look similar. In the present study, we fit mixtures of exponential distributions to several data sets of inter-event times. The idea of approximating empirical distributions resembling power-law distributions by mixtures of exponential distributions is not a new idea [13]. We perform model selection to determine the number of component exponential distributions fitted to data and also compare the mixture of exponential distributions with two types of conventionally used power-law distributions. A technical challenge is that the maximum likelihood estimator of the mixture of exponential distributions does not satisfy the asymptotic normality such that it is not allowed to use the AIC (Akaike information criterion) or BIC (Bayes information criterion) (e.g. p. 203 in Ref. [14]). Therefore, based on the minimum description length criterion [15, 16], we perform the procedure so-called latent variable completion [17, 18, 19, 20] to derive variants of minimum description length justified for mixtures of exponential distributions. We find that mixtures of up to three exponential distributions are the best performer in many cases. Python codes for estimating and selecting the mixture of exponential distributions are available at Github (https://github.com/naokimas/exp_mixture_model).

2 Methods

2.1 Mixture of exponential distributions and its maximum likelihood estimation

We fit mixtures of exponential distributions, which we refer to them as the exponential mixture models, abbreviated as EMMs, to distributions of inter-event times, denoted by $\tau\in[0,\infty)$, and compare the fit with the case of power-law distributions. We denote by $k$ the number of the exponential distributions to be mixed, which we refer to as the number of component distributions. The mixing weight, i.e., the probability that the $j$th exponential distribution is used, is denoted by $\pi_{j}$ ($1\leq j\leq k$). The probability density function for an EMM is given by

[TABLE]

where $\bm{\pi}\equiv\{\pi_{1},\ldots,\pi_{k}\}$, $\bm{\mu}\equiv\{\mu_{1},\ldots,\mu_{k}\}$, and $\mu_{j}$ is the mean of the $j$th exponential distribution.

Given a series of inter-event times, $\bm{\tau}\equiv\{\tau_{1},\ldots,\tau_{n}\}$, and the number of components, $k$, we estimate the values of the model parameters, $\theta=\left(\bm{\pi},\bm{\mu}\right)$, as follows.

First, we run the expectation-maximization (EM) algorithm [22]. As an initial condition, we set $\pi^{(0)}_{j}=1/k$ ($1\leq j\leq k$), where the superscript indicates the iteration number in the EM procedure. We also draw $\mu^{(0)}_{j}$ ($1\leq j\leq k$) such that $\log_{10}\mu^{(0)}_{j}$ independently obeys the uniform density on $[\log_{10}\tau_{\min},\log_{10}\tau_{\max}]$, where $\tau_{\min}$ and $\tau_{\max}$ are the minimum and maximum of the inter-event time in the given data set, respectively. In this manner, we intend to generate a sufficiently broad initial distribution of $\mu_{j}^{(0)}$ while a majority of the $\mu^{(0)}_{j}$ values is relatively small. Intuitively, a large value of $\mu_{j}^{(0)}$, which sometimes occurs, reflects the long tail of the distribution of inter-event times. We iterate the following EM steps $t_{\max}$ times by alternating the expectation (E) and maximization (M) steps. The E step is given by

[TABLE]

where $t=0,1,2,\ldots,t_{\max}-1$. One can interpret $\gamma_{ij}^{(t)}$ as the probability that $\tau_{i}$ is generated from the $j$th exponential distribution. The M step is given by

[TABLE]

After the iteration of the E and M steps $t_{\max}$ times, one obtains $\pi^{(t_{\max})}_{1},\ldots,\pi^{(t_{\max})}_{k}$ and $\mu^{(t_{\max})}_{1},\ldots,\mu^{(t_{\max})}_{k}$, which is an approximate maximum likelihood estimator. This estimator corresponds to the case in which one does not estimate the latent variables that explicitly indicate which exponential distribution out of the $k$ exponential distributions has generated each $\tau_{i}$.

When the values of the latent variables are required, we compute them as follows. We denote the latent variable corresponding to each $\tau_{i}$ ($1\leq i\leq n$) by $z_{i}$ ($1\leq z_{i}\leq k$). In other words, $\tau_{i}$ is estimated to be generated from the $z_{i}$th exponential distribution. We first estimate the latent variables by

[TABLE]

We then estimate the other parameters by

[TABLE]

where $\bm{\hat{}}{\bm{z}}=\{\hat{z}_{1},\ldots,\hat{z}_{n}\}$, $n_{j}$ is the number of inter-event times whose estimated latent variable $\hat{z}_{i}$ is equal to $j$, and

[TABLE]

The $\hat{\pi}_{j}$ and $\hat{\mu}_{j}$ values given by Eqs. (5) and (6) are a maximum likelihood estimator when the joint distribution of the inter-event times (i.e., $\tau_{1}$, $\ldots$, $\tau_{n}$) and the latent variables (i.e., $\hat{z}_{1}$, $\ldots$, $\hat{z}_{n}$) is estimated. They are different from the values estimated by the EM algorithm (i.e., $\pi^{(t_{\max})}_{1},\ldots,\pi^{(t_{\max})}_{k}$ and $\mu^{(t_{\max})}_{1},\ldots,\mu^{(t_{\max})}_{k}$).

To cope with the problem of local maxima, we estimate the values of $\hat{\pi}_{j}$, $\hat{\mu}_{j}$, and $\hat{z}_{i}$ ($1\leq j\leq k$, $1\leq i\leq n$) ten times starting from different initial conditions. Then, among the ten maximum likelihood estimators, we employ the one that has yielded the largest likelihood of the joint distribution of $\bm{\tau}$ and $\hat{\bm{z}}$.

2.2 Model selection criteria for mixtures of exponential distributions

We consider the following six model selection criteria.

2.2.1 AIC and BIC

For an EMM with $k$ components without an explicit consideration of the latent variables, the AIC [23] and the BIC [24] are given by

[TABLE]

and

[TABLE]

respectively. In Eqs. (7) and (8), $\hat{\theta}_{\rm EM}(\bm{\tau})$ is the EM estimator of the parameters of the EMM, i.e., $\pi^{(t_{\max})}_{1},\ldots,\pi^{(t_{\max})}_{k},\mu^{(t_{\max})}_{1},\ldots,\mu^{(t_{\max})}_{k}$.

2.2.2 AIC and BIC with latent variable completion

An EMM is nonidentifiable, which by definition dictates that a distribution does not correspond to a parameter set in a one-to-one manner [14]. For example, if $k=2$ and $\mu_{1}=\mu_{2}$, the distribution of $\tau$ is the same exponential distribution regardless of the value of $\pi_{1}(=1-\pi_{2})$. When a distribution is nonidentifiable, the maximum likelihood estimator is not asymptotically normal, and the conventional AIC and BIC, which are derived on the basis of asymptotic normality, lose justification [14].

A method for overcoming this theoretical barrier is to apply model selection criteria to a complete variable model, in which one completes the values of the latent variables, rather than to use a marginalized model [17, 18, 19, 20]. We call this method the latent variable completion [20] throughout this paper. Under latent variable completion, one explicitly incorporates the likelihood of the latent variables, $z_{i}$ $(1\leq i\leq n)$, without marginalizing them, into a model selection criterion. Then, the joint distribution of $\tau_{i}$ and $z_{i}$ $(1\leq i\leq n)$ is an identifiable distribution such that the use of the AIC and BIC is justified. Also in practice, latent variable completion is better at describing some empirical data (for example, at estimating an appropriate number of components, $k$) than model selection criteria that do not use latent variable completion [17, 18].

We refer to the corresponding AIC or BIC as the AIC or BIC with latent variable completion [17, 19] and denote them by $\text{AIC}_{\text{LVC}}$ and $\text{BIC}_{\text{LVC}}$, respectively. We obtain

[TABLE]

and

[TABLE]

In Eqs. (9) and (10), $\hat{\theta}(\bm{\tau},\bm{\hat{}}{\bm{z}})$ = $\{\hat{\pi}_{1}(\bm{\hat{}}{\bm{z}})$, …, $\hat{\pi}_{k^{*}}(\bm{\hat{}}{\bm{z}})$, $\hat{\mu}_{1}(\bm{\tau},\bm{\hat{}}{\bm{z}})$, …, $\hat{\mu}_{k^{*}}(\bm{\tau},\bm{\hat{}}{\bm{z}})\}$, derived in Eqs. (5) and (6), is the maximum likelihood estimator for the joint distribution of inter-event times $\bm{\tau}$ and latent variables $\bm{\hat{}}{\bm{z}}$. The likelihood in Eqs. (9) and (10) is given by

[TABLE]

In Eqs. (9), (10), and (11), $k^{*}$ is the number of components that are used at least once under the estimated latent variable values. The latent variable specifies which of the $k$ exponential distributions has produced $\tau_{i}$ for each $i$. Therefore, as a result of estimating the latent variables, a $j$th exponential distribution may not have any inter-event times $\tau_{i}$ belonging to it (i.e., $n_{j}=0$). We exclude such exponential distributions, i.e., those with $n_{j}=0$, and only consider the remaining $k^{*}$ exponential distributions. To calculate $\text{AIC}_{\text{LVC}}$ and $\text{BIC}_{\text{LVC}}$, we use $k^{*}$ instead of $k$ because $k^{*}$ is the actual number of components given the estimated latent variable values. For the same reason, we will use $k^{*}$ instead of $k$ in the two model selection criteria introduced in the following sections, which also explicitly use the estimated latent variables and latent variable completion.

2.2.3 Normalized maximum likelihood codelength

Here we introduce model selection criteria on the basis of the minimum description length (MDL) principle [16]. This principle asserts that the best model should be the one that attains the shortest codelength required for encoding the data as well as the model itself. The codelength is the number of bits required for encoding the data into a binary sequence under the information-theoretic requirement that each codeword can be uniquely decodable even without commas. In fact, MDL approaches to mixture modeling have a long history, as represented by the Snob program [25, 26].

The normalized maximum likelihood (NML) codelength is a type of minimum description length [27, 16]. The NML codelength minimizes the worst-case (i.e., in terms of data $\bm{x}$) regret value, which is equal to the the actual codelength minus the ideal codelength, i.e., $\min_{\theta\in\Theta}\left(-\log p(\bm{x};\theta)\right)$ [28]. The other two model selection criteria, which we will explain in the next two sections, are based on the NML. Therefore, we explain the NML codelength in this section. Although we assume that inter-event times $\tau_{i}$ are continuous-valued, we first explain the NML codelength for a sequence of discrete-valued variables $\bm{x}=\{x_{1},\ldots,x_{n}\}$ for expository purposes.

To shorten the codelength for the given data $\bm{x}$, one should in principle minimize the Shannon information given by $-\log p(\bm{x};\theta)$, where we assume throughout the present paper that $\log$ is in base $e$ unless we specify the base. The minimizer is obviously given by the maximum likelihood estimator, $\hat{\theta}(\bm{x})$. However, the maximum likelihood estimator generally yields

[TABLE]

because $\hat{\theta}$ depends on $\bm{x}$. In Eq. (12), the summation is taken over all possible values of $\bm{x}$. Therefore, we use the normalized probability distribution given by

[TABLE]

to encode the data [29]. In other words, we set

[TABLE]

where

[TABLE]

is called the parametric complexity. In Eq. (14), the first term on the right-hand side is small when the model fits the data well, and the second term represents the complexity of the model. In general, $C(n)$ tends to increase as $k$ increases because an increase in $k$ leads to the expansion of the parameter space in which one searches $\hat{\theta}$, which leads to an increase in $C(n)$ [30].

One can similarly derive the NML codelength for a sequence of continuous-valued variables by replacing the summation by the integral. For example, the equivalent of Eq. (15) in the case of continuous-valued variables is given by

[TABLE]

In the remainder of this section, we explain the NML codelength for an exponential distribution [31, 19], not for an EMM, for two reasons. First, analytically calculating the NML codelength for an EMM seems formidable. Second, we will use the NML codelength for single exponential distributions in deriving the two types of the NML-based codelengths for EMMs in the following sections.

Consider an exponential distribution given by

[TABLE]

For this distribution, we obtain

[TABLE]

where $\mathcal{T}(\mu_{\min},\mu_{\max})$ is a region of $\bm{\tau}^{\bm{\prime}}\equiv\{\tau_{1}^{\prime},\ldots,\tau_{n}^{\prime}\}\in\mathbf{R}^{n}$ such that the maximum likelihood estimator given by

[TABLE]

is contained in $[\mu_{\min},\mu_{\max}]$. We have to confine the domain of integration to $\mathcal{T}$ in this manner because the integral in Eq. (18b) diverges if $\mathcal{T}$ is replaced by $\mathbf{R}^{n}$. In practice, if we take $\mu_{\min}$ $(>0)$ small enough and $\mu_{\max}$ large enough, $\hat{\mu}\in[\mu_{\min},\mu_{\max}]$ would be satisfied for empirical sequences of inter-event times $\bm{\tau}^{\bm{\prime}}$. The standard method to evaluate the integral in Eq. (18b) is to transform the integral to that in the parameter space [32, 19]. This method yields an integral of the form similar to $\int_{\mu_{\min}}^{\mu_{\max}}(1/\hat{\mu}){\rm d}\hat{\mu}$, which is manageable [19]. Region $\mathcal{T}$ is a useful heuristic for avoiding such a divergence in NML-related calculations [32].

To incorporate the codelength necessary for encoding the value of $\mu_{\min}$ and $\mu_{\max}$ into the NML codelength, we rewrite $\mu_{\min}=\exp(m_{\min})$ and $\mu_{\max}=\exp(m_{\max})$, where $m_{\min}$ and $m_{\max}$ are integer. Then, we encode $m_{\min}$ and $m_{\max}$ instead of $\mu_{\min}$ and $\mu_{\max}$.

The first term on the right-hand side of Eq. (18a) is given by

[TABLE]

where $\hat{\mu}$ is given by Eq. (19). To obtain the second term on the right-hand side of Eq. (18a), we substitute Eq. (17) into Eq. (18b), which leads to

[TABLE]

where $\Gamma(n)$ is the gamma function.

By substituting Eqs. (20) and (21) into Eq. (18a), one obtains

[TABLE]

To carry out model selection, we add the codelength for $m_{\min}$ and $m_{\max}$, denoted by $\ell(m_{\min})$ and $\ell(m_{\max})$, to $L_{\rm NML}(\bm{\tau})$ to obtain

[TABLE]

The derivation of $\ell(m_{\min})$ and $\ell(m_{\max})$ is given in Supplementary Materials.

2.2.4 Normalized maximum likelihood codelength with latent variable completion

As a first NML type of the model selection criterion for EMMs, we consider a latent variable completion of the NML codelength. We refer to the resulting criterion by $\text{NML}_{\text{LVC}}$ and the codelength by $L_{\text{LVC}}$. Although an analytical expression for the NML codelength for an EMM without latent variable completion is unavailable, we can analytically calculate the NML codelength for the joint distribution of inter-event times and the latent variables.

Similar to the derivation in the case of a mixture of Gaussian distributions [19], we derive the codelength for EMMs as follows:

[TABLE]

where

[TABLE]

We remind that $p(\bm{\tau}^{\bm{\prime}},\bm{\hat{}}{\bm{z}}^{\bm{\prime}};\hat{\theta}(\bm{\tau}^{\bm{\prime}},\bm{\hat{}}{\bm{z}}^{\bm{\prime}})$ is given by Eq. (11). Region $\mathcal{T}^{\prime}(\mu_{\min}^{\prime},\mu_{\max}^{\prime})$ of $\bm{\tau}^{\bm{\prime}}\in\mathbf{R}^{n}$ is defined such that each of $\hat{\mu}_{1}$, $\ldots$, $\hat{\mu}_{k^{*}}$ in the maximum likelihood estimator $\hat{\theta}(\bm{\tau}^{\bm{\prime}},\bm{\hat{}}{\bm{z}}^{\bm{\prime}})$ when the latent variables are explicitly estimated (i.e., Eq. (6)) is contained in $[\mu_{\min}^{\prime},\mu_{\max}^{\prime}]$. Region $\mathcal{T}^{\prime}(\mu_{\min}^{\prime},\mu_{\max}^{\prime})$ coincides with region $\mathcal{T}(\mu_{\min}^{\prime},\mu_{\max}^{\prime})$ when $k^{*}=1$. Similarly to the case of Eq. (18b), we confine the domain of integration to $\mathcal{T}^{\prime}$ to avoid the divergence of the integral in Eq. (25). We set

[TABLE]

where

[TABLE]

and

[TABLE]

where

[TABLE]

for two reasons. First, $\log C_{\rm EMM}(n,k^{*})$ is small when $\mu_{\min}^{\prime}$ is large or $\mu_{\max}^{\prime}$ is small. Second, encoding of $\mu_{\min}^{\prime}$ and $\mu_{\max}^{\prime}$ is facilitated by the introduction of integer variables, as we did for the NML codelength (section 2.2.3).

By substituting

[TABLE]

into Eq. (24), where $H(q_{1},\ldots,q_{k^{*}})=\sum_{j=1}^{k^{*}}-q_{j}\log q_{j}$ represents the entropy, one obtains

[TABLE]

By substituting Eq. (11) in Eq. (25), one obtains the parametric complexity as follows:

[TABLE]

where

[TABLE]

Note that we have used Eq. (21) to derive Eq. (33) and that $C_{\rm EMM}(n,1)$ coincides with $C_{\rm exp}(n)$. Because we cannot calculate $C_{\rm EMM}(n,k^{*})$ using Eq. (32) due to combinatorial explosion, we calculate $C_{\rm EMM}(n,k^{*})$ using the recursive relationship [19] given by

[TABLE]

The calculation of $C_{\rm EMM}(n,k^{*})$ using this recursive relationship, which dominates the computation time for calculating $L_{\text{LVC}}(\bm{\tau},\bm{\hat{}}{\bm{z}})$, requires $O(n^{2}k^{*})$ time.

Finally, we add the codelength to encode integers $m_{\min}^{\prime}$ and $m_{\max}^{\prime}$ to obtain

[TABLE]

2.2.5 Decomposed NML codelength

The second type of the NML codelength with latent variable completion, which we call the decomposed NML (DNML) codelength [20, 21], also completes the latent variable and then calculates the NML codelength. The DNML does so after decomposing the joint distribution into the distribution of the latent variables, $\bm{z}$, and that of the observables, $\bm{\tau}$, conditioned on the value of $\bm{z}$ [20]. The DNML codelength was originally formulated for the case in which both observables and latent variables were discrete [20]. Here we regard the inter-event time, $\tau_{i}$, as continuous-valued and the latent variable, $\hat{z}_{i}$, as discrete-valued.

We start by considering

[TABLE]

Equation (22) yields

[TABLE]

where $m_{\min}^{\prime}$ and $m_{\max}^{\prime}$ are given by Eqs. (27) and (29), respectively. It should be noted that we could prepare integers $m_{\min}$ and $m_{\max}$ for each $j$ $(1\leq j\leq k^{*})$ to code $\{\tau_{i}\}_{i=1,\hat{z}_{i}=j}^{n}$ through $\hat{\mu}_{j}$, similarly to Eq. (23), to derive an alternative of Eq. (37). However, that coding method requires $2k^{*}$ integers and therefore less efficient than using only two integers $m_{\min}^{\prime}$ and $m_{\max}^{\prime}$.

The NML codelength $L_{\text{NML}}(\bm{\hat{}}{\bm{z}})$ on the right-hand side of Eq. (36) is for a multinomial distribution of $\hat{z}_{i}$ and is given by

[TABLE]

In Eq. (38), $C_{\rm mult}(n,k^{*})$ is the parametric complexity for the multinomial distribution having $k^{*}$ elements. Using Eq. (5), one obtains

[TABLE]

where $n_{j}^{\prime}$ is the number of inter-event times whose latent variable $\hat{z}^{\prime}_{i}$ is equal to $j$. One can recursively calculate $C_{\rm mult}(n,k^{*})$ [33] by

[TABLE]

The computational time for solving the set of recursive equations is given by $O(n+k^{*})$ and dominates the computational time for $L_{\text{DNML}}(\bm{\tau},\bm{\hat{}}{\bm{z}})$.

Finally, similarly to the case of $\text{NML}_{\text{LVC}}$, we add the codelength to encode integers $m_{\min}^{\prime}$ and $m_{\max}^{\prime}$ to obtain

[TABLE]

2.3 Power-law distributions

We also fitted two types of power-law distributions to the empirical distributions of inter-event times. The first type is the Pareto distribution whose probability density function is given by

[TABLE]

which is defined for $\tau\geq b$. We use the maximum likelihood estimator given by [34, 35]

[TABLE]

and

[TABLE]

The second estimator was the one proposed by Clauset et al. [36]. In this method, which we refer to as the PLFit algorithm, one selects the $\hat{b}$ value such that the distribution of the data points satisfying $\tau_{i}\geq\hat{b}$ is as close as possible to a power law. The power-law exponent given the $\hat{b}$ value is estimated by Eq. (44). It should be noted that the PLFit does not intend to fit a power-law or different distribution to the data points whose values are less than $\hat{b}$. We used a publicly available Python implementation of the PLFit [37].

3 Results

3.1 Data

We used the following seven data sets of time-stamped event sequences obtained from human behavior. Basic properties of each data set are shown in Table 1.

First, we used the data that the SocioPatterns project collected from individuals in an office building that hosted scientific departments [38]. We refer to this data set as the Office data set. The individuals wore a sensor on their chest, with which physical proximity between pairs of individuals was detected. The recording lasted for two weeks. A time-stamped event was defined as a contact that lasted at least 20 seconds, which was the time resolution of the data. For each individual, we ignored the partner of the contact and used the beginning time of each contact as the time of the event. In this manner, we examined a sequence of time-stamped events for each individual in this and the following data sets. Note that, because a contact is a symmetric relationship between two individuals, any event between nodes $i$ and $j$ is considered for both $i$ and $j$.

The inter-event time was defined as the difference between the starting time of the two consecutive events. There was no event during the night time due to the circadian rhythm. The effect of the circadian rhythm on analysis of inter-event times may be considerable but beyond the scope of the present study. Therefore, we excluded the inter-event times when the two events belonged to different days unless otherwise stated. When there were multiple events in consecutive 20-second time windows between the same pair of individuals, we regarded that they constitute a single event lasting over 20 seconds rather than a sequence of events with inter-event times equal to 20 seconds. In practice, in this case, we discarded all but the first event in each sequence of the multiple events in consecutive 20-second time windows and then calculated the inter-event times. We also aggregated events for a focal node that occurred at the same time with different partners into one event. Therefore, the minimum inter-event time is equal to 20 seconds. For the individuals that have at least 100 inter-event times, the mean and standard deviation of the inter-event times are shown in Table 1.

Second, we used the Hypertext 2009 dynamic contact network (Hypertext for short) [39]. Under the SocioPatterns project, the data were collected from approximately 75% of the participants in the HT09 conference in Torino over three days. As was the case for the Office data set, we ignored inter-event times that span multiple days.

The third data set originates from Reality Mining Project (Reality for short), which provides time-stamped physical proximity relationships between the participants of an experiment, who are students or faculty members of MIT. The data were obtained from Bluetooth logs of mobile phones [40]. We did not skip inter-event times overnight because the data were collected over months and the mean inter-event time was much longer than that for the Office and Hypertext data sets (Table 1). We used the data from September 2004 to May 2005, which were the months in each of which there were at least $10^{3}$ events in total.

Fourth, Bitcoin OTC is an over-the-counter marketplace where users trade with bitcoin [41]. In Bitcoin OTC, users rate other users regarding the trustworthiness. We neglected the rate score and who a focal user $i$ rated. For each user $i$, we examined a series of time-stamped events, where an event was defined by rating behavior by user $i$ toward anybody. We refer to this data set as Bitcoin.

The Email-Eu-core (Email for short) data set was collected from a large European research institution between October 2003 and May 2005 [42]. A node was an email address. An edge was defined between two nodes if and only if they sent email to each other at least once. A time-stamped event for node $i$ in the network was defined by an email that $i$ sent to anybody, disregarding the identity of the recipient of the email.

The CollegeMsg (College for short) data were collected from students belonging to an online community at the University of California, Irvine, between April and December 2004 [43]. The event was defined as the message sent from a user to another user. We downloaded the Bitcoin, Email, and College data from the Stanford Network Analysis Project (SNAP) website (http://snap.stanford.edu/).

The sexual contact (Sexual for short) data set provides times of online sexual encounters between escorts and sex buyers [44]. We used event sequences for each sex buyer by regarding a commercial sexual activity with any escort as an event.

3.2 Fitting of the EMM to the individual with the largest number of events and comparison with power-law distributions

We fitted EMMs with different numbers of components, $k$, to the sequence of inter-event times obtained from each individual in each data set. We examined $k=1$, $2$, …, $9$, $10$, $20$, $50$, and $100$ in each case. Then, according to each of the six model selection criteria, AIC, BIC, $\text{AIC}_{\text{LVC}}$, $\text{BIC}_{\text{LVC}}$, $\text{NML}_{\text{LVC}}$, and DNML, we selected the best $k$ value.

Results of model selection for the individual with the largest number of events in each data set are shown in Table 2. Table 2 indicates that the effective number of components finally chosen, $k^{*}$, is at most four in all but the three out of the 42 combinations of the data set and criterion. Furthermore, for the four out of the seven data sets, $k^{*}$ is at most three regardless of the criterion. These results suggest that a mixture of a small number of exponential distributions may be a reasonable approximation to empirical distributions of inter-event times in many cases.

Next, we compared the performance between the estimated EMMs and power-law estimators in approximating the empirical distributions. We have fitted power-law distributions because they have widely been applied to empirical distributions of inter-event times [1, 2, 3, 4]. For the same individuals as those used in Table 2, the empirical and estimated distributions of inter-event times are compared in Fig. 1. When the different model selection criteria yielded different EMMs (i.e., EMMs with different values of $k$ or $k^{*}$), we showed all the selected EMMs overlaid in the figure. We compared the different distributions in terms of the survival probability (i.e., the probability with which the inter-event time is larger than $\tau$, i.e., $\int_{\tau}^{\infty}p(\tau^{\prime}){\rm d}\tau^{\prime})$ and the odds ratio defined by \text{OR}(\tau)=\left[1-\int_{\tau}^{\infty}p(\tau^{\prime}){\rm d}\tau^{\prime}\right]\big{/}\int_{\tau}^{\infty}p(\tau^{\prime}){\rm d}\tau^{\prime}. Note that the odds ratio is particularly good at capturing differences between distributions at small values of $\tau$ [45].

Figure 1 elicits the following observations. First, in five data sets (i.e., except Bitcoin and Email) the right tails of the distribution seem to be more accurately approximated by an EMM than the power-law distributions. This may be because the two power-law estimators employed here assume a strictly power-law tail, whereas empirical data usually have exponential cutoffs. For the Bitcoin and Email data sets, the power-law estimate obtained by the PLFit algorithm seems to be roughly as accurate as EMMs in approximating right tails of the distribution. Second, the odds ratio plots show that, at small $\tau$ values, the approximation looks the most accurate with the Pareto distribution in all but the College data set.

To examine whether these casual observations extend to the entire set of inter-event times from the same individuals, we investigated the likelihood of the empirical data in each model. The likelihood for the three models, i.e., the selected EMM, maximum likelihood estimator of the Pareto distribution, and PLFit is compared in Table 3. Because the PLFit discards small inter-event times, we start by comparing the EMM and Pareto distribution when all inter-event times are used. The table indicates that the EMM realizes a larger likelihood value than the Pareto distribution in five data sets and vice versa in the other two data sets (i.e., Hypertext and Reality). The EMM has more parameters than the Pareto distribution, which one may suspect as a reason why the EMM fits the data better than the Pareto distribution in more data sets than vice versa. However, the results are qualitatively the same when the two models are compared in terms of the AIC and BIC (Table S1). We remark that the maximum likelihood estimator for the Pareto distribution is not asymptotically normal such that the application of the AIC and BIC to the Pareto distribution as well as to EMMs is not justified. It should be noted that the other four criteria are not relevant to the Pareto distribution because it is not a latent variable model.

For the Hypertext and Reality data sets, for which the likelihood is larger for the Pareto distribution than the optimal EMM, there are relatively many data points with $\tau_{i}=\tau_{\min}$ (Table 3). The Pareto distribution beats the EMM probably because it avoids devoting the probability mass to values less than $\tau_{\min}$ by definition and it is accurate at many data points that have $\tau_{i}=\tau_{\min}$. To clarify this point, we compared the likelihood when the data points with $\tau_{i}=\tau_{\min}$ were excluded. The results are shown in Table 3. To simplify the discussion, for the EMM we employed the value of $k$ that minimized the AIC. For the Reality data, the EMM fits the data better than the Pareto distribution when we only consider the inter-event times larger than $\tau_{\min}$. For the Hypertext data set, the Pareto distribution is still better than the EMM, but the difference in the likelihood is much smaller than when all the inter-event times are used. Note that the Reality data set has a larger fraction of inter-event times with $\tau=\tau_{\min}$ than the Hypertext data set (Table 3).

The PLFit only fits a power-law distribution to relatively large inter-event times. Therefore, we compared the likelihood of the EMM, Pareto distribution, and PLFit for the set of inter-event time values satisfying $\tau_{i}\geq\hat{b}$, where $\hat{b}$ is the threshold estimated by the PLFit. The PLFit is the best performer among the three fitting schemes (Table 3). This is probably because the threshold is optimized for the PLFit in this comparison and the PLFit is allowed to discard small inter-event times. In contrast, the EMM intends to fit the entire set of inter-event times. Furthermore, in this comparison, the EMM is disadvantageous by the amount of probability mass that it has to devote to inter-event times smaller than $\hat{b}$. It should be noted that the PLFit discards a considerable fraction of data points, ranging between 27–94%, except for the College data set (Table 3). It should also be noted that the EMM yields a larger likelihood value than the Pareto distribution for all the data sets with this $\hat{b}$ value.

We do not conclude whether or not the EMM outperforms the Pareto distribution or PLFit. Both the Pareto distribution and PLFit benefit from carefully setting the lower bound on the inter-event times, $\hat{b}$. When all inter-event time data have to be modeled and we do not know the smallest possible value of the inter-event time, EMMs may be preferred to the Pareto distribution and the PLFit.

So far, we have excluded the inter-event times across consecutive days in the Office and Hypertext data sets. We also ran the same analysis for these two data sets without excluding such long inter-event times. The results were qualitatively the same except that the EMM estimated the number of components larger by one (Supplementary Materials). The new component (i.e. a single exponential distribution in the EMM) corresponded to the across-day inter-event times, with a considerably longer mean inter-event time than that for the other estimated components. Furthermore, the mean inter-event time of the constituent exponential distributions that exist both when across-day inter-event times are included and when excluded are almost the same in the two cases (Supplementary Materials). We conclude that inter-event times on the timescale of a day or longer contribute one exponential distribution with a large mean to the entire EMM without interfering with the constituent exponential distributions on smaller time scales.

3.3 Population results

In the previous section, we showed that mixtures of a small number of exponential distributions tended to be selected according to the different model selection criteria. To assess the generality of this result, we carried out the model selection for all individuals having at least 100 inter-event times. Because the Sexual data set has only two individuals with at least 100 inter-event times, we instead used the 56 individuals that had at least 50 inter-event times for this data set. The histogram of the optimal $k^{*}$ value for each combination of the data set and criterion is shown in Fig. 2. Consistently with the results for the individuals with the largest number of inter-event times in each data set, the effective number of components, $k^{*}$, was estimated to be at most three in four out of the seven data sets (i.e., Office, Hypertext, Reality, Sexual). For the Bitcoin data set, the distributions of $k^{*}$ had the mode at 3, and we obtained $2\leq k^{*}\leq 4$ in most cases. For the other two data sets (i.e., Email and College), the AIC and BIC yielded distributions of $k^{*}$ whose mode was four and the largest value was seven. Under the other four criteria, which are the ones justified for EMMs, the mode was equal to three, and the largest $k^{*}$ was equal to five.

To exclude the possibility that longer data (i.e., larger $n$) tend to yield a larger $k^{*}$ value, we calculated the Pearson correlation coefficient between $n$ and the optimal $k^{*}$ for each combination of the data set and criterion. To calculate the Pearson correlation coefficient, we only used the individuals with at least 100 inter-event times (and at least 50 inter-event times for the Sexual data set), similar to Fig. 2. The Pearson correlation coefficient values, denoted by $r$, and its $p$ value are shown in Table 4. The correlation coefficient was small and insignificant except for the Email data set combined with all criteria and the College data set combined with the AIC or BIC. When we confined the analysis to the individuals with larger numbers of inter-event times, i.e., at least 500 and 200 inter-event times for the Email and College data sets, respectively, the Pearson correlation coefficient was smaller. In particular, for the $\text{AIC}_{\text{LVC}}$, $\text{BIC}_{\text{LVC}}$, $\text{NML}_{\text{LVC}}$, and DNML, the correlation was close to zero and insignificant when we only considered the individuals with many inter-event times. Note that the largest number of inter-event times for these two data sets were sufficiently larger than the new threshold values (Email: $n=3948$, College: $n=1090$) and that there remained sufficiently many individuals with the new threshold values (Email: 117 individuals, College: 61 individuals). The results suggest that our main result that the optimal number of components, $k^{*}$, tends to be small is expected to remain true for long sequences of inter-event times.

4 Discussion

We showed that EMMs (i.e., mixtures of exponential distributions) with a small number of components, up to three in many cases and four in some cases, were a reasonably good fit to various empirical distributions of inter-event times. In general, there are various mechanisms behind power-law distributions of an observable, and many of them are believed to be relevant as generative mechanisms of long-tailed degree distributions of empirical networks [46]. In contrast, the present results suggest that empirical long-tailed distributions can be also approximated by EMMs. Furthermore, EMMs can be obtained as an outcome of stochastic point processes in which the system’s state switches between a small number of discrete states in each of which the process generates events as a Poisson process. In practice, a human or animal individual may maintain a small number of relatively discrete states, such as active and rest. The event rate may depend on the discrete state. It should be noted that, within each state, the mechanism to generate event sequences is maximally simple: Poisson process with a constant rate. To the best of our knowledge, Poissonian views of long-tailed distributions of inter-event times were first introduced in [5]. The contribution of the present work is to have introduced a formal model selection framework to directly compare EMMs with different numbers of components and having carried out some comparisons with power-law distributions.

There is a recent debate regarding the abundance of scale-free (i.e., power-law) degree distributions in empirical networks. When a pure power-law tail is assumed for the degree distribution, there are empirically few scale-free networks [47]. In contrast, scale-free degree distributions are abundant when impurity in the power-law distributions via the regularly varying distributions is considered [48]. In the present study, we fitted the pure power-law distributions via the PLFit algorithm [36]. The PLFit was more accurate than EMMs while the PLFit discarded small inter-event times, which occupied a non-negligible fraction of inter-event times in each data set. We did not examine the method presented in Ref. [48]. This is because the purpose of the present study was not to find an accurate estimator or to test if a power-law distribution of inter-event times was a good fit. Rather, our purpose was to test the hypothesis that the system generating time-stamped events, such as humans, transit between a small number of discrete states each of which is a Poisson process generator.

A reasonably accurate fit of the EMMs to the empirical data revealed in the present study mainly comes from a high accuracy at small $\tau$ values because there are many data points with small $\tau$. Such short inter-event times yield a component exponential distribution with a small mean inter-event time and may be produced as a result of the visit of the individual to an active state. Given these considerations, looking at the entire distribution of the data rather than focusing on the distribution’s tail may help us to understand mechanisms generating the data.

The authors of [10] reached a similar conclusion to ours, where they split the inter-event times into two groups by thresholding. Then, they fitted an exponential distribution to the inter-event times exceeding the threshold and another exponential distribution to the inter-event times less than the threshold. Their models are EMMs with two components. In contrast to their study, we systematically fitted EMMs and inferred the number of components using model selection criteria valid for EMMs. Testing our method on various other distributions of inter-event times to clarity the reason why EMMs with a small number of components apply accurately to some data sets and not others warrants future work.

Data accessibility

The data are open resources. Python codes for estimating and selecting the mixture of exponential distributions are available at Github (https://github.com/naokimas/exp_mixture_model).

Authors’ contributions

KY and NM conceived the study. NM designed the study. MO carried out the theoretical analysis and performed the data analysis. MO, KY and NM discussed the results and drafted the manuscript. All authors read and approved the manuscript.

Competing interests

The authors declare that they have no competing interests.

Funding

K.Y. and N.M. acknowledge the support provided through JST CREST Grant Number JPMJCR1304, Japan.

4.1 Acknowledgments

We thank the SocioPatterns collaboration (http://www.sociopatterns.org) for providing the two data sets used in the present paper.