Seeing through noise in power laws

2.1. Methods overview

Here, we provide a minimal description of the methods and notation, reserving a detailed description for the electronic supplementary material, S1, while we discuss analysis and validation in §3. We estimate α using MLEs and test goodness of fit using the K-S statistic, in accordance with widely accepted best practices [15]. We bin the input data by truncation, replacing each data point x_i with the binned value $\lfloor x_i{\rfloor }_{\lambda }=x_m{\lambda }^{\lfloor {\log }_{\lambda }{x}_i/x_m\rfloor }$, so that the data takes on the discrete set of values x_mλ^k for k = 0, 1, 2, …. The likelihood for power-law data binned this way is also a power law [26],

Binned data nonetheless entail different MLEs and null distributions of the K-S statistic [34]. We therefore notate the MLEs ${\hat{\alpha }}_{\lambda }$ to emphasize that the MLE depends on the binning ratio λ and we use an appropriate bootstrapped null distribution of K-S statistics [14,35] (electronic supplementary material, S1). The limit λ → 1⁺ corresponds to no binning, and so we use the symbol ${\hat{\alpha }}_1$ for the usual continuous MLE [15] and λ = 1 to denote the raw data. Naively applying ${\hat{\alpha }}_1$ to binned data results in substantially biased estimates even for very fine binning with λ close to 1 [26,34].

2.2. Earthquakes

Earthquake magnitudes are historically believed to follow a power law, where the exponent is the Gutenberg–Richter b-value [5,19,31]. The earthquake dataset [15] contains 17 450 samples ranging over seven orders of magnitude from 0.5 to 7.8 on the Richter scale. Traditional Richter magnitudes inherently constitute a case of logarithmic binning, because they record the base-10 logarithm of wave amplitudes to two digits of precision [36]. Exponentiating these magnitudes onto their natural scale results in a dataset with discrete values proportional to integer powers of λ = 10^0.1 ≈ 1.26 ranging over seven orders of magnitude. Partial censorship is known to complicate inference below the magnitude of completeness m_c [37], which depends on the density of the earthquake sensor network. We choose x_m = 10^3.5, or magnitude 3.5, to visibly exceed m_c, retaining 5910 data points.

The results in figure 1a and table 1 substantially depend on binning. We conduct inference with no binning (λ = 1), λ = 10^0.1 and λ = 10. Choosing λ = 10^0.1 matches the quantization of the Richter scale. This ‘binning’ preserves the input data exactly but entails a different MLE and null distribution of K-S statistics. The raw and fine-binned empirical distributions therefore overlap in figure 1a but imply different interpretations of the data. The coarse binning λ = 10 represents increments of 1.0 in magnitude and divides the data into five bins.

Differences in the results point to features of the data rather than the underlying mechanisms. The effect of quantization noise on MLEs is evident from the significant difference between the continuous MLE ${\hat{\alpha }}_1=0.837\pm 0.006$, which falsely assumes no quantization, and the MLE that accounts for the quantization ${\hat{\alpha }}_{{10}^{0.1}}=0.77\pm 0.01$. Quantization noise therefore biases ${\hat{\alpha }}_1$ by approximately 0.07 or 9%, comparable to the 12% bias observed in other earthquake data [26]. Because errors in α occur in the exponent, however, this 9% correction predicts a nearly twofold higher rate of earthquakes larger than 7.8, the largest that appear in the dataset.

Such quantization noise would be sufficient reason for K-S tests to reject the pure power law, yet the test also rejects the power law when accounting for quantization with λ = 10^0.1 (p < 0.01). This is because special magnitudes 3.0, 3.5, 4.0, 4.5, …, contain an excess of observations, as can be observed in the histogram (electronic supplementary material, figure S1) and quantile–quantile (QQ) plot (electronic supplementary material, figure S2). This excess creates a statistically significant association between ‘round’ numbers and higher-than-adjacent frequencies (p < 0.0001, χ² test, 1 d.f.). The dataset therefore appears to combine low-precision data quantized at magnitude increments of 0.5 or 1.0 with high-precision data at increments of 0.1. The K-S test then understandably rejects the power law when the preponderance of data at increments of 0.5 cannot be explained by chance. Coarse binning with λ = 10 ignores these excesses and the test then accepts the power law with p = 0.73.

