New Statistical Algorithm Sharpens Estimates of Extreme Earthquake Risk

Earthquakes are among the most stubbornly unpredictable phenomena in nature, and the mathematics used to describe them has long struggled to keep pace with their extremes. Now, two Algerian mathematicians have developed a practical computational method that makes a powerful but awkward probability distribution far easier to use, potentially improving how seismologists quantify the risk of the largest tremors. In a study published in the journal Advances in Statistical Climatology, Meteorology and Oceanography, Fatah Benatia of Mohamed Khider University in Biskra and Mohammed Ridha Kouider of the Higher School of Social Security in Algiers present an efficient algorithm for computing maximum likelihood estimates of the parameters of the q-Pareto distribution, a heavy-tailed distribution that has attracted growing attention for modeling the magnitudes of large earthquakes.

The q-Pareto distribution, or q-PD, belongs to a broader family of q-exponential distributions that emerged from the work of Constantino Tsallis, who introduced the q-logarithmic and q-exponential functions in the 1990s as generalizations of the ordinary logarithm and exponential. These functions arise naturally in the study of complex systems, where they maximize a generalized entropy known as Tsallis entropy, extending the classical Boltzmann-Gibbs framework that governs systems at equilibrium. The single parameter q controls how slowly the tail of the distribution decays, and this is precisely what makes the q-PD attractive for phenomena dominated by rare, extreme events. When q equals one, the familiar exponential behavior is recovered; when q differs from one, power-law-like tails emerge, allowing small changes in parameters to translate into enormous changes in the predicted frequency of catastrophic outcomes.

Heavy-tailed distributions of this kind are not merely academic curiosities. Earthquake magnitudes are measured on logarithmic scales, meaning that each incremental increase in magnitude corresponds to roughly a tenfold increase in shaking amplitude and about a thirty-two-fold increase in the energy released. The Richter scale, developed in 1935, works reasonably well for small to moderate local events but loses accuracy for the largest earthquakes, exactly where accurate statistical modeling matters most. Distributions related to the Pareto family, including the log-Pareto, tapered Pareto, and generalized Pareto distributions, have therefore become standard tools for describing the energy released by major seismic events. The q-PD, proposed in its current form by de la Barra and Vega-Jorquera in 2021, has proven especially effective for modeling the magnitudes of the largest, most extreme earthquakes in global catalogues, and it also finds use in insurance and actuarial statistics for pricing reinsurance premiums against catastrophic losses.

The theoretical foundation of the new work rests on a deep connection between the q-PD and extreme value theory. A cornerstone result of extreme value theory, established independently by Balkema and de Haan in 1974 and by Pickands in 1975, states that for data drawn from almost any unknown distribution, the values that exceed a sufficiently high threshold converge, in distribution, to the generalized Pareto distribution. This is the mathematical justification for the peak-over-threshold method used across hydrology, climatology, finance, and geophysics. Benatia and Kouider prove a theorem showing that the generalized exponential distribution they define, obtained by transforming the q-exponential function, overlaps exactly with the generalized Pareto distribution under a simple reparameterization. Because the q-PD is linked to this generalized exponential distribution through an exponential transformation of the data, the q-PD inherits the same status as a limiting distribution for threshold exceedances, making it a legitimate and principled model for extreme earthquakes rather than an ad hoc curve fit.

Estimating the two parameters of the q-PD, a positive shape parameter q and a positive scale parameter b, is where the practical difficulty lies. Maximum likelihood estimation is the gold standard in statistics because it produces consistent estimators with asymptotically normal distributions and is asymptotically efficient in many settings. However, the likelihood equations for the q-PD do not admit closed-form solutions. Setting the partial derivatives of the log-likelihood to zero yields a nonlinear equation in a transformed parameter that must be solved numerically, and the solution space is constrained: the shape parameter must remain strictly positive, and the admissible range of the data depends on the unknown parameters themselves. The authors show that when the shape parameter is negative, the density becomes unbounded and the standard regularity conditions underpinning maximum likelihood theory fail, so the search must be confined to a carefully defined feasible region.

