Testing statistical laws in complex systems

TL;DR

This paper examines how correlations in complex system data affect the testing of statistical laws, proposing a conservative method to improve the reliability of such tests.

Contribution

It introduces a new testing approach that accounts for correlations, reducing false rejections of statistical laws in complex systems analysis.

Findings

Correlations cause standard tests to falsely reject statistical laws.

The proposed method reduces false rejections and increases confidence in law validation.

Accounting for correlations leads to more accurate statistical law assessments.

Abstract

The availability of large datasets requires an improved view on statistical laws in complex systems, such as Zipf's law of word frequencies, the Gutenberg-Richter law of earthquake magnitudes, or scale-free degree distribution in networks. In this paper we discuss how the statistical analysis of these laws are affected by correlations present in the observations, the typical scenario for data from complex systems. We first show how standard maximum-likelihood recipes lead to false rejections of statistical laws in the presence of correlations. We then propose a conservative method (based on shuffling and under-sampling the data) to test statistical laws and find that accounting for correlations leads to smaller rejection rates and larger confidence intervals on estimated parameters.

Figures

d of the manuscript. Correlations in the sequence $x_{i}$ exist even forbidding local steps because when we sample edges we sample two nodes and the selection of the two is not random, e.g., degree-degree correlation in networks. Another distinction of our approach is that we work in the rank-frequency picture so that node $x$ with degree $k_{x}$ contributes with $k_{x}$ symbols $x$ in the sequence $x_{i}$ (in contrast to the analysis of the distribution $P(k)$ of nodes of degree $k$).

Equations

$$p(x;\alpha )=C{x}^{-\alpha }.$$

$$\log \mathcal{L}=\sum _{i=1}^N\log p(x_i;\vec{\alpha }).$$

$$p(x;\vec{\alpha })=p(x,\alpha )=\alpha e^{-\alpha (x-x_{\min })},$$

$$F(x;\vec{\alpha })\equiv P(X\ge x;\vec{\alpha })=e^{-\alpha (x^{\beta }-x_{\min }^{\beta })}.$$

$$p(x;\alpha )=C{x}^{-\alpha }$$

