Bulletin of Social Informatics Theory and Application  ISSN 2614-0047 

Vol. 1, No. 1, March 2017, pp. 34-40  34 

https://doi.org/10.31763/businta.v1i1.23         

Seismic analysis using maximum likelihood of gutenberg-

richter 

Arum Handini Primandari 1,*, Khusnul Khotimah 2 

Department of Statistics, Islamic University of Indonesia, Indonesia 
1 primandari.arum@uii.ac.id *; 2 khusnul.khotimah@student.uii.ac.id 

* corresponding author 

 

1. Introduction 

An earthquake is natural disaster which occur as result of seismic wave released by earth crust [1]. 
It happens mostly triggered by the rupture of geological faults and volcano activities which cause mild 
to severe tremor. Other causes such as landslide, mine blast, and nuclear test may result mild tremor. 

To measure the velocity of seismic wave, an instrument called seismograph is used [2]. This 
instrument is installed under the earth’s surface and is connected to seismogram (a device to record 
ground motion). Seismogram works by recording the ground motion which is visualized through 
vertical waves [3]. Unlike seismograph, seismogram is placed in an observation station. To define the 
size of the magnitude of earthquake, a scale named Richter scale is developed. This scale is determined 
from the base-10 logarithm of the amplitude of seismic waves [4]. 

The latitude of D.I. Yogyakarta, Indonesia is -7.797068, and the longitude is 110.370529. D.I. 
Yogyakarta is located in Java which lies on convergent plateau zone where the Indo-Australian Plate 
is sub ducted under the Eurasian Plate. The offshore plate movement will suppress the onshore plates 
which is located near D.I. Yogyakarta. Therefore, D.I. Yogyakarta has highly possibility to encounter 
earthquakes. 

There are many considerable researches conducted to predict the place and time the earthquake 
likely happens. However, until now there is no definite method that can precisely predict the 
occurrence of earthquake. The least the researcher can do is to conjecture the potential tremor to occur. 
Thus, we can make any effort to reduce the impact of the catastrophe.  

Aki and Utsu proposed estimator obtained from maximum likelihood method [5]. This estimator 
is called b-value which used to determine seismicity rate. Since it is an estimator, there are the 
uncertainty factor (bias) [6]–[8]. The probability of earthquake may be obtained since the seismicity 
rate had been calculated [6], [9]. 

A R T I C L E  I N F O   A B S T R A C T  

 

Article history 

Received December 21, 2016 

Revised January 14, 2017 

Accepted February 4, 2017 

 An earthquake is one of catastrophe which often claim numerous lives and cause 
great damage to infrastructure. Multiple studies from various field have been 
conducted in order to make a precise prediction of earthquake occurrence, such 
as recognizing the natural phenomena symptoms leading to the shaking and 
ground rupture. However, up till now there is no definite method that can 
predict the time and place in which earthquake will occur. By assuming that the 
number of earthquake follow Gutenberg-Richter law, we work b-value derived 
using Maximum Likelihood Method to calculate the probability of earthquake 
happen in the next few years. The southern sea of D.I. Yogyakarta was divided 
into four areas to simplify the analysis. As the result, in the next five years the 
first and second area have high enough probability (>0.3) to undergo more than 
6.0-magnitude earthquake.  

 
This is an open access article under the CC–BY-SA license. 

    

 
Keywords 

Earthquake 

Gutenberg-richter 

Maximum likelihood 

Probability prediction 

 

 

http://creativecommons.org/licenses/by-sa/4.0/
http://creativecommons.org/licenses/by-sa/4.0/


ISSN 2614-0047 Bulletin of Social Informatics Theory and Application 35 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

The Gutenberg-Richter  (GR)  distribution  is  widely  used  to  show  statistical relation, 

𝑙𝑜𝑔
10

 𝑁(𝑀) = 𝑎 − 𝑏𝑀  

Where N(M) is the number of earthquake with magnitude equal to or greater than M, a and b are 
constants which show level and the character of seismicity in the region of concern. Using Maximum 
Likelihood Method, the estimator called b-value is derived. Since b-value is obtained, we also can get 
the value of a. 

Seismic activity of D.I. Yogyakarta is able to be quantified by parameters that are seismicity  index,  
probability  of  tremor  occurrence,  and  return  period  of  the highest tremor within the area. The 
probability describes the risk of earthquake including only seismic information, without others factor 
such as geology structure, population density, quality of infrastructure, etc. 

The hazard risk of tremor is useful for construction planning. Therefore, we can build the 
earthquake resistant buildings.  Earthquakes are classified into some categories from minor to great. 
This research calculated the probability of tremor with 4.0 magnitude or more which potentially 
destructive. 