To navigate these complications, the researchers built their algorithm around a modified bisection algorithm for multiple roots, a numerical technique previously developed by Kouider. Compared with the more commonly used combination of bisection and Newton’s method, the modified bisection approach offers distinct advantages: it does not require computing derivatives of the objective function at reference points, and it does not depend on an initial solution, a factor that can cause Newton’s method to converge to spurious points if the starting values are poorly chosen. The authors prove a theorem that provides explicit lower and upper bounds for the transformed parameter, derived using Jensen’s inequality and properties of the order statistics of the sample. These bounds guarantee that the root-finding procedure operates within a known interval, and a small tolerance parameter, scaled to the sample, keeps the iterative solutions from straying outside the feasible set. Once candidate roots are located, the algorithm evaluates the log-likelihood at each corresponding parameter pair, distinguishes genuine local maxima from boundary maxima, and selects the pair that achieves the highest likelihood.

The algorithm also delivers something that raw point estimates cannot: uncertainty quantification. By computing the second-order partial derivatives of the log-likelihood, the authors assemble the observed Fisher information matrix, whose inverse provides the asymptotic variance-covariance structure of the estimators. For samples larger than thirty observations, confidence intervals are constructed using quantiles of the standard normal distribution; for smaller samples, the algorithm switches to Student’s t-distribution with the appropriate degrees of freedom. This built-in machinery for confidence intervals is essential for applied seismologists, who need to know not just a best guess for the tail behavior of earthquake magnitudes but also how reliable that guess is given limited data.

To validate the method, the authors ran a simulation study using the R statistical software. They generated a sample of fifteen values from a q-Pareto distribution with known parameters, choosing a shape parameter of 0.5, a scale parameter of 1, and a threshold of 1, by inverting the distribution function applied to uniform random numbers. Applying their algorithm to this synthetic data recovered maximum likelihood estimates of approximately 0.504 for the shape parameter and 1.024 for the scale parameter, remarkably close to the true values used to generate the data. The corresponding 95 percent confidence intervals, spanning roughly 0.039 to 0.969 for the shape parameter and 0.369 to 1.680 for the scale parameter, comfortably contained the generating values, demonstrating that the procedure is both accurate and honest about its uncertainty even with a very small sample.

The real-world test came from tectonic earthquake activity in and around the British Isles. Using a catalogue of 38 earthquake exceedances recorded over a sixty-day window from late March to late May 2024, drawn from data published by the British Geological Survey, the authors applied their algorithm to the logarithms of the recorded values. The analysis yielded a maximum likelihood shape parameter estimate of about 0.477 and a scale parameter estimate of about 0.557, with 95 percent confidence intervals of roughly 0.211 to 0.743 and 0.345 to 0.768 respectively. A shape parameter well below one indicates a heavy but ultimately manageable tail, quantifying how quickly the probability of larger events decays with magnitude in this region. The authors note that earthquake intensity measurement itself carries subtleties, since intensity scales such as the Modified Mercalli Inventory assess shaking and damage at specific locations rather than the energy emitted at the source, and they treat the threshold selection accordingly in their analysis.

Beyond the immediate application to seismology, the study lays groundwork that other researchers can build upon. Because the generalized exponential distribution the authors define coincides with the generalized Pareto distribution under their reparameterization, the same estimation machinery can be used to estimate the extreme value index, the fundamental quantity that governs how heavy the tail of any underlying distribution is. The approach also connects to a growing body of work on estimating generalized Pareto parameters under difficult conditions, including random censoring, which the authors and collaborators have explored in earlier publications. As extreme events, from earthquakes to floods to financial crashes, increasingly dominate public risk assessments, tools that make sophisticated heavy-tailed modeling fast, robust, and accessible could prove indispensable. The q-Pareto distribution, once a theoretical curiosity from the mathematics of complex systems, is now a practical instrument for reading the seismic signature of the planet’s most violent moments.

Subject of Research: Maximum likelihood estimation of q-Pareto distribution parameters and its application to modeling extreme earthquake magnitudes

Article Title: On computing the maximum likelihood estimates for the q-Pareto distribution and application to earthquakes

Article References: On computing the maximum likelihood estimates for the q-Pareto distribution and application to earthquakes. (n.d.). https://doi.org/10.5194/ascmo-12-211-2026

Image Credits: AI Generated

DOI: 10.5194/ascmo-12-211-2026

Keywords: q-Pareto distribution, maximum likelihood estimation, extreme value theory, earthquakes, seismology, heavy-tailed distributions, Tsallis entropy, generalized Pareto distribution, numerical algorithms, peak-over-threshold, statistics, risk modeling