$$p(x_{i+1}|x_i)=\{\begin{array}{ll}p(x;\alpha )& \text{ w/ probab. }\mu \\ \frac{1}{|B(x_i)|}\delta (x_i\in B(x_i))& \text{ w/ prob. }1-\mu ,\end{array}$$

$$A(x\mapsto x^{\prime })=\frac{g(x^{\prime }\mapsto x)p(x^{\prime })}{g(x\mapsto x^{\prime })p(x)},$$

$$g(x\mapsto x^{\prime })=\mu p(x^{\prime })+(1-\mu )/|B(x)|.$$

$${\tilde{\tau }}^{\ast }=\sum _{\tau =-\infty }^{\infty }(C(\tau ){)}^2,$$

Full text

Introduction

Statistical regularities collected in the form of “universal laws” play a central role in complex systems Newman ; Mitzenmacher2004 ; Laws . Zipf’s law of word frequencies Zipf , the Gutenberg-Richter law of earthquake magnitudes Sornette , scale-free degree distributions in networks BarabasiAlbert1999 , and inter-event time distributions between bursty events Bunde ; Barabasi ; bursts ; burstsC are prominent examples that triggered entire research lines devoted to explaining the origin and to exploring the consequences of these laws.

Recently, the empirical support of such laws has been heavily questioned. The best known example is the case of scale-free degree distribution of networks: after the seminal work of Barabasi and Albert in 1999 BarabasiAlbert1999 , the early 2000’s were marked by findings of power-law distributions in various network datasets, while in the last five years the trend has reversed and it is now common to read that networks with power-law degree distribution are rare Khanin ; Broido (see Ref. Quanta2018 for a journalistic account). This recent shift in conclusions, which appears in the analysis of Zipf’s law in language PRX ; Laws ; Francesc and also in other areas CriticalTruth ; CSN , is partially due to new (larger) datasets but mostly due to the improved statistical methods: least-squared fitting and visual inspection of double-logarithmic plots (used since Zipf) have been replaced by maximum likelihood methods made popular in the influential article by Clauset, Shalizi, and Newman CSN , see Refs. Goldstein04 ; Bauke07 ; Deluca13 ; Hanel for variations. A point often ignored in the interpretations of the recent findings is that these methods rely on two hypotheses:

• H1:

The observations $x$ are distributed as $p(x;\vec{\alpha})$, where $\vec{\alpha}$ are parameters, e.g. for a power law

[TABLE]

• H2:

The empirical observations $x_{i},i=1\ldots N$ are independent (e.g., of $i$ or $x_{i-1}$).

While the statistical laws correspond to H1, the statistical tests rely also on H2 (implicitly assumed, e.g., when the log-likelihood is computed as $\sum_{i=1}^{N}\log p(x_{i})$ StumpfEPL ; Khanin ; CSN ; Hanel ; Sornette ). Complex systems are characterized by strong (temporal and spatial) inter-dependencies Eisler.2008 and it is thus not clear whether the recent claims Khanin ; CriticalTruth ; Broido of violation of the statistical laws arise from systematic deviations of the law itself (H1) or, instead, whether they are due to the well-known fact that observations are not independent (H2).

In this Letter we show that dependencies in the data (violation of H2) have a strong impact on the empirical analysis of statistical laws, leading to rejections even in processes that satisfy the law (H1), and to over-confident selection of models and parameters. We then propose an alternative method that distinguishes between H1 and H2, yielding an upper bound on the degree of correlations for which the statistical law is rejected.

General setting.

Let $\{x_{i}\}=x_{1},x_{2},\ldots,x_{N}$ be an ordered sequence obtained from a measurement process that asymptotically has a well defined distribution $p(x)=\frac{\#x_{i}=x}{N}$ as $N\rightarrow\infty$. In observations of dynamical systems (or time series), $x_{i}$ will typically depend on the observations at previous times so that for all times $\tau$ smaller than some (relaxation) time $\tau^{*}$ we find $p(x_{i}|x_{i-\tau})\neq p(x)$. Violations of H2 happen also when data is not measured as a time-series. In the case of Zipf’s law of word frequencies, syntax restrict the valid sequences of word tokens, in violation of H2 (both in the rank-frequency and frequency distribution pictures Laws ; Francesc ). In the case of networks, H2 can be violated because of the generative process or because of the sampling employed to observe the nodes and links (typically a subsample of an underlying network). In fact, it has been shown that the degree distribution of networks is sensitive to the sampling procedure StumpfPNAS ; StumpfPRE ; Lee . Moreover, the hypotheses H1 and H2 of the standard tests for power-law distribution do not build a proper probabilistic network model Note1 , are thus not suitable to a rigorous statistical analysis Crane , and the analysis of the degree distribution of networks requires further assumptions about the sampling/generative process.

More generally, strong correlations are ubiquitous in complex systems Eisler.2008 and it is hard to imagine a case for which H2 holds. In Fig. 1 we show how previously proposed statistical laws and correlations appear together in paradigmatic complex systems: the Gutenberg-Richter law for earthquakes (exponential Sornette ), interevent times of words (stretched exponential Bunde ; bursts ; burstsC ), Zipf’s law for word frequencies (power-law Zipf ), and scale-free distribution for the node-degree in networks (power-law BarabasiAlbert1999 ). While earthquake events and interevent times naturally occur as time series data, we mapped word frequencies in texts and the network data into ordered sequences $\{x_{i}\}$ based on a simple sampling process (see caption of Fig. 1) in order to illustrate and quantify the violation of H2 in an unified framework.

Constructed example.

We now show that the traditional methods CSN lead to a rejection of a power-law distribution (1) even for data which are power-law distributed for $N\rightarrow\infty$. This is done by building a Markov process B1 ; B2 in which H1 is satisfied but H2 is violated (i.e., $x_{i}$ depends on $x_{i-1}$ and $p(x)=Cx^{\alpha}$ for $N\rightarrow\infty$, see SM Sec. III).

In Fig. 2 we show that the violations of H2 have a strong influence on the analysis of statistical laws formulated in H1. In particular, the application of the traditional recipes CSN lead to the wrong conclusion that the data is not compatible with a power law distribution: the probability of rejecting the null hypothesis at a $5\%$ significance level is much larger than $5\%$ even for small sample sizes $N$ (inset of Fig. 2b). This corresponds to a type-I error because, by construction, the data satisfies H1. The origin of this failure thus originates from the fact that correlations lead to an effective reduction of the number of independent observations implying larger fluctuations which lead to larger deviations from the fitted model. Specifically, we recall that the test employed in Ref. CSN consists of comparing the Kolmogorov-Smirnov (KS) distance between the correlated data and the fitted curve, $KS_{correlated}$ (blue curve), and the KS distance between independent samples of the model (H1+H2) and the fitted curve, $KS_{model}$ (orange curve). More precisely, the statistical law is rejected at $5\%$ significance level if $KS_{correlated}>KS_{model}$ in $95\%$ realizations (samplings) of the model. While in our artificial data $KS_{correlated}\propto 1/\sqrt{N}$ (as expected) and thus $KS_{correlated}\rightarrow 0$ for $N\rightarrow\infty$, this convergence is shifted from the convergence of $KS_{model}$ (Fig. 2b) due to the correlations. This shift leads to an increased rejection rate ($\approx 1$, p-value $\approx 0$).

Violations of H2 are important not only in the hypothesis-testing setting discussed above, they also lead to increased systematic and statistical errors (bias and fluctuations) in the fitting of the parameter $\hat{\alpha}$ (Fig. 2c) and, thus, in the selection between models Hastie09 ; Burnham02 .

Real data.

In order to confirm that the results discussed above are also relevant in real datasets – which have a fixed size $N$ – we consider two types of undersampling of data to sizes $n<N$: taking $n$ points either randomly or preserving the structures/correlations by taking consecutive portions of the time series (the network and word-frequency databases are first mapped to a time series, as in Fig. 1). In order to distinguish between the effect of the shape of the distribution (H1) and correlations (H2) we compare the distribution of the $n$ points with i) the proposed statistical law and ii) the empirical distribution (i.e., the one obtained for $n=N$). Our results (see SM, Sec. IV) confirm that correlated data show higher rejection-rate and fluctuations of parameters.

