Statistical characteristics of earthquake magnitude based on the composite model

1.   Introduction

Earthquakes with a magnitude greater than five cause many casualties, and destructive earthquakes are a threat to human life and property. Many small and medium-sized earthquakes cause more disasters than people expect, even affecting social stability, and are of great concern to the government and the public. Therefore, small and medium-magnitude seismic data study is also of practical importance. The Peaks Over Threshold (POT) model based on the generalized Pareto distribution (GPD) is widely used in data analysis in various industries because it can fully use the information of extreme value contained in the data ^[1]. In recent years, there have been some findings of earthquake hazards, especially in the analysis of large-magnitude earthquakes ^[2]. The classical fixed threshold modeling approach uses graphical diagnostics. Coles ^[3] outlined the common graphical diagnostics for threshold choice. Scarrott and MacDonald ^[4] gave a long review of tail and non-tail estimation methods. However, this model falls short in practical data processing in two ways as data analysis with such a model is generally done in two steps. In the first, the threshold u is chosen graphically by looking at the mean excess plot ^[5] or simply setting it as some high percentile of the data ^[6]. Both methods are subjective and the parameter estimates depend on the threshold. In such areas, the main problem is the scarcity of data or, more specifically, modeling with a fairly small amount of observations.

To address this issue, Arnoldo et al. ^[7] proposed a new dynamically weighted mixture model, where one of the terms is the GPD and the other is a light-tailed density function. Mendes and Lopes ^[8] proposed a procedure to fit a mixture model by maximum likelihood where the tails are GPD and the distribution center is normal. The normal spliced with GPD tail developed by Carreau and Bengio ^[9] named the 'hybrid Pareto' model was further developed to include constraints on parameters to ensure continuity up to the first derivative of the density. Behrens ^[10] proposed a model to analyze data characterized by extreme events where a threshold is estimated directly. Nadarajah and Bakar ^[11] proposed new composite models based on the lognormal distribution. Caladerin-Ojeda ^[12,13] analyzed claim data and all French settlements from 1962 to 2012 using the composite Weibull-Burr model and the composite lognormal-Pareto distribution separately. Extreme likelihood estimation is widely used in composite models ^[14]. Carreau & Bengio ^[15] and Carreau et al. ^[16] implemented a neural network learning approach in nonstationary and bivariate modeling situations. Bayesian inference was used by Behrens et al. ^[10] with sensible prior forms for the bulk, tail and threshold parameters. Tancredi et al. ^[17] was the first to propose an extreme value mixture model that combined a nonparametric estimator for the bulk distribution spliced with an extreme value tail model. Li et al. ^[18] extended the model by fitting the tails with GPD and the central part with gamma, lognormal, mixed gamma and Weibull distributions to build four combined models for the evaluation of post-earthquake loss in Yunnan Province, respectively. Due to the complexity of the combined models, the Bayes estimates of the parameters were calculated using the Markov chain Monte Carlo (MCMC) method used in the previous papers. There are other methods to estimate the parameters. Dupuis ^[19] proposed a robust procedure for fitting GPD, including statistics to guide threshold choice. Thompson et al. ^[20] developed an automated procedure for threshold estimation and uncertainty quantification. They set a uniformly spaced grid of possible threshold values (between the median and the 98% empirical quantile). For each potential threshold, the GPD is fitted (using maximum likelihood (ML) estimation) and the differences in the modified scale parameters for neighboring thresholds are calculated.

The composite model allows fitting all the data and considering the thresholds as parameters for direct estimation. We propose an empirical distribution function-based parameter estimation method to estimate the parameters in the composite model. The correctness of the parameter estimation method is tested by generating simulated data. Finally, the model and parameter estimation method is used to analyze the seismic data of the fracture zone with the composite model.

2.   Proposed composite models

Consider that ${X_1}, {X_2}, \cdots, {X_n}$ are independent and identically distributed observations and $u$ is the threshold over which these observations are exceedances.

Then, $({X_i}|{X_i} \geqslant u)\sim G(\xi, \sigma, u)$, where

is the distribution function of the GPD and the probability density function of the GPD is

Observations below the threshold obey the $H$ distribution, which can be estimated.

Assume that is any distribution such as a Weibull, Gamma or Lognormal distribution. The distribution function of the composite distribution is as follows:

Composite model one Gamma-GPD composite model:

In formula (3), if , , where

is the probability density function of the Gamma distribution. The distribution function of the composite model is

Composite model two Weibull-GPD composite model:

In formula (3), if , , where

is the probability density function of the Weibull distribution. The distribution function of the composite model is

Composite model three Lognormal-GPD composite model:

In formula (3), if , , where

is the probability density function of the Weibull distribution. The distribution function of the composite model is

3.   Parameter estimation of composite models

3.1.   Parameter estimations based on the empirical distribution functions

Definition 1. A statistical estimation of $F(x)$ based on a random sample ${X_1}, {X_2}, \cdots, {X_n}$ is the empirical distribution function defined by

