Adaptive estimation of the Gutenberg–Richter b value using a state space model and particle filtering

Earthquakes follow an exponential distribution referred to as the Gutenberg–Richter law, which is characterized by the b value that represents a ratio of the number of large earthquakes to that of small earthquakes. Spatial and temporal variation in the b value is important for assessing the probability of a larger earthquake. Conventionally, the b value is obtained by a maximum-likelihood estimation based on past earthquakes with a certain sample size. To properly assess the occurrence of earthquakes and understand their dynamics, determining this parameter with a statistically optimal method is important. Here, we discuss a method that uses a state space model and a particle filter, as a framework for time-series data, to estimate temporal variation in the b value. We then compared our output with that of a conventional method using data of earthquakes that occurred in Tohoku and Kumamoto regions in Japan. Our results indicate that the proposed method has the advantage of estimating temporal variation of the b value and forecasting magnitude. Moreover, our research suggests no heightened probability of a large earthquake in the Tohoku region, in contrast to previous studies. Simultaneously, there is the potential of a large earthquake in the Kumamoto region, emphasizing the need for enhanced monitoring.


3/17
Supplementary Figure 4.The procedure for evaluating a method for predicting the distribution of magnitude is as follows: (a) infer the b value one period ahead and calculate the q-th quantile of the predictive magnitude.(b) Determine where an actual magnitude exceeds the q-th quantile or not and count the number of events exceeding the threshold.N exec.(n) represents the cumulative number of events until the n-th event.(c) Calculate the maximum distance between the two lines N exec.(n) and the expected value nq.This maximum distance is denoted as Loss.(d) Add a point at (q, Loss).Moving average kernels for one window.In the conventional moving average-based method, as illustrated, weights are applied to the values between sliding windows to compute the average.The Simple Moving Average (SMA) uses a constant weight within the window, the Exponential Moving Average (EMA) employs weights that decrease exponentially for past values, and the Weighted Moving Average (WMA) utilizes weights that decrease linearly.Supplementary Figure 8.The normalized number of times the cumulative count of events has surpassed the q-th percentile of the predictive distribution of magnitude by the time of the n-th earthquake, denoted as N exc.(n) in the Method (see Evaluation Method for Magnitude Prediction Distribution section for details).The panel displays, horizontally, the q-quantiles of the predicted distribution of magnitudes, and vertically, the results of our proposed method (SSM+PF) alongside three moving averages (SMA, EMA, WMA).The blue line represents the cumulative values of actual N exc.(n), while the orange line indicates theoretical values.Results in all panels are of the Tohoku dataset.

1 .
The same kind of plots as Fig.2in the main text.Black line and filled area indicate median and 50% interval of the posterior distribution of the b value estimated by a particle filter (model 2), respectively.

2 .
The same kind of plots as Fig.2in the main text.Black line and filled area indicate median and 50% interval of the posterior distribution of the b value estimated by a particle filter (model 3), respectively.

3 .
The same kind of plots as Fig.2in the main text.Black line and filled area indicate median and 50% interval of the posterior distribution of the b value estimated by a particle filter (model 4), respectively.

Supplementary Figure 5 . 6 .
Weights for applying the weighted average method.(a) Tohoku and (b) Kumamoto cases.This chart represents the weights for applying a weighted average to the results of the traditional moving average method.The stacked graph displays the weights for the results based on the settings of SMA-50, SMA-75, SMA-100, SMA-125, SMA-150, SMA-175, SMA-200, EMA-50, EMA-75, EMA-100, EMA-125, EMA-150, EMA-175, EMA-200, WMA-50, WMA-75, WMA-100, WMA-125, WMA-150, WMA-175 and WMA-200 respectively, from the bottom up.SMA, EMA and WMA represent simple moving average, exponential moving average and weighted moving average, respectively.The numbers, e.g.50, 75, etc., represent the window width of moving average.The early period in the Kumamoto (b) region is a period in which the number of earthquakes was small and the b value tended to increase (Fig.2b), and the weight of the long moving average window was large.Results of b value estimates for artificial data using Model 4. As described in the Method, although Model 4 is the most complex among Models 1-4, it is able to estimate b values most adaptively without optimizing any parameter.Therefore, we used Model 4 to demonstrate its behavior with test data, showcasing its performance and utility in a synthetic case.The median and 95% confidence intervals of estimated b values using particle filtering are represented by black lines and blue regions, respectively.The results of conventional moving average-based methods SMA-200, EMA-200, WMA-200 and WAC-opt are indicated with a blue, orange, green and red line, respectively.The true b value of artificial data is displayed by a blue dashed line (0 < t ≤ 400, b = 0.8; 401 < t ≤ 700, b = 0.5; 701 < t ≤ 1000, b = 0.7).The stem plot shows magnitude generated by GR law with the true b value as a function of time.As the estimates are updated based on the data, the early estimates deviate significantly from the true values according to the initial distribution setting.

Supplementary Figure 9 .
The same kind of plots as Supplementary FigureS8.Results in all panels are of the Kumamoto dataset.