The earthquake is the supreme terrifying and harsh phenomena of nature that can do significant damages to infrastructure and cause the death of people. Nepal is one of the paramount catastrophe prone countries in the world. Nepal situated in the center of the Himalayan range, lies in between 80˚4' to 88˚12' east longitude and 26˚22' to 30˚27' north latitude (MoHA & DP Net, 2015) . Less than 10% of earthquakes happen within seismic plates, but remaining 90% are commonly found in the plate periphery (Lamb & Jones, 2012) . The entire region of Nepal is likely to experience devastating earthquakes as it lies between two seismically energetic Indian and Eurasian tectonic plates (MoUD, 2016) .
Nepal has a long history of numerous earthquakes and has experienced great earthquakes in the past two centuries with moment magnitudes Mw = 7 and greater. The recorded earthquake in the history of Nepal was on 7th June 1255 AD with magnitude Mw = 7.7. On 16th January 1934 AD, an earthquake called Nepal Bihar Earthquake, hit Nepal and its surrounding regions with Mw = 8.4 magnitude. The fatality figures were the highest for any recorded earthquake in the history of Nepal (MoHA & DP Net, 2015; MoUD, 2016) . The latest earthquake experienced in Nepal was on 25th April 2015 at 11:56 am local time. The earthquake of magnitude 7.8 Mw, called Gorkha Earthquake, hit at Barpark located 82 kilometers northwest of Nepal’s capital of Kathmandu affecting millions of citizens (USGS, 2016) . This event has been the most powerful earthquake disaster to strike Nepal since the earthquake in 1934, tracked by many aftershocks, the largest being Mw = 7.3 magnitude on 12th May 2015. The devastating earthquake included about 9000 fatalities, 23,000 injuries, more than 500,000 destroyed houses, and 270,000 damaged houses (Lamb & Jones, 2012; NPC, 2015) . People worldwide desire to know the likelihood of earthquakes but neither physical nor statistical models are adequate for predictions and other analysis of seismic pattern (Konsuk & Aktas, 2013; Vere-Jones, Ben-Zion, & Zuniga, 2005) . There is a little evidence of failure of earthquake prediction, but this does not deny the need to look forward and decrease the hazard and loss of life (Nava, Herrera, Frez, & Glowacka, 2005) .
This paper anticipated to deal with the questions 1) What is the frequency-magnitude relationship of earthquake in this region? 2) Every how many years (in average) an earthquake occurs with magnitude ≥ M? 3) What is the probability of an occurrence of at least one earthquake of magnitude ≥ M in the next “t” years? The aim of the earthquake prediction is to aware people about the possible devastating earthquakes timely enough to allow suitable reaction to the calamity and reduce the loss of life and damage from the earthquake occurrence (Vere-Jones et al., 2005; Nava et al., 2005) .
Some researchers believed that the most analysis of seismic hazards is sensitive to inaccuracies in the earthquake catalogue. It is assumed that the long-term earthquake catalogue is not homogeneous and the regular earthquakes, which might include foreshocks and aftershocks of characteristic events, follow Gutenberg-Richter frequency magnitude relationship (Wyss, Shimazaki, & Ito, 1999; Kagan, 1993) . In seismology, the Gutenberg-Richter relation is mainly used to find the association between the frequency and magnitude of the earthquake occurrence because the distributions of earthquakes in any areas of the planet characteristically satisfy this relation (Gutenberg & Richter, 1954; Gutenberg & Richter, 1956) .
The earlier research papers have applied the generalized linear models (GLM), which included Poisson regression, negative-binomial, and gamma regression models, for an earthquake hazard analysis. Konsuk and Aktas (2013) analyzed that the magnitude random variable is distributed as the exponential distribution. The relation between magnitude and frequency is characterized using the Gutenberg Richter function. Shrey and Baker (2011) fitted logistic regression model by maximum likelihood method using generalized linear model for predicting the probability of near fault earthquake ground motion pulses and their period. Turker and Bayrak (2016) estimated an earthquake occurrence probability and the return period in ten regions of Turkey using the Gutenberg Richter model and the Poisson model.
Several studies mentioned that the generalized linear model is used to include a common method for computing parameter estimates, and it also provides significant results for the estimation probabilities of earthquake occurrence and recurrence periods, which are considered as significant parameters of seismic hazard related studies (Nava et al., 2005; Shrey & Baker, 2011; Turker & Bayrak, 2016) . This study is noteworthy on its own from the Statistical and Geoscience perspectives on fitting the models to the earthquake data of Nepal. The important seismic parameters (a and b values) of Gutenberg Richter (GR) relationship and generalized linear models are examined by studying the past earthquake data. It is also intended to estimate the probability of an earthquake occurrence and its return periods of occurring earthquakes in the future “t” years using GR relationship and compared with the Poisson model.
2. Materials and Methods
The earthquake data are obtained from the National Seismological Centre, Department of Mines and Geology, Kathmandu, Nepal, which covers earthquakes from 25th June 1994 through 29th April 2019. It includes epicenter, latitude, longitude, stations, reporting time, and date. The most important factors affecting the seismic hazard in this region are taken into account such as frequency, magnitude, probability of exceedance, and return period of earthquake (Sebastiano, 2012) . The statistical analysis has been accomplished using IBM SPSS 23.0 for Mac OS. In this study, the magnitude values, measured in local magnitude (ML), 4.0 or greater are used for earthquake data. The local magnitude is the logarithm of maximum trace amplitude recorded on a Wood-Anderson seismometer, located 100 km from the epicenter of the earthquake (Sucuogly & Akkar, 2014) . The available data are tabulated for the frequency distribution of magnitude 4 ≤ M ≤ 7.6 and the number of earthquakes for “t” years. The Kolmogorov Smirnov goodness of fit test and the Anderson Darling test is used to check the normality assumption of the data (Gerald, 2012) .
The Kolmogorov Smirnov test statistics is defined by
where, F is the theoretical cumulative distribution of the distribution being tested.
The Anderson Darling test statistics is defined by
Here, F is the cumulative distribution function of the specified distribution and n is the sample size. The Durbin Watson test is used to measure the autocorrelation in residuals from regression analysis.
The Durbin Watson test statistics is calculated using
where, ei are residuals from ordinary least squares regression (Gerald, 2012) .
2.1. Gutenberg-Richter (GR) Relation
Gutenberg and Richter (1954) have suggested an expression for the magnitude and frequency of earthquake events larger than magnitude (M). It states that the logarithm of the frequency is linearly dependent on the magnitude of the earthquake. The relation is generally fitted to the data that are available for any region of the globe. The Gutenberg Richter relation is
where, N is a number of earthquakes having magnitude larger than M during a time period “t”, logN is a logarithm of the number of earthquakes with magnitude M, “a” is a constant that measures the total number of earthquakes at the given source or measure of seismic activity, and “b” is a slope of regression line or measure of the small versus large events.
The annual frequency of exceeding the M event magnitude is computed dividing the number of events N by the “t” years,
Taking logarithm on both sides of Equation (5) we get,
The probability of occurrence of at least one earthquake of magnitude ≥ M in the next “t” years is
The number of years, in an average, an earthquake occurs with magnitude ≥ M is given by
2.2. Generalized Linear Models
In the present study, generalized linear models (GLM) are applied as it basically eliminates the scaling problem compared to conventional regression models. An important characteristic of GLM is that it assumes the observations are independent. The other assumption about the error structure is that there is, a single error term in the model. The normality and constant variance properties are not a compulsion for the error component. GLM is most commonly used to model count data. The selection of measurement scale is a significant feature of model selection; for example, in this study, transformed scale, such as logN and lnN are assumed to be better for additivity of systematic effects (McCullagh & Nelder, 1989) .
Let be the independent response observations with mean respectively. The random element Y has an independent normal distribution with constant variance σ2 and E(Y) = μi. The systematic component: covariates produce a linear predictor . The link between the random and systematic components is .
The generalized linear model is made up of a linear predictor
and two functions 1) a link function that describes how the mean, E(Y) = μi, depends on the linear predictor and 2) a variance function that describes how the variance, Var(Y) depends on the mean, Var(Y) = ϕV(μi), where the dispersion parameter ϕ is a constant (McCullagh & Nelder, 1989; Dobson & Barnett, 2008) .
The dependent variable yi is a count (number of earthquake occurrence), such that . Hence, a rational probability model for count data is frequently the Poisson distribution.
The probability function of a Poisson distribution is given by
where, the parameter μi > 0. The mean and variance of Poisson distribution are equal to the parameter μ. Nevertheless, this statement may not be true and occasionally over dispersion or under dispersion conditions can be observed. When the observed variance is greater than the variance of a theoretical model, over dispersion happens. In the existence of over dispersion, the generalized negative binomial regression model (GNBR) offers an alternative to the generalized Poisson regression model (GPR). If the observed variability is significantly smaller than the predicted variance or under dispersion, Gamma models are more appropriate. GLM allows choosing the suitable model fit on the basis of dispersion parameters and model fit criteria. The procedures of model fitting are 1) model selection 2) parameter estimation and 3) prediction of future values (McCullagh & Nelder, 1989; Kokonendji, 2014) .
2.2.1. Model Selection
The model selection information criteria that are based on likelihood functions and applications to the parametric model based problems are 1) Akaike information criterion (AIC): AIC procedure is generally considered to select the model that minimizes AIC = −2LL + 2d, where LL is the maximized log likelihood of the model given “n” observation, “d” is the dimension of a model. Since the likelihood function’s value is multiplied by −2, ignoring the second component, the model with the minimum AIC is the one with the highest value of the likelihood function. 2) Bayesian information criterion or Schwarz information (BIC): It is also a widespread model selection principle. It selects the model that minimizes . The best model is the one that provides the minimum AIC and BIC (Fabozzi, Focardi, Rachev, Arshanapalli, & Markus, 2014) .
2.2.2. A Goodness of Fit
The goodness of fit of a statistical model is continued to explain how well it fits a set of observed values “y” by a set of fitted values derived from the model. It tests the hypothesis as H0: The model fits, and H1: The model does not fit. It also reviews the inconsistency between observed values and the expected value because a small discrepancy may be acceptable, but not the larger one (McCullagh & Nelder, 1989) . The deviance residual is considered for the generalized measure of discrepancy. The residual sum of squares is the deviance for Normal distribution and is given by . For Poisson regression, the deviance is G2, which is minus twice the log likelihood ratio.
G2 is also called likelihood ratio statistic and is defined as
where, yi is the observed values and is the fitted value. The small value of G2 indicates that the model fits well (Bishop, Fienberg, & Holland, 2007) .
The other significant measure of discrepancy is the generalized Pearson Chi Square statistics, which is given by
where, is the estimated variance function for the distribution concerned. The Pearson Chi square statistics for the Normal distribution is the residual sum of squares, where as for the Poisson distribution it is the Pearson Chi square statistics, and is given by
where, yi is the observed value, and is the expected value under the assumption that null hypothesis is true, i.e. the assumed model is a good one. X2 and G2 are both measure how closely the model fits the observed data. The null hypothesis is rejected if the values of X2 and G2 are large enough. After selecting the model, the unknown parameters have to be estimated. The probability of exceedance in a time period “t”, described by a Poisson distribution, is given by the relationship: . The return period of earthquake is a statistical measurement representing the average recurrence interval over an extensive period of time and is calculated using the relation (Madsen & Thyregod, 2010; Raymond, Montgomery, Vining, & Robinson, 2010; Shroder & Wyss, 2014) .
3. Results and Discussion
The number of occurrence of earthquakes (n) is a count data and the parametric statistics for central tendency, mean = 26 and median = 6 are calculated. It is observed that the most of the values are less than 26; hence, the average value cannot be deliberated as the true representation of the data. It can also be perceived that the data is positively skewed and lacks symmetry; and thus the normality assumption has been severely violated. Therefore, to convert the non-normal data to the normal log transformation of cumulative frequency of earthquakes logN is used.
The hypothesis for normality testing is
H0: The data follow a specified distribution and
H1: The data do not follow a specified distribution.
In order to check the distribution of the transformed variable, first of all Kolmogorov Smirnov test is applied.
Table 1 displays the Kolmogorov Smirnov test statistics for testing specified distribution of data. The p-value is not significant (0.147 > 0.05) and failed to accept H1 for logN, which displayed that normality, exists in the data. In addition, lnN also statistically fitted to the Poisson distribution, the p-values is not significant (0.629 > 0.05). The very severe limitation of the Kolmogorov Smirnov test is that the distribution must be fully specified, i.e. the parameters are known. If location, scale and shape parameters are estimated from the available data, the critical region of this test is no longer valid (Gerald, 2012) . Therefore, the Anderson Darling test is used to observing normality of the data.
The Anderson Darling test is not available in SPSS version 23 and hence it is calculated using Anderson Darling normality test calculator for excel. The result is displayed in Table 2. The p-value = 0.09505 > 0.05 indicates normality.
The Durbin-Watson test is used to determine whether there is evidence of first order autocorrelation in the data and result presented in Table 3. The hypothesis for the Durbin Watson test is H0: There are no first order autocorrelation and H1: The first order correlation exists.
The small value of the D-W score (0.596 < 2) indicates a positive first order autocorrelation, which is assumed to be a common occurrence in this case. Hence, it can be concluded that the observations are linearly independent. The correlation value R = 0.995 specifies that there is a very high degree of association between the magnitude and occurrence of the earthquake. The higher value
Table 1. Kolmogorov Smirnov (K-S) test.
Table 2. Anderson Darling (AD) test.
Table 3. Durbin Watson (D-W) test.
Predictors: (Constant), M. Dependent Variable: logN.
of coefficient of determination (R2 = 0.991) portrayed, the magnitude of earthquake explained 99.1% of the variation in occurrence of earthquake while 0.9% were due to other variables that were not included in the model. According to the results, it is observed that logN and lnN can be considered as dependent variables for Gutenberg-Richter model and generalized Poisson regression model or negative binomial regression model respectively.
The model selection criterion for generalized linear models is illustrated in Table 4. The cumulative frequency of earthquake (N) is divided by the time period (t) and used as a response variable in generalized linear models to select a suitable model. It demonstrates the values of AIC, and BIC for model selection which are reasonably smaller for the GPR model than the normal and GNBR. Hence, the generalized Poisson regression model is considered as the suitable model to fit the data. The significant measures of discrepancy for the Poisson regression model is deviance residual (value/df = 0.170) and generalized Pearson Chi square statistics (value/df = 0.110). These values measure how diligently the model fits the observed data. Furthermore, the generalized Poisson regression model is detected to be the best model to fit the data because 1) it was suitable for count data of earthquake occurrences, 2) model information criterion AIC and BIC are fewer, and 3 deviance and Pearson Chi square statistics are less than one. After selecting the model, the unknown parameters are estimated.
The estimated parameters of the Gutenberg Richter relationship are demonstrated in Table 5. The GR relationship of the earthquakes that had occurred in time period t = 25 years is expressed as logN = 6.532 − 0.887M, where, N is the number of earthquakes ≥ M, logN is the dependent variable, M is the predictor. The model provides the important parameters of the earthquake such as
Table 4. Model selection criterion for GLM.
Table 5. Parameter estimation for Gutenberg Richter model.
a = 6.532, b = −0.887, a' = a − log(bln10) = 6.22, a1= a − log(t) = 5.13, and = a' − log(t) = 4.82.
Table 6 displays the estimated parameters in the generalized Poisson regression model and is given by lnN = 15.06 − 2.04M, where, lnN is the response variable. The other significant parameters of the earthquake are obtained: a = 15.06, b = −2.04, a' = 13.513, a1 = 11.84, and = 10.29.
The relationship between frequency and magnitude of an earthquake ≥ 4 using GR model and GPR model is shown in Figure 1.
The seismic risk expressed in percentage and the return period of the earthquake in years in the Gutenberg Richter model is illustrated in Table 7. The GR relation is logN(M) = 6.532 − 0.887M. The annual frequency of exceeding the M event magnitude is N1(M) = N(M)/t = N(M)/25. Taking logarithm on both sides, logN1(M) = logN(M) − logt = logN(M) − log25 = 6.532 − 0.887M − 1.398 = 5.134 − 0.887*M. For magnitude 7.5, logN1(M ≥ 7.5) = 5.134 − 0.887*7.5 = −1.5185. Now, N1(M ≥ 7.5) = 10(−1.5185) = 0.030305. Hence, the return period for 7.5 magnitude is given by TR(M ≥ 7.5) = 1/N1(M) = 32.99 years. Similarly, the return period for magnitude 6 and 7 are calculated as 1.54 and 11.88 years. The probability of occurrence of at least one earthquake of magnitude ≥ M in the next “t” years, is obtained by the relation, . For illustration, when M = 7.5 and t = 50 years, P(t) = 1 − e(−0.030305*50) = 78%, which is the probability of exceedance in 50 years.
The probability of exceedance expressed in percentage and the return period of an earthquake in years for the Poisson regression model is shown in Table 8. The GPR relation obtained is lnN = 15.06 − 2.04M. The annual frequency of exceeding the M event magnitude for 7.5 ML is calculated as N1(M) = exp(a − bM − lnt) = 0.031. The probability of occurrence of at least one earthquake of magnitude ≥ 7.5 within 50 years is obtained as 79% and the return period is 31.78
Table 6. Parameter estimation for generalized Poisson regression model.
Table 7. Probability of exceedance (%) and return period using GR model.
Table 8. Probability of exceedance (%) and return period using GPR Model.
years. The theoretical return period is the reciprocal of the probability that the event will be exceeded in any one year. There is a statistical statement that on an average, a 10 years event will appear once every ten years and the same process may be true for 100 year event. The theoretical values of return period in Table 8 are slightly greater than the estimated return periods.
Figure 2 demonstrates the probability of earthquake occurrence (%) for different time periods in years using GR and GPR models. In GR model, the
Figure 1. Magnitude (ML)-frequency relation using GR and GPR models.
Figure 2. The probability of exceedance (%) for t years using GR and GPR models.
probability of an earthquake incident of magnitude less than 6 is almost certainly in the next 10 years and more, with the return period 1.54 years. In GPR model, the probability of the earthquake event of magnitude less than 5.5 is almost certainly in the next 5 years and more, with the return period 0.537 years (196 days).
The probability of exceedance using the GR model is found to be less than the results obtained from the GPR model for magnitude higher than 6.0. Likewise, the return periods obtained from both the models are slightly close to each other. The return period values of GPR model are comparatively less than that of the GR model. The probability of exceedance in 10 years with magnitude 7.6 for GR and GPR models is 22% and 23% and the return periods are 40.47 years and 38.99 years respectively. The estimated values depict that the probability of exceedance increases when the time period increases. It can also be noticed that the return period of the earthquake is larger for the higher magnitudes.
The frequency magnitude relationship of the earthquake data of Nepal modelled with the Gutenberg Richter (GR) model is logN= 6.532 − 0.887M and with generalized Poisson regression (GPR) model is lnN = 15.06 − 2.04M. The parameters a and b values for GR and GPR models are (a = 6.532, b = −0.887) and (a =15.06, b = −2.04) respectively. The return periods from GPR model are moderately smaller than that of GR model. In GR model, the return period for 7.5, 7 and 6 magnitudes are 32.99 years, 11.88 years and 1.54 years respectively. In GPR model, the return period for 7.5, 7 and 6 magnitudes are 31.78 years, 11.46 years, and 1.49 years respectively.
This study suggests that the probability of earthquake occurrence produced by both the models is close to each other. In GR model, the probability of earthquake occurrence of at least one earthquake of magnitude ≥ 7.5 in the next 10 years is 26% and the magnitude ≥ 6.5 is 90%. The probability of exceedance of magnitude 6 or lower is 100% in the next 10 years. Similarly, in GPR model, the probability of earthquake occurrence of at least one earthquake of magnitude ≥ 7.5 in the next 10 years is 27% and the magnitude ≥ 6.5 is 91%. The probability of exceedance of magnitude 6 or lower is 100% in the next 10 years.
However, it is very important to understand that the estimated probability of an earthquake occurrence and return period are statistical predicted values, calculated from a set of earthquake data of Nepal. All the parameters required to describe the seismic hazard are not considered in this study. The earthquake catalogue has 25 years of data so the predicted values of return period and the probability of exceedance in 50 years and 100 years cannot be accepted with reasonable confidence. Actually, nobody knows that when and where an earthquake with magnitude ≥ M will occur with probability 1% or more. Nevertheless, the outcome of this study will be helpful for the preparedness planning to reduce the loss of life and property that may happen due to earthquakes because Nepal lies in the high seismic region. Further research can be conducted considering other rational earthquake hazard parameters for different regions that are prone to earthquake occurrence.
The data studied in this paper is the earthquake data from the National Seismological Centre, Department of Mines and Geology, Kathmandu, Nepal, which covers earthquakes from 25th June 1994 through 29th April 2019. It is an open access data available on the website http://seismonepal.gov.np/earthquakes.
 Kokonendji, C. C. (2014). Over-and Under Dispersion Models. In N. Balakrishnan (Ed.), Methods and Applications of Statistics in Clinical Trials: Planning Analysis and Inferential Methods (pp. 506-526). New York: John Wiley & Sons.
 Lamb, M. R., & Jones, B. K. (2012). United States Geological Survey Natural Hazards Response: Geological Survey Fact Sheet 2012-3061. Reston, VA: US Geological Survey.
 Nava, F. A., Herrera, C., Frez, J., & Glowacka, E. (2005). Seismic Hazard Evaluation using Markov Chains: Application to the Japan Area. Pure and Applied Geophysics, 162, 1347-1366.
 Shrey, S. K., & Baker, J. W. (2011). Regression Models for Predicting the Probability of Near-fault Earthquake Ground Motion Pulses, and Their Period. In Proceedings of 11th International Conference on Applications of Statistics and Probability in Civil Engineering (8 p.). Zurich, Switzerland.
 Turker, T., & Bayrak, Y. (2016). A Poisson Method Application to the Assessment of the Earthquake Hazard in the North Anatolian Fault Zone. Challenge Journal of Structural Mechanics, 2, 109-121.
 USGS (2016). M 7.8 Nepal Earthquake 2015—A Small Push to Mt. Everest. Reston, VA: United States Geological Survey.
 Wyss, M., Shimazaki, K., & Ito, A. (1999). Seismicity Patterns, Their Statistical Significance, and Physical Meaning. Pure and Applied Geophysics, 155, 203-205.