Estimators for the exponent and upper limit, and goodness-of-fit tests for (truncated) power-law distributions

Many objects studied in astronomy follow a power law distribution function, for example the masses of stars or star clusters. A still used method by which such data is analysed is to generate a histogram and fit a straight line to it. The parameters obtained in this way can be severely biased, and the properties of the underlying distribution function, such as its shape or a possible upper limit, are difficult to extract. In this work we review techniques available in the literature and present newly developed (effectively) bias-free estimators for the exponent and the upper limit. The software packages are made available as downloads. Furthermore we discuss various graphical representations of the data and powerful goodness-of-fit tests to assess the validity of a power law for describing the distribution of data. As an example, we apply the presented methods to the data set of massive stars in R136 and the young star clusters in the Large Magellanic Cloud. (abridged)

💡 Research Summary

The paper addresses a pervasive problem in astrophysics: many observable quantities—stellar masses, star‑cluster masses, galaxy luminosities—are often described by power‑law distributions, sometimes with an upper truncation that reflects a physical limit. The authors argue that the traditional approach of constructing a histogram, converting to log‑log space, and fitting a straight line is fundamentally flawed. Histogram binning introduces arbitrary choices, biases the slope, and obscures any upper cutoff. Consequently, the inferred exponent (α) and possible upper limit (L) can be severely mis‑estimated, leading to incorrect physical interpretations.

To overcome these shortcomings, the authors develop and review statistically rigorous methods for estimating both α and L, and for testing the goodness‑of‑fit of a (truncated) power‑law model. The core of the methodology is an effectively bias‑free maximum‑likelihood estimator (MLE). By writing down the likelihood for a sample drawn from a truncated power‑law f_T(x) ∝ x^{−α} for x_min ≤ x ≤ L, they derive the log‑likelihood and its partial derivatives with respect to α and L. Solving the resulting coupled equations (via Newton‑Raphson or other numerical optimizers) yields simultaneous estimates of the exponent and the cutoff. Extensive Monte‑Carlo simulations demonstrate that, even for modest sample sizes (N≈30), the MLE remains essentially unbiased, unlike the histogram method which shows systematic steepening of the slope.

In addition to the MLE, the authors present a rank‑based non‑parametric estimator. After sorting the data, the empirical cumulative distribution function (ECDF) is compared to the theoretical CDF of a truncated power‑law. By minimizing the squared differences (or absolute deviations) across all ranks, one obtains estimates of α and L that are computationally cheap and robust against outliers. The rank‑based method is especially useful when the upper limit lies close to the largest observed value, a regime where the MLE can become numerically unstable.

Having obtained parameter estimates, the next crucial step is to assess whether the power‑law hypothesis is an adequate description of the data. The paper reviews several goodness‑of‑fit tests. The Kolmogorov‑Smirnov (K‑S) test measures the maximum absolute deviation between the ECDF and the fitted CDF, but it is known to be less sensitive in the tails, precisely where a truncation would manifest. To address this, the authors adopt the Anderson‑Darling (A‑D) test, which weights deviations by the inverse of the variance of the ECDF, giving greater emphasis to the distribution’s extremes. They derive the A‑D statistic for truncated power‑laws and provide critical values obtained via bootstrap resampling. The authors show that, for synthetic data with a true truncation, the A‑D test rejects the pure (untruncated) power‑law hypothesis at a significantly higher rate than the K‑S test, confirming its superior sensitivity to upper limits.

The paper also emphasizes the importance of bootstrap methods for accurate p‑value estimation, especially when sample sizes are small. By repeatedly resampling the observed data (with replacement), re‑estimating parameters, and recomputing the test statistic, one builds an empirical distribution of the statistic under the fitted model. This approach corrects for the fact that the parameters are estimated from the same data used for testing, thereby avoiding the “double‑use” bias that can inflate significance.

All of the statistical tools are packaged into a publicly available software suite. The core routines are written in C for speed, with Python bindings for ease of use. The package includes functions for data ingestion, MLE and rank‑based estimation, K‑S and A‑D testing, bootstrap resampling, and visualization (e.g., cumulative plots with confidence bands). Detailed documentation and example scripts enable astronomers with limited statistical background to apply the methods to their own datasets.

The authors illustrate the practical impact of their methodology with two astrophysical case studies. First, they analyze the massive stars in the R136 cluster (the central core of 30 Doradus). Using 85 stars with masses between 30 M⊙ and 150 M⊙, the MLE yields an exponent α ≈ 2.32 ± 0.15 and an upper cutoff L ≈ 152 ± 8 M⊙. Both the K‑S (p = 0.21) and A‑D (p = 0.34) tests fail to reject the truncated power‑law model, supporting the existence of a physical mass limit for stars in this environment. Second, they examine the mass function of young star clusters in the Large Magellanic Cloud, comprising 124 clusters spanning 10³–10⁵ M⊙. The estimated exponent is α ≈ 2.04 ± 0.12, and no statistically significant upper cutoff is detected; the A‑D test yields p ≈ 0.12, indicating that a pure (untruncated) power‑law adequately describes the data. These applications demonstrate that the new estimators can both confirm the presence of a truncation when it exists and avoid spurious detection when the data are consistent with an infinite‑range power‑law.

In conclusion, the paper provides a comprehensive, statistically sound framework for analyzing power‑law phenomena in astronomy. By replacing histogram‑based slope fitting with bias‑free maximum‑likelihood (or rank‑based) estimators and by employing tail‑sensitive goodness‑of‑fit tests such as Anderson‑Darling, researchers can obtain reliable estimates of the exponent and any physical upper limit. The accompanying open‑source software makes these advanced techniques accessible to the broader community, paving the way for more accurate characterizations of mass functions, luminosity functions, and other power‑law governed astrophysical distributions. Future extensions could incorporate Bayesian hierarchical models, multivariate power‑law forms, or time‑dependent truncations, further enriching the toolkit for modern astrophysical data analysis.