Alternative approach.

In the vast literature of statistical methods for dependent data, two general approaches can be identified. The first approach is to incorporate the violation of independence in more sophisticated (parametric) models, e.g., in time series one could consider Gaussian/Markov processes vankampen . This is of limited use in our case because statistical laws aim to provide a coarse-grained description (stylized facts) valid in many systems, instead of different detailed models of particular cases. The second (non-parametric) approach, which we pursue here, is to de-correlate or de-cluster the data, leading to a dataset with an “effective” sample size $N^{*}\leq N$ Gasser ; Weiss78 ; Chicheportiche11 . In practice, the analysis consists of multiple realizations of the following three steps:

• (i)

Randomize (shuffle) the original sequence and select randomly $n$ points, for different $n\in[1,N]$.

• (ii)

Apply the traditional statistical analysis (i.e., the hypothesis test, model comparison, and fitting based on H1+H2) to the randomized dataset obtained in (i), investigating their dependence on $n$.

• (iii)

Estimate the correlation $\tau^{*}$, defined as the time after which two observations (in the time series) are independent from each other. Out of the total $N$ samples we thus estimate $N^{*}=N/\tau^{*}$ to be the number of independent samples and therefore we select the results from step (ii) for $n\approx N^{*}$.

The determination of $\tau^{*}$ – or the effective sample size $N^{*}$ – in step (iii) requires knowledge or assumptions about how the data was generated. For the case of temporal sequences we propose to compute the auto-correlation and take as $\tau^{*}$ the lag for which it reaches an interval around zero (1-percentile of the random realizations, as in Fig. 1). In the constructed example (Fig. 2), we obtain $\tau^{*}=407\Rightarrow N^{*}=245$ which leads to a rejection-rate (at p-value$=0.05$) equal to $5\%$ for all $n<N^{*}$. For the case of networks, the determination of the effective sample size $N^{*}$ depends on the generative process and/or the sampling used to measure the data (here we assumed a specific edge sampling method, as described in Fig. 1.) In Fig. 3 we show evidence of the effectiveness of our approach through a systematic analysis of the p-value distribution as a function of $n$ for both the constructed and empirical datasets. This is further corroborated in artificial data (see SM, Sec. V) showing that i) our method for the selection of $\tau^{*}$ is superior to the one proposed in Ref. Gasser (sum of the autocorrelation function) and ii) can be equally applied to data with other types of correlation: a Markov process with negative correlation and a Gaussian process with long-range correlations. In all cases our approach shows an uniform distribution of p-values under the null hypothesis.