The outline of this paper was arranged as follow, first introduction described the background of 
research. Second, problem formulation states what this research would do. Method, the third 
outline, explains the source of data and the tools which is used to analyze the data. The seismic 
analysis of D.I. Yogyakarta was provided in fourth section. The conclusion and open problem, 
successively was written in fifth and sixth section. 

2. Problem Formulation 

To describe the seismicity condition in D.I. Yogyakarta from 2011 – 2015, we need a series 
of measurable parameters. This study is to quantify these parameters that are seismicity index, 
earthquake probability, and return period. The b-value and a which are obtained from Maximum 
Likelihood Method, used to draw the seismicity index. This index represents the number of 
earthquake occurrence with magnitude equal to or more than M0. 

Since seismicity index is given, the problem is how to calculate the probability for an area to 
encounter a number of earthquakes with magnitude equal to or more than M0  and the return 
period. Return period indicates the repetitive interval (in years) of earthquake to occur. In this 
research we analyzed the tremor which have magnitude equal to 4.0 (on the Richter scale) or more. 

3. Method 

3.1. Data 

Earthquake data from 2011 – 2015 is obtained from Subsection Data and Information of 
Meteorology Climatology and Geophysics Agency (BMKG) of D.I. Yogyakarta. The epicenter of 
earthquake located in southern offshore of D.I. Yogyakarta divided into six areas. 

1. Area 1       :  -8.26° – -9,31° latitude and 110,01° – 110,47° longitude 

2. Area 2       :  -8,26° – -9,31° latitude and 110,47° – 110,93° longitude 

3. Area 3       :  -9,37° – -10,37° latitude and 110,01° – 110,47° longitude 

4. Area 4       :  -9,37° – -10,37° latitude and 110,47° – 110,93° longitude 

5. Area 5       : -10,37° – -11,42° latitude and 110,47° – 110,93° longitude 

6. Area 6       : -10,37° – -11,42° latitude and 110,01° – 110,47° longitude 



36 Bulletin of Social Informatics Theory and Application   ISSN 2614-0047 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

 

Fig. 1. The division of earthquake zone 

We did not include area 5 and 6 into analysis because it has few tremors. There was even no 
tremor with magnitude more than 5.0 Richter. Moreover, area 5 and 6 is far from D.I. Yogyakarta 
onshore. 

3.2. Gutenberg-Richter Law 

In seismology, the Gutenberg-Richer (G-R) law holds [10], 

𝑙𝑜𝑔
10

𝑁 (𝑀) = 𝑎 − 𝑏 (𝑀 − 𝑀0; 𝑀 ≥ 𝑀0)  

where N ( M ) is the number of earthquake in particular magnitude M, M0 is the minimum 
magnitude above which all earthquake within a certain region are recorded. Parameter b, called b-
value, and parameter a are constant related to seismicity rate. 

1) Maximum Likelihood Estimator (MLE) 

Equation (1) implies that magnitude distributed exponentially [10]: 

𝑝(𝑀) = 𝛽𝑒 −𝛽(𝑀−𝑀0); 𝑀 ≥ 𝑀0  

where, β = bln(10), p(M) is the probability density function of M. MLE of b-value was derived by: 

a) The likelihood function defines as follows: 

𝐿(𝑀) = ∏ 𝛽𝑒 −𝛽(𝑀𝑖−𝑀0)𝑛𝑖−1 

           = 𝛽𝑛𝑒 −𝛽 ∑ (𝑀𝑖−𝑀0)
𝑛
𝑖=1   

b) The natural logarithm of likelihood function is given below: 

ln 𝐿 (𝑀) = 𝑛 𝛽 −  𝛽 ∑ (𝑀𝑖 − 𝑀0)
𝑛
𝑖=1    







ISSN 2614-0047 Bulletin of Social Informatics Theory and Application 37 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

c) Define the parameter β  which optimized the log logarithm likelihood function by taking the derivative 
and setting it equal to 0. 

𝜕 ln 𝐿 (𝑀)

𝜕𝛽
=  

𝑛

𝛽
−  ∑ (𝑀𝑖 −  𝑀0)  = 0

𝑛
𝑖=1   

𝛽 =  
1

�̅�− 𝑀0
   

By substituting β = bln(10) to equation (7), we obtain the MLE of β : 

�̂� =  
1

�̅�−𝑀0
  

Thus we have MLE of parameter β which called b-value. Instead of equation (8), this research used 
the modified b-value below. 

�̂� =  
log 𝑒

�̅�−(𝑀0− ∆𝑀/2 )
  