Taking λ = 10 still offers ample power to reject alternative distributions, as we reject the power law given fewer bins and far less data in other examples. Therefore, by taking its flaws into account, we conclude that the dataset provides strong evidence that the earthquake magnitudes follow a power law, in agreement with long-held conventional wisdom in geology (table 1).

2.3. Wealth

Vilfredo Pareto—for whom the continuous power law distribution is named—observed power laws in wealth and income as early as 1895 [28]. Measuring wealth remains a difficult and important problem in economics, social science and government [38]. The Forbes 400 list (subsequently The Forbes World’s Billionaires) famously tracks the wealth of the world’s richest people. Prior work, however, found the 2003 dataset inconsistent with a power law [15]. We analyse the dataset’s 261 billionaires, using x_m = 10⁹ and binning by powers of λ = 2 and λ = 4 (table 1 and figure 1b).

The Forbes lists are subject to limitations on accurate measurements of extreme wealth as well as data collection artefacts. The quantization scheme is inconsistent, with values truncated to the nearest 5 million or 100 million depending on whether net worth exceeded one billion. Among the billionaires, all three MLEs and regression give mutually consistent α estimates, which suggests overall conformity to a power law, but K-S tests disagree depending on binning.

Bumps evident in the tail distribution (figure 1b) and QQ plot (electronic supplementary material, figure S2) between 4 billion and 10 billion cause the K-S test on the raw data to reject the power law with p = 0.03, suggesting further inconsistency in representation of values. Indeed, many billionaires are listed has having either 4 or 5 billion while others are distinguished between 4.4, 4.7 or 4.9. In a case of heaping or attraction to particular values, higher-precision values occur noticeably more often just shy of round numbers, for example 6.9, 8.9, 9.7 and 9.8. Any larger amount of binning such as λ = 2 or 4, however, effectively removes these artefacts (figure 1b, dotted lines) and causes the test to accept the power law, vindicating Pareto’s seminal observation.

2.4. Wildfires

A wildfire dataset demonstrates a contrasting example in which binning leads hypothesis tests to reject the power laws when they otherwise would not, despite removing information from the sample. Some models of wildfire spread exhibit power-law scaling, for example through self-organized criticality [39]. However, the power law distribution is not borne out in this empirical data.

The data include 203 784 precise observations of wildfire area over 10 years in the USA recorded consistently to ±0.1 acre with up to six digits of precision over a range from 0.1 to 412 050 acres. A basic visual inspection of the distribution in figure 1c reveals smooth and substantial curvature clearly inconsistent with a power law, even allowing for error and statistical variation. Methods for fitting x_m have nonetheless been devised to locate a fraction of the tail distribution that is consistent [15]. This method finds power-law behaviour only in the highest 0.3% of the data—520 data points.

This x_m = 6324 acres indicates either a cusp beyond which the distribution follows a power law, or else the limit at which there is insufficient statistical power to reject the power law due to too little data. We conduct the inference with this x_m for λ = 1, 2 and 4. The MLEs ${\hat{\alpha }}_1$, ${\hat{\alpha }}_2$ and ${\hat{\alpha }}_4$ give mutually consistent results but differ from the regression estimate, which hints at persistent curvature beyond x_m.

The tests, however, are more likely to reject the power law the coarser the binning. The test accepts the power law in the raw data ($p=0.26$) as stipulated by the calibration of x_m. Binning with increasing λ, however, drives the test to reject the power law, with fine bins giving p = 0.09 and coarse bins ultimately rejecting the power law with p = 0.01. The coarse binning with λ = 4 corresponds to four bins, the last with only one data point. K-S statistics for binned data are therefore capable of rejecting the power law with 10-fold less data and fewer bins than in the earthquake data with λ=10.

While the earthquake and wealth examples show how raw K-S statistics ascribe undue weight to smoothness in the data, conversely, if smoothness is removed from consideration, the test must ascribe more weight to the shape of the distribution over powers of λ. Binning therefore causes K-S statistics to reveal curvature in the empirical distribution by allowing across-scale curvature to stand out from within-scale noise.

3.1. Effects of normal and lognormal noise