An important message of our analysis is that conclusions about the statistical law can be obtained even when the precise value of $\tau^{*}$ (or the effective sample size $N^{*}$) is unknown in step (iii). By shuffling and undersampling the sequence at different sizes $n$ – steps (i) and (ii) – we can investigate how the results depend on $n$ and obtain the range in $\tau^{*}$ for which the different conclusions hold. For instance, in the case of earthquakes – Fig. 3c – we see that the rejection increases dramatically around $n\approx 10^{3}$. We thus conclude that, in this dataset of size $N\approx 10^{5}$, we falsify the Gutenberg-Richter law if $\tau^{*}\leq N/n\approx 10^{2}\text{ observations }\approx 20$ days Note3 . The conservative estimate of $\tau^{*}$ in Fig. 1 was $\tau^{*}=791>10^{2}$ and therefore we conclude that based on this data we cannot reject the Gutenberg-Richter law, contrary to the conclusion obtained assuming independent observations. We find similar results – Fig. 3d – for the stretched exponential distribution of inter-event times between words, while for Zipf’s law – Fig. 3c – the outcome is uncertain, and the power-law degree distributions in networks – Fig. 3d – is rejected even in the correlated case.

Discussion and Conclusion

Statistical laws in complex systems are typically formulated (as in H1) without reference to the generative process of the data. Therefore, ideally, the empirical test of these laws should be designed to account for a large class of processes generating $\{x_{i}\}$. Traditional methods CSN based on the hypothesis of independent data (H2) are weak tests because they include a strong hypothesis that is easily violated, therefore favoring rejection. In fact, here we have shown how these methods: (i) lead to wrong rejections of the laws because of correlated data; and (ii) are over-optimistic regarding uncertainties of the estimated parameters. Stronger tests of statistical laws should make weaker assumptions about the generative process so that rejections of the compound hypothesis provide much stronger evidence of the rejection of the law (H1). Here we proposed a methodology which allows us to identify the strongest assumption about correlations of the data $\tau^{*}$ for which the law can be rejected. Being conservative in the choice of $\tau^{*}$ (i.e., choosing large values for which we are confident that $x_{i}$ and $x_{i+\tau^{*}}$ are uncorrelated) overcomes the main shortcoming of the traditional approach CSN and ensures that when we reject the law this is not happening due to correlations in the data (failing to reject the law is never a confirmation of its validity). In this sense, our approach is similar in spirit to the Bonferroni correction to account for multiple hypothesis testing bonferroni (both aim to avoid overconfident or spurious rejections of hypotheses). Our approach is even applicable in cases with no well-defined mixing time $\tau^{*}$ (e.g. long-range correlations) because it yields very large values of $\tau^{*}\approx N$ (no two points are independent).

Instead of directly testing whether the statistical law is valid (hypothesis testing), often the best we can do is to compare different alternatives (model comparison) StumpfEPL ; Broido ; CSN ; Laws ; PRX ; Hastie09 ; Burnham02 . Also in this case violations of the hypothesis of independence are important and have been mostly ignored in the analysis of statistical laws in complex systems (see Refs. Lee ; Khanin ; Sornette for exceptions). As shown above, due to correlations (and violations of H2) actual data show much larger fluctuations than expected under the hypothesis of independent observations. By using a shuffled and undersampled dataset we obtain larger uncertainties in the estimated parameters; we expect similar lack of certainty in the choice of best models. The need to account for violations of the independence assumption, shown in this Letter, applies much more broadly than the cases treated above. Correlations should be accounted for whenever testing statistical laws in complex systems, such as linguistic laws Laws , scaling laws with system size – maximum likelihood methods based on H2 have been applied to biological allometric laws dodds and to city data us – and different distributions of inter-event time (burstiness) Bunde ; Barabasi ; bursts ; burstsC .