where $I(\cdot)$ denotes the indicator function of the event in the brackets.

From the Glivenko-Cantelli theorem ^[21], it follows

Therefore, when the sample size $n$ is sufficiently large, ${F_n}(x)$ approximates well the distribution function , where is the parameter of the distribution. Based on the above conclusions, we devise the following method for estimating the unknown parameters of the distribution. Define the loss function as

where ${F_n}(x)$ is the empirical distribution function.

Find that satisfies and is an estimate of the parameter . This method allows for the numerical estimation of complex distributions that contain several parameters. For the composite model in this paper, which contains several parameters and cannot be solved analytically, we use the above empirical distribution function-based method to estimate the parameters.

To illustrate the feasibility of the method, we compare the traditional maximum likelihood estimation with the estimation based on the empirical distribution function.

3.2.   Maximum likelihood estimation of parameters

The likelihood function of the composite model is

where

is the probability density function of GPD. Find that satisfies and is the maximum likelihood estimate of the parameter . Since there is no analytical solution to the maximum likelihood estimate, the numerical solution is solved here using the optimization algorithm.

3.3.   Numerical simulation

The numerical simulation process of parameter estimations based on the empirical distribution functions is as follows and the flowchart is shown in Figure 1.

(1) For the given parameters, define the distribution function of the composite model and generate random numbers of [0, 1]. The distribution functions of the three composite models are given in section two. Apply the inverse transformation method to obtain samples obeying the composite models.

(2) Calculate the empirical distribution function ${F_n}(x)$ of the samples. Define the loss function as

(3) Apply the Nelder-Mead simplex algorithm in Matlab to find the minimum of the loss function and get the parameter estimations.

(4) The number of samples generated is n = 100, n = 200, n = 300, n = 500 and n = 1000 and the mean of the estimates obtained by repeating m = 100 times, respectively, is the parameter estimate.

(5) Apply nonparametric bootstrap methods to obtain confidence intervals for parameters.

The numerical simulation process of parameter estimations based on the maximum likelihood method is similar to the empirical distribution function based method.

In the composite model ${\bf{x}} = ({x_1}, {x_2}, \cdots, {x_n})$ for the given parameters ${\bf{ \pmb{\mathsf{ θ}} }}$, the numerical simulation results are given in Tables 1–3 as follows. The histograms of the parameters of each model are shown in Figures 2–4.

From Tables 1–3, it can be seen that for all three composite models, the proposed method based on empirical distribution functions is similar to the maximum likelihood estimation and is a feasible parameter estimation method. Both parameter estimation methods are relatively robust. Except for the parameter $\xi$, the deviation of each parameter estimate from the true value is small. The standard deviation and 95% confidence intervals for the parameters are shown in Tables 4–8.

The standard deviations of the parameters $\lambda$ are smaller in model 2, which is consistent with Figure 3. The standard deviation of the parameters $\xi$ is greater in model one. The method appears to have a variability with $\xi$ but is robust for the other parameters.

A comparison of the empirical distribution functions with the distribution function of the three composite models is shown in Figures 5–7. The two curves are very close to each other in each figure. The confidence interval for the empirical distribution function is also shown in each figure. The distribution function falls between the upper and lower confidence levels of the empirical distribution function.

3.4.   Seismic risk based on the composite model

GPD in the composite model can describe the tail data. The quantile is important for extreme value analysis. Scholars ^[22,23,24] proved that if the distribution function $F$ belongs to the attraction domain of a generalized extreme value distribution (GEV), the limiting distribution of excesses is the GDP. The distribution function of the excess is

where,

and

From Eq (15) ${\widehat {\overline F }_u}(x) = {(1 + \xi \frac{x}{\sigma })^{ - \frac{1}{\xi }}}$, the frequency ${N_u}/n$ above the threshold $u$ is the estimation of $\overline F(u)$, where ${N_u}$ is the number of samples above the threshold, then $\widehat {\overline F }(u) = \frac{{{N_u}}}{n}$. We can get the estimation of $\overline F(u + x)$ as follows:

then we can get the estimation of $F(x)$:

From Eq (17), we get the estimation of the p-quantile ${x_p}$:

where $F({x_p}) = p, 0 < p < 1$.

If $\xi < 0$ and $p \to 1$, the estimation of the upper limit point ${x^*}$ of the support $F$ is

$P\{ X > {x_p}\} = 1 - F({x_p}) = 1/T$, $T = 1/\overline F({x_p})$ is the theoretical return period with the return level ${x_p}$ and $F({x_p}) = p{\text{ }}(0 < p < 1)$. The return period estimate can be obtained:

If the given return period is $T$ years, bring $p = 1 - 1/(365T)$ into Eq (18), and we can estimate the return level

In the following part, the earthquake magnitude data is analyzed using the above three composite models to study the seismic risk.