We investigate particular cases of normal and lognormal noise where samples from the power-law equation (1.1) with x_m = 1 are convolved with unbiased additive normal noise with variance ${\sigma }_{+}^2$ or multiplicative lognormal noise with variance ${\sigma }_{\times }^2$, then truncated again at x_m for inference. This scheme represents several common sources of error: additive and proportional measurement error, censorship of small values, error in estimating x_m (the tail fraction [24,41]), and contamination from a bulk distribution [42]. The convolution creates measurement error, while the truncation creates both a case of censorship—because data reduced below x_m by noise are removed from the sample—and a bulk distribution distinct from its tail—because applying normal or lognormal error to power law samples yields a unimodal distribution peaked near x_m. Normal error is furthermore a pessimistic case for binning. The quantization and heaping observed in the empirical cases can be completely removed by coarser binning, whereas normal error is always capable of large adjustments with some probability, and so no λ can completely prevent normal error from moving data between bins. In this way, our noisy power laws are chosen as representative and difficult cases of dirty data.

Figure 2 shows the effect of normal and lognormal noise with σ = 0.2 on MLEs and K-S statistics in samples of size 500 with α = 1.5. The noise is barely visually discernible on log–log plots of the tail distribution, but nonetheless biases the raw MLE ${\hat{\alpha }}_1$ by 10–13% on average depending on noise treatment. This bias substantially exceeds the approximately 3.5% standard deviation of the estimator, causing the 95% confidence intervals (CIs) to exclude the true value 68% of the time, and corresponds to a twofold higher prediction for the frequency of extreme events larger than those observed in the dataset.

• Download figure
• Open in new tab
• Download PowerPoint

The theory of maximum likelihood guarantees that ${\hat{\alpha }}_1$ is the most efficient [22] and indeed the distribution of ${\hat{\alpha }}_1$ is the most sharply peaked (figure 2b). However, ${\hat{\alpha }}_1$ is correspondingly the most sensitive to errors while all other estimates are more robust to noise. Coarser binning, with larger values of λ, brings the MLEs ${\hat{\alpha }}_{\lambda }$ closer to the true value on average. For additive noise, binning with λ = 2 approximately halves the bias relative to the standard deviation of the estimator (z-score of the true value) to 1.56. Using λ = 4 further reduces the z-score of the true value to 0.84 and CIs exclude the true value only 15% of the time. In this example, regression minimizes the total root mean squared (RMS) error due to bias and variability, ${({\mathrm{\Sigma }}_i{(\alpha -{\hat{\alpha }}_i)}^2)}^{1/2}$, although regression is known to produce biased estimates in other cases [20,21] such as earthquakes above.

The raw K-S p-values are sensitive detectors of noise (figure 2c). For pure samples, the distribution of p-values is uniform. Uniform p-values control the false positives rate because, for a given stipulated allowable false positive rate ρ, uniformity implies that p < ρ exactly ρ fraction of the time. With noise, however, the p-values reject the power law at ρ = 5% nearly all of the time (additive: 98%, multiplicative: 99%). Binning with sufficiently large λ restores a uniform distribution of p-values (figure 2c), making the hypothesis test robust to noise. With λ = 2, rejection rates fall to 23% and 34%, while with λ = 4 the rejection rates are statistically consistent with $\rho =5\mathrm{\%}$ over 1000 trials (additive noise: 5.4%, multiplicative noise: 6.1%).

The validity of the p-values with and without binning depend on the scientific question. If the question is whether the data follow the continuous power law exactly, then the raw p-values with λ = 1 are valid and the test is adept at rejecting power laws with even small amounts of noise. If the question is whether the data follow a power law allowing for some experimental noise, then the raw p-values fail to control the false positives rate and nearly always falsely reject the power law even when the data do follow a power law. The p-values for binned data are valid for either question by controlling the false positives rate when λ is sufficently large. The distribution of binned data (equations (1.1), electronic supplementary material, equation (S3)) is therefore a better null distribution for the more practical latter question.