Acknowledgments:

We thank F. Font-Clos and J. Moore for the careful reading of the manuscript.

Code and data availability:

see Ref. codes .

I Fitting

Given a set of data $\left\{x_{i}\right\}=x_{1},x_{2},\ldots,x_{N}$ we find the best parameters $\hat{\vec{\alpha}}$ by maximizing the log-likelihood:

[TABLE]

Exponential.

For the earthquake data, the observable $x$ is the magnitude. We thus fit a continuous exponential distribution with support $x\in[x_{\min},\infty)$

[TABLE]

where $x_{\min}$ is a constant. We choose $x_{\min}=2$ according to Ref. [29] main text.

Stretched Exponential.

For the interevent times data, the observable $x$ is the number of words between two consecutive occurrences of a given word (e.g. “the”) in a book. We fit a discrete stretched exponential with two parameters $\vec{\alpha}=(\alpha,\beta)$ and support $x=x_{\min},x_{\min}+1,\ldots,\infty$ defined by its cumulative distribution

[TABLE]

From this we obtain the probability distribution as $P(x;\vec{\alpha})=F(x;\vec{\alpha})-F(x+1;\vec{\alpha})$. We choose $x_{\min}=3$ as the minimum value for which estimation of the parameters is approximately independent of $x_{\min}$ (not shown).

Power law.

For the text and network data, the observable $x$ is a rank. We thus fit a discrete power-law with support $x=1,2,\ldots,V$ with $V$ the maximum (observed) rank and $\alpha>0$ such that

[TABLE]

with $C=\frac{1}{H_{V,\alpha}}$, where $H_{V,\alpha}\equiv\sum_{x=1}^{V}x^{-\alpha}$ is the generalized harmonic number of order $\alpha$ of $V$.

II Edge-sampling networks

Here we describe the procedure used to obtain an ordered sequence of nodes $x_{i},i=1,\ldots,2L$ starting from a simple graph of $V$ nodes and $L=N/2$ edges (or links). We select nodes by choosing the edges $x_{i}\leftrightarrow x_{j}$ of a network in a particular order (without replacement) and then adding the two nodes linked by this edge ($x_{i},x_{j}$) in a random order to the list of observations $x_{i}$. A node with degree $k$ thus appears $k$ times in the sequence. If the node labels correspond to their ranks (in degree), $p(x)$ corresponds to the rank frequency distribution.

We start from a randomly selected edge. The next edge is selected randomly from the remaining list of edges involving $x_{i}$ or $x_{j}$ (with probability $0.5$) or randomly from the complete list of remaining edges in the network (with probability $0.5$). If there are no remaining edges involving $x_{i}$ or $x_{j}$ we choose an edge randomly in the complete network. This corresponds to a combination of random sampling of edges with a random walk (using local steps) in the network. This procedure is repeated until all edges $L$ of the network are sampled, leading thus to a sequence $x_{i}$ for $i=1,\ldots,2L$.

Lee et al. [SM4] showed that random edge sampling is equivalent to node sampling (previously discussed in Ref. [SM5]), both leading to deviations of the degree distribution from power-law. The local steps used here enhance the extent into which the measured $x_{i}$ deviate from an independent sample (H2), as shown in Fig. 1d of the manuscript. Correlations in the sequence $x_{i}$ exist even forbidding local steps because when we sample edges we sample two nodes and the selection of the two is not random, e.g., degree-degree correlation in networks. Another distinction of our approach is that we work in the rank-frequency picture so that node $x$ with degree $k_{x}$ contributes with $k_{x}$ symbols $x$ in the sequence $x_{i}$ (in contrast to the analysis of the distribution $P(k)$ of nodes of degree $k$).

III Artificial Data