4.   Earthquake magnitude analysis of the eastern Bayan Hara Block

The data in this paper is from the National Seismic Science Data Sharing Center (https://data.earthquake.cn/). 19,221 records before December 2019 are the research samples in the scope of reference ^[25]. The annual magnitudes of the earthquake are shown in Figure 8.

The magnitude data is analyzed with the above three composite models. The final results of the parameters obtained by the empirical distribution function-based method are summarized in Tables 9–11.

Although the distributions describing the central part are different, the parameter estimations of the GPD that describe the tail are very close to each other in the three composite models. The values of the loss functions are 1.0212, 1.4330, 1.5874 in the three models, respectively, and the residual of composite model one is the smallest. From Figure 9 we can see that the Gamma distribution function is closer with the first half of the empirical distribution function, which is the same as the loss function value. Akaike information criterion (AIC) values are also reported in Tables 9–11 to compare the performance of different models. AIC values are 10.388 in composite model one and 14.265 in composite model three. The AIC value of 311.5178 for model two is much bigger than that for models one and three. Thus, the performance of models one and three are comparable, and model two performs the least. The value of the loss function for model one is the smallest; hence, model one is selected as the reasonable model.

Nowadays, there are non-parametric and semi-parametric approaches for the distribution of the non-tail, which make modeling more flexible. Mohd et al. ^[26] used the semi-parametric for measuring income inequality in Malaysia. Based on this, we developed a semi-parametric model for comparison with the composite model.

Assume that $p$ denotes the proportion of top earthquake magnitude data. In the semi-parametric model, a generalized Pareto distribution is fitted to $100p\%$ of the upper tail data.

The $100(1 - p)\%$ lower tail data is modeled by an empirical distribution. The empirical distribution is given by

where $I(\cdot)$ denotes the indicator function of the event in the brackets, and k is the number of observations in the upper tail of the earthquake magnitude.

The full semi-parametric distribution can be written as

where $G(\xi, \sigma, u)$ is the distribution function of the GPD.

We used the same empirical-distribution function-based method to calculate parameters. The result is in Table 12 below.

The parameter $u$ is the threshold of earthquake magnitude; therefore, the estimation value of 7.0469 is too large and unreasonable. The data above the threshold is only 15. The sample size is reduced and, consequently, the variance of the parameter estimates increases.

Therefore, model one can be chosen to study the distribution of earthquake magnitude. The earthquake magnitude data that is smaller than the threshold obeys the Gamma distribution with parameters $\alpha = 5.7666, \beta = 1.4296$, while the tail data larger than the threshold obeys a GPD. This conclusion can be used as a complement to probabilistic seismic hazard analysis (PSHA).

The estimation of $\widehat \xi$ in Tables 9–11 is less than zero, indicating that the return level of magnitude predicted by the model has a theoretical upper limit. The theoretical maximum magnitude is about ${M_s} = 8.4$ according to Eq (19) in composite model one.

In the Pareto distribution, the most important part of the tail analysis is the estimation of the quantile. The quantile is the return level of the magnitude of earthquake. We obtain the $1 - p$ quantile of the distribution according to Eq (18). Bring the estimates of $u, \sigma$ and $\xi$ into Eq (21), and the return level for the given return period can be summarized in Table 13.

The 100-year return level is ${M_s} = 8.24$, indicating that earthquakes with a magnitude of about ${M_s} = 8.2$ will occur in the east of Bayan Har every 100 years. This result is consistent with the fact that there was an earthquake of ${M_s} = 8.0$ in 2008. It shows that the estimated return levels are all within the confidence intervals of the return level in Figure 10. The model is reasonable.

5.   Conclusions

In response to the difficulty of selecting a threshold when analyzing data using the GPD and the fact that the model can only study data above the threshold, this article studied all the data using composite models. Composite models were developed by applying the Gamma distribution, the Weibull distribution, and the lognormal distribution combined with the GPD separately. A parameter estimation method was proposed based on the empirical distribution function. The parameter estimation method based on the empirical distribution function is suitable for the loss function, which is complex and can not be derived to optimize the function. The method appears robust for most of the measurements in the three composite models in this paper. For some parameters, the method appears to have some variability. This method is suitable for solving numerical solutions of parameter estimates in complex distributions. The composite models and the empirical distribution function-based method are applied to earthquake magnitude data to provide a reference for earthquake hazard analysis. The composite model can analyze all the data and the threshold can be directly estimated as a parameter. Composite models in seismic magnitude analysis expand the scope of their use and provide new ideas for magnitude analysis. Composite models and parameter estimation methods can also be used in other industries.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This research was funded by the Fundamental Research Funds for Central Universities, Grant number ZY20215140 and Self-Financing Project of Scientific Research and Development Plan of Lang Fang Science and Technology Bureau, Grant number 2022011019.

The authors thank all the reviewers for their helpful comments.

Conflict of interest

The authors declare no conflict of interest.