In practice, magnitude is rounded into ∆M, usually ∆M = 0.1 . The uncertainty of b-value defines 
as follow [11]: 

�̂�𝑏
2 =  

�̂�

√𝑁
  

where N is the number of earthquakes occurrence. Shi and Bolt provide formula to estimate the 
error of b-value that is [11]: 

�̂�𝑏 = 2.30�̂�
2√

∑ (𝑀𝑖− �̅�)
2𝑁

𝑖=1

𝑁 (𝑁−1)
  

Compared equation (10), equation (11) have reliable estimation for b-value error. Parameter a and 
b in G-R distribution are constant. The value of a is formulated as follow [7]: 

�̂� = 𝑙𝑜𝑔10  (
𝑁(𝑀)

𝑇
) + 𝑏𝑀0 ; 𝑀 ≥ 𝑀0   

where T is time interval in which the observation recorded. Holding the estimation of parameter a 
and b, the seismicity index (or rate) of an area with magnitude M ≥ M0 can be written bellow:  

𝑁1(𝑀) = 10
�̂�−�̂�𝑀  

The formula to determine return period is given bellow: 

𝜃 =  
1

𝑁1(𝑀)
  

2) Probability 

The probability of earthquake occurrence with magnitude equal to or more than M within T period, 
usually called cumulative distribution function, given as follow: 

𝑃(𝑀, 𝑇) = 1 − 𝑒 −𝑁1(𝑀)𝑇  









38 Bulletin of Social Informatics Theory and Application   ISSN 2614-0047 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

4. Seismicity Analysis of D.I. Yogyakarta 

We define M0 = 4, since the data recorded was all the earthquakes with magnitude more than or 
equal to 4.0 occurred in different depth. The summary of the data is explained in Table 1. 

Table 1.  Summary of earthquake data of D.I. Yogyakarta offshore 2011-2015 

No Area Mmin Mmax �̅� N(Mmin) 
1 Area 1 4.0 5.8 4.435 31 

2 Area 2 4.0 5.5 4.484 19 

3 Area 3 4.0 4.5 4.300 6 

4 Area 4 4.1 4.5 4.217 6 

 

In area 1, earthquake with 4.0-magnitude until 5.8-magnitude happened in 10 km until 50 km 
depth, while the 5.8-magnitude earthquake happened in 50 km depth. The distribution of earthquake 
based on its depth is showed in Fig. 1. Meanwhile, others areas’ description was provided in appendix. 

 

Fig. 2. The distribution of earthquakes based on its depth in area 1 

Earthquake magnitudes are classified into some categories that are light (4.0-4.9), moderate (5.0-
5.9), strong (6.0-6.9), major (7.0-7.9), and great (>8.0). Since the data recorded have more than or 
equal to 4.0-magnitude, we calculate the seismicity index for 4.1-magnitude. 

Table 2.  Seismicity index for M  ≥ 4.1 

No. Area M0 b-value a-value N1  M

1 Area 1 4 0.895 4.374 5.045 

2 Area 2 4 0.813 3.833 3.151 

3 Area 3 4 1.241 5.043 0.902 

4 Area 4 4 1.627 6.585 0.825 

 

According to Table 2, area 1 has seismicity index at 5.045 which means every year area 1 
encounters about 5 earthquakes with magnitude more than or equal to 4.1. In area 3 and 4, tremor 
happen about once for every year. Meanwhile, in area 2, tremors happen four times within a year. 
Estimation of b-value error in every area is shown at Table 3. 

 

 

 

 



ISSN 2614-0047 Bulletin of Social Informatics Theory and Application 39 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

Table 3.  Estimation of b-value error 

No. Area M0 b-value ˆ 

1 Area 1 4 0.895 0.151 

2 Area 2 4 0.813 0.147 

3 Area 3 4 1.241 0.259 

4 Area 4 4 1.627 0.290 

 

The probability having tremor with magnitude more than or equal to 4.1 in the next of five years 
is showed in Table 4. 

Table 4.  Probability of earthquake with magnitude more than or equal to 4.1 

No. Area N1(M) P(  M ≥  4.1)

1 Area 1 5.045 1.000 

2 Area 2 3.151 1.000 

3 Area 3 0.902 0.989 

4 Area 4 0.825 0.984 

 

In five years ahead, every area is almost certain for having earthquake with given magnitude. D.I. 
Yogyakarta encountered earthquake with 5.9-magnitude for about 57 seconds. This natural disaster 
caused a lots of deaths, about 5700, and injuries. The epicenter was located near area 1 and 2 in 33 
km depth. Thus, we took the strong earthquake with minimum magnitude at 6.0 and calculate both 
seismicity index and the probability for experiencing this earthquake in five years ahead. 