Here we describe how to generate a correlated sequence $x_{1},x_{2},\ldots,x_{i},\ldots,x_{N}$ with variable correlation time $\tau^{*}$ and an arbitrary marginal distribution $p(x)$ for $N\rightarrow\infty$. We construct a Markov process of order one $p(x_{i+1}|x_{i},x_{i-1},\ldots)=p(x_{i+1}|x_{i})$ such that the next value $x_{i+1}$ is proposed from $p(x)$ with probability $\mu$ and from a pre-defined set of neighbours $B(x_{i})$ with probability $1-\mu$:

[TABLE]

where $|B(x_{i})|$ is the number of elements of $B(x_{i})$. The strict validity of $p(x)$ is ensured by accepting or rejecting this proposal using the Metropolis-Hasting method [SM1]. Starting from a random value $x$, a new value $x^{\prime}$ is proposed according to Eq. (5). This proposal is accepted with a probability given by the Metropolis-Hasting condition

[TABLE]

where $g(x\mapsto x^{\prime})$ is the probability of proposing $x^{\prime}$ from $x$, and $p(x)$ is the probability distribution we wish to impose for $x$ (e.g., the power-law distribution). If $A>1$, $x^{\prime}$ is accepted. If $x^{\prime}$ is accepted, $x_{i+1}=x^{\prime}$, otherwise $x_{i+1}=x_{i}$. Each $x$ has neighbours $B(x)$ that we assume to be reciprocal, i.e. $x^{\prime}\in B(x)\Rightarrow x\in B(x^{\prime})$ for all $x,x^{\prime}$. This implies that $A=1$ if $x$ and $x^{\prime}$ are not neighbours of each other. If they are neighbours, the probabilities $g$ obtained from the rule (5) are

[TABLE]