A classical strategy for removing the effects of censorship of small values, contamination from a bulk distribution, and error in estimating the tail fraction or x_m is to simply increase x_m, as more extreme values are more likely to follow the limiting tail distribution and larger values are less likely to be contaminated by error [15,24,41,42]. This strategy also applies here, and so we compare its effectiveness against binning with x_m fixed at 1. Fitting x_m by minimizing K-S distance, the power law [15,41] selects x_m = 1.4 ± 0.8 ($\pm 1.96\hspace{0.17em}\mathrm{s}.\mathrm{d}.$) and ${\hat{\alpha }}_1=1.52\pm 0.21$ for additive error and x_m = 1.5 ± 0.9, ${\hat{\alpha }}_1=1.50\pm 0.23$ for multiplicative error (figure 2b). This strategy minimizes bias by throwing away the most contaminated data, typically 20–25%, but comes at the cost of nearly doubling the variability of the estimator. Binning with λ = 4 and x_m = 1 performs slightly better in representing the true value—with RMS error of 0.09 and 0.11 for additive and multiplicative noise in contrast to 0.11 and 0.12 when fitting x_m. Fitting x_m, on the other hand, minimizes bias. Neither method alone outperforms linear regression by either standard in this example.

3.2. Power-accuracy trade-off in maximum-likelihood estimators

Much as increasing x_m removes potentially erroneous data at the cost of statistical power, binning also removes information from the sample, inducing a trade-off between accuracy and power in estimation. Increasing λ increases variability on the estimator (electronic supplementary material, S1),

but potentially reduces unknown bias, depending on errors in the data. The binning ratio λ can therefore trade unknown bias for a calculable increase in variability. Better trade-offs are achievable with larger sample size, as larger samples both reduce variability on the estimator and expand the range of the data, allowing larger λ. Increasing λ with sample size, by holding variability constant, therefore allows larger samples to decrease bias due to error, in contrast to the typical cases of estimation in which bias is independent of sample size. We investigate the trade-off numerically in samples of size n = 1 000 000 (200 replicates; α = 1.5; σ₊, ${\sigma }_{\times }=0.02,0.2$). Figure 3 shows a clear trade-off between accuracy and precision as a function of λ in our normal and lognormal noisy power-law data.

• Download figure
• Open in new tab
• Download PowerPoint

We then compute the λ_RMS that minimizes the RMS error in representing the true value for each noise treatment. This λ_RMS ranged from 20 to 330, corresponding to dividing the data into two to four bins. In comparison with the proportional range of the data r := max_i x_i/x_m, the fraction of log-range covered by the first bin log λ_RMS/log r varied only from 0.32 to 0.61 (median r: 12 500, $\sqrt{r}\hspace{0.17em}:\hspace{0.17em}111$). That is, the λ_RMS that best represents the true value typically divided the data into two or three bins despite a 10-fold difference in error magnitudes.

3.3. Sensitivity and specificity of Kolmogorov–Smirnov tests

Hypothesis testing generally involves a trade-off between sensitivity and specificity. The test should reject the null hypothesis when it is false (sensitivity) and accept it when it is true (specificity). Statistical variability makes perfect performance on both measures impossible, whereas sensitivity often can be increased at the cost of specificity and vice versa. While greater λ generally increases specificity by ignoring errors in power-law data, it is possible for greater λ to increase sensitivity as well, as observed in the wildfires example. Here, we show in simulations of power law, noisy power law and lognormal distributions that λ may increase specificity with no cost to sensitivity, or else tune a trade-off between sensitivity and specificity (figure 4).

• Download figure
• Open in new tab
• Download PowerPoint

We evaluate the performance of hypothesis testing using lognormal samples truncated below x_m = 1 as an alternative, non-power-law distribution. The lognormal distribution is heavy-tailed and so its extremes can be difficult to distinguish from a power law [15]. Both tail distributions can appear straight on a log–log plot, but the lognormal has curvature that depends on its parameters (figure 4a). The truncated lognormal probability density function is

where l(x) is the usual lognormal density function and Z is the normalization ${\int }_{x_m}^{\mathrm{\infty }}l(u)\hspace{0.17em}\mathrm{d}u$. Its three parameters are the mean μ, variance ${\sigma }_{\ln }^2$ and tail threshold x_m. As σ_ln goes to infinity, the curvature diminishes and f(x) converges to a power law with a fixed α = 1. The lognormal tail can therefore approximate a power law with any α > 1 arbitrarily well for sufficiently large but finite σ_ln. We set x_m to 1, use σ_ln to adjust curvature, and finally choose μ so that the slope is equal to α = 1.5 at x_m to approximate the small power law samples of figure 2 where n = 500.