Table 5.  Seismicity index and probability for M ≥ 6.0 

No. Area N1  M  P  M  6.0

1 Area 1 0.100 0.395 

2 Area 2 0.090 0.362 

3 Area 3 0.004 0.020 

4 Area 4 0.001 0.003 

 

The earthquake with magnitude more than 6.0-magnitude may cause a lot of damage. Area 1 has 
the highest probability among the four area for going through this earthquake. 

5. Conclusion 

Employing Maximum Likelihood Method, we obtain MLE called b-value which is used to 
determine seismicity index. The b-value of earthquake with 4.1 magnitude or more of the four area 
ranged from 0.895 – 1.627. While seismicity index obtained is about 0.825 – 5.045, which means 
there are about 1 – 5 times earthquakes happen within a year. Since seismicity index provided, the 
probability was calculated. The probability for having more than or equal 4.1-magnitude tremor in 
five years ahead is almost certain for all areas.  

The minimum magnitude for hazardous earthquake was given at 6.0-magnitude. The highest 
probability for encountering this earthquake is 0.395 which occurred in area 1. Meanwhile, the lowest 
probability is 0.003 which occurred in area 3. Whereas, area 1 is closer D.I. Yogyakarta onshore than 
area 3. 

6. Open Problem 

The formula for determining seismicity rate should consider the depth where the earthquake 
happen. The shallower the epicenter of earthquake, the more hazardous to the population. 



40 Bulletin of Social Informatics Theory and Application   ISSN 2614-0047 
 Vol. 1, No. 1, March 2017, pp. 34-40 

 Seismic analysis using maximum likelihood of gutenberg-richter  

References 

[1] K. D. Marano, D. J. Wald, and T. I. Allen, “Global earthquake casualties due to secondary effects: a 

quantitative analysis for improving rapid loss analyses,” Nat. hazards, vol. 52, no. 2, pp. 319–328, 2010. 

[2] I. Koulakov, S. El Khrepy, N. Al-Arifi, I. Sychev, and P. Kuznetsov, “Evidence of magma activation 

beneath the Harrat Lunayyir basaltic field (Saudi Arabia) from attenuation tomography,” Solid Earth, vol. 

5, no. 2, pp. 873–882, 2014. 

[3] Incorporated Research Institutions for Seismology, “How Does a Seismometer Work?,” 2007. [Online]. 

Available: https://www.iris.edu/hq/inclass/fact-sheet/how_does_a_seismometer_work. [Accessed: 17-

May-2016]. 

[4] D. M. Boore, “The Richter scale: its development and use for determining earthquake source parameters,” 

Tectonophysics, vol. 166, no. 1–3, pp. 1–14, 1989. 

[5] A. Kijko and A. Smit, “Extension of the Aki‐ Utsu b‐ value estimator for incomplete catalogs,” Bull. 
Seismol. Soc. Am., vol. 102, no. 3, pp. 1283–1287, 2012. 

[6] A. Barth, F. Wenzel, and C. Langenbruch, “Probability of earthquake occurrence and magnitude 

estimation in the post shut-in phase of geothermal projects,” J. Seismol., vol. 17, no. 1, pp. 5–11, 2013. 

[7] M. Lindman, “Physics of Aftershocks in the South Iceland Seismic Zone: Insights into the earthquake 

process from statistics and numerical modelling of aftershock sequences,” Uppsala University, 2009. 

[8] C. Godano, E. Lippiello, and L. de Arcangelis, “Variability of the b value in the Gutenberg–Richter 

distribution,” Geophys. J. Int., vol. 199, no. 3, pp. 1765–1771, 2014. 

[9] A. Budiman, R. Nandia, and M. T. Gunawan, “Analisis Periode Ulang Dan Aktivitas Kegempaan Pada 

Daerah Sumatera Barat Dan Sekitarnya,” J. Ilmu Fis., vol. 3, no. 2, pp. 55–61, 2011. 

[10] A. De Santis, G. Cianchini, P. Favali, L. Beranzoli, and E. Boschi, “The Gutenberg–Richter law and 

entropy of earthquakes: two case studies in Central Italy,” Bull. Seismol. Soc. Am., vol. 101, no. 3, pp. 

1386–1395, 2011. 

[11] W. Marzocchi and L. Sandri, “A Review and New Insights on The Estimation of The b-value and Its 

Uncertainty." Annals Of Geophysic 46 (6): 1271-1281.,” Ann. Geophys., vol. 46, no. 6, pp. 1271–1282, 

2003.