The procedure above is repeated $N$ times. The autocorrelation of $x$ decays exponentially and $\frac{\#x}{N}\rightarrow p(x)$ for $N\rightarrow\infty$. For $\mu=1$ the data becomes independently distributed in agreement with H2. For $\mu<1$ the sequence of $x_{i}$ is correlated with a correlation time $\tau\lessapprox 1/\mu$. In the example shown in the manuscript we consider: $p(x)$ to be a power-law distribution – Eq. (1) of the paper – with support $x=1,\ldots,1000$ (with qualitatively same results if other maximum cutoffs are chosen, see Fig. 1), $\mu=0.01$, and the neighbours $B(x_{i})$ of $x_{i}$ in the local step $x_{i+1}$ to be the $k-$nearest neighbours $x_{i}-k,\ldots,x_{i}-1,x_{i},x_{i}+1,\ldots x_{i}+k$ (i.e., we choose one of the $k$-nearest neighbours by chance with equal probability) setting $k=5$.

The process described above is a generic process that aims to capture some of the features of different systems. In a text, the neighbours of $x_{i}$ could correspond to the word types which are syntactically allowed to follow $x_{i}$. In a network for which nodes (and their degrees $x_{i}$) are discovered through a random walk through its links, the step of choosing neighbours of $x_{i}$ mimics the cases in which nodes connected to a given node have a degree more similar to it than a random node (positive degree assortativity of the underlying network) leading to $p(x_{i+1}\approx x_{i}|x_{i})>p(x_{i};\alpha)$.

IV Separating correlations and models in fitting real data

Here we show how correlations in real data affect the fit of different models similar to the analysis on synthetic data. In contrast to synthetic data, for real data we do not know the generating process and thus we don’t know the underlying distribution and we are unable to generate datasets of arbitrary length $n<N$. Therefore, when fitting a model to the full data, it is impossible to know whether rejection of the hypothesis is due to correlations or due to deviations from the underlying distribution. In order to disentangle the two effects, isolating the effect of correlations, we consider two approaches. First, we shuffle and undersample the data randomly (as argued in the main manuscript). Second, we compare the data not only to the parameterized statistical law but also to a 0-parameter function defined by the empirical distribution of the full data set.

The results for the four datasets used in the main manuscript are shown in SM-Figs. 2, 3, 4, and 5. We calculate the KS distance between the empirical distribution and correlated and uncorrelated datasets, respectively, of different size obtained from subsampling (panel a). We observe that the KS-distance for the correlated data is much larger than for the uncorrelated data. While the uncorrelated data is by construction rejected at a rate of $0.05$, the correlated data is rejected with a much higher probability. This confirms that violations of the hypothesis of independent data (H2) lead to rejections. Fitting the parametric statistical law leads to qualitatively similar results (panel b). The correlated data yields much larger values of KS-distance to the fitted distribution than the uncorrelated data. This in turn leads to an increase in the rejection rate of the model. Beyond hypothesis testing, we also observe that the correlations have a substantial impact on the estimation of the parameters of the respective model (panel c) – not only in terms of an increase in the fluctuations but also in the average value.

V Testing our method in artificial data

Here we report numerical tests of the method we propose in the main manuscript applied to artificially-generated data. We compare the outcome of our method to the method proposed by Gasser [SM2]

The artificial data (described in Sec. III above) corresponds to a Markov process with stationary distribution $p(x;\alpha)$. Starting from an arbitrary initial condition $x_{i}$, the probability that after $t$ times the state $x_{i+t}$ equals $x$ approaches the stationary distribution $p(x;\alpha)$ for long times $t$, with the typical time for this to occur given by the mixing time $\tau_{m}$ of the Markov Chain [SM3]. Points separated by times $t>\tau_{m}$ can thus be considered as independent samples from $p(x;\alpha)$. The method we propose selects randomly $n$ out of the $N$ points in the data, implying that these points are separated by a distance $N/n$. In our method, the specific value of $n$ chose as the effective sample size $n=N^{*}$ is determined by the correlation time $\tau^{*}$ – computed from the autocorrelation function of $\{x_{t}\}$ – so that $N/N^{*}=\tau^{*}$. Our goal of having an effective sample size of independent points is achieved if $\tau^{*}>\tau_{m}$, while ideally $\tau^{*}\approx\tau_{m}$.

The reasoning above indicates that the identification of the time $\tau^{*}$ is a critical element of our method. The alternative method we consider here, by Gasser [SM2], proposes an alternative estimation of $\tau^{*}$ as

[TABLE]

where $C(\tau)$ is the auto-correlation function of the series $x_{i}$ at lag $\tau$.

We first consider the artificial data described in the previous section for different values of the free-parameter $\mu$. SM-Figure 6 shows the comparison between the traditional approach (ignoring correlations), our method, and Gasser’s method. The superiority of our method is confirmed by the fact that we obtain a flat distribution of p-values (panel a-c), expected from the fact that by construction the stationary distribution of the Markov Chain equals the fitting distribution. In order to quantify the deviation from a flat distribution, we calculate the entropy over the distribution of p-values varying the correlation strengths $\mu$ (panel d). Our method consistently yields an entropy virtually indistinguishable from the maximum value expected from a flat distribution (dotted line). In contrast, the other methods yield substantially smaller entropy values, reflecting a deviation from a flat distribution. This in turn translate into much smaller or larger fraction of realizations being rejected than expected from the imposed significance threshold, i.e. p-value $<0.05$ (panel e). The main reason for the best performance of our method is that the estimated correlation time $\tau^{*}$ is larger than the one estimated by the alternative method in Eq. (6) (see panel f).

Our method applies also to datasets showing different types of auto-correlation functions. To confirm this, we performed comparisons – similar to the ones reported above – using other types of artificial data that have a known marginal distribution $p(x;\alpha)$ but different types of autocorrelation. For instance, in SM-Fig. 7 we test power-law data with negative correlations (for odd time lags) and in SM-Fig. 8 we test Gaussian data with long-range correlations. Even though both samples follow by construction the tested statistical law, the distribution of p-values is peaked at very small values, leading to a much larger rate of rejection than expected by chance. Applying our methodology (i.e., estimating a correlation time $\tau^{*}$ from the auto-correlation, shuffling, and subsampling the data to a smaller size $n=N^{*}=N/\tau^{*}$) yields an approximately uniform distribution of p-values consistent with the validity of the null hypothesis.