We find that binning entails no loss of sensitivity or statistical power, provided λ divides the data into at least four bins (figure 4b). The sensitivity itself depends strongly on the curvature of the lognormal. At σ_ln = 1, the most curved, sensitivity is approximately 90% whereas at σ_ln = 2, 28% is the best attainable. However, sensitivity was independent of λ for λ < 2 for any given curvature. The median range of all samples was 12.8, and so λ = 2 typically puts the maximum data point in the fourth bin, [8,16) (figure 4a).

Increasing λ removes noise, and so increases the specificity of the test, when we take rejections of the power law due to noise as false positives. The specificity given pure continuous power-law data is always 1 − ρ = 0.95, whereas specificities on noisy data are not always 0.95, but instead depend on λ due to noise biasing K-S p-values. The range 2 < λ < 4 shows a trade-off between sensitivity and specificity (figure 4b), in which increasing λ decreases sensitivity to reject lognormal tails and increases specificity for accepting comparable noisy power-law samples with ${\sigma }_{+},{\sigma }_{\times }=0.2$ (cf. electronic supplementary material, section (a)).

Overall, increasing λ increases specificity at no cost to sensitivity up to a point. This much is always desirable. Further increases in λ sacrifice sensitivity for specificity, which may be desirable, up to such λ as the test becomes useless because it is unable to reject the power law more often than ρ under any circumstances. In general, these thresholds of λ depend on sample size, type and magnitude of errors in the data, and possible alternative distributions, but reasonable guesses can be made based on number of bins.

4.1. Choosing values of λ

Considerations involved in choosing values of λ include the nature of the errors in the data, possible alternative distributions, the sample size and the scientific question. For example, binning should always respect the quantization in the data, if present. If the data are logarithmically quantized at some λ₀ as in the earthquake data, only integer powers of the quantization, $\lambda ={\lambda }_0^k$, should be considered. If the data are quantized to a linear or inconsistent scale, λ should be chosen large enough to make the differences negligible, or more complicated binning schemes should be devised (cf. [34]).

Estimating α requires λ to be within the proportional range of the data r, so that the MLE ${\hat{\alpha }}_{\lambda }$ is well defined (electronic supplementary material, equation (S4)). At the two extremes of the range 1 < λ ≤ r, ${\hat{\alpha }}_1$ maximizes bias while minimizing statistical variability, whereas ${\hat{\alpha }}_r$ minimizes bias but maximizes variability. Neither extreme is desirable unless the data is perfect, in which case there is no bias and ${\hat{\alpha }}_1$ minimizes variability. On imperfect data, equally dividing the logarithmic range of the data using $\lambda =\sqrt{r}$ balances bias against variability to some degree, and is reasonable unless more detailed information is available. Indeed, in the numerical experiments of figure 3, increasing noise by a factor of 10 only doubles log λ_RMS/log r, so that $\lambda =\sqrt{r}$ achieves comparable performance to λ_RMS over a wide range of error magnitudes. This λ corresponds to a minimum of bins, however, which offers the hypothesis test no information, and so these estimates are predicated on other belief that the data do indeed follow a power law. Other reasonable prescriptions for estimation are to minimize bias subject to an allowable statistical error using equation (2.1) or choosing λ to greatly exceed a known magnitude or proportional error.

Hypothesis testing imposes different limits on λ because more bins are required to distinguish the shapes of alternative distributions. Increasing λ reduces the chance of a false rejection due to trivial errors, but also might reduce capacity to reject alternative distributions. Numerical experiments indicate that choosing λ < r^1/4 entails no loss of power to reject lognormal data, whereas choosing λ > r^0.4 offers no power to reject alternatives. Binning the data into at least five bins by choosing λ up to r^1/4 is therefore a prescription for detecting curvature over the range of the data, whereas detecting deviations over a limited range such as roll-off [31,44] or censorship of small values [37] would require more bins within the interval of interest. Coarser binning than r^1/4 may sacrifice sensitivity for a possible increase in specificity, but binning beyond λ = r^0.4 is likely to have no sensitivity and is therefore useless.