## 3 Results and Discussion  
A descriptive summary of the data by each spatial block is shown in Table 1. In this table, we can see that the magnitude of the maxima for each block is varied and heterogeneous. In the area of Block 4, whose central point is in the parallel located at 15\({}^{\circ}\) N and in the meridian located at 94.5\({}^{\circ}\) W,  
Figure 3: Random censoring. X, event; M, maxima; U, censored observation. The blue dotted line shows the limits of the blocks. the maximum intensity recorded was \(M_{w}8.2\) regarding to the magnitude of the moment, in contrast to the maximum intensity measured in Block 22, where the maxima was \(M_{w}4.6\). In six of the twenty-two blocks, the maximum magnitude exceeded \(M_{w}7.0\). An important characteristic that we observed in all the blocks is that there is a marked difference between the third and fourth quartiles, which indicates that a distribution of heavy tails should be appropriate for the study of this type of data. Clearly, tectonic and physical conditions in southern Mexico vary from region to region, causing seismic patterns that are reflected in extreme events of earthquakes in surrounding areas. Furthermore, these trends may also vary from one time period to another. We verified this behavior again in the statistics of means and quantiles presented in Table 1, which show that the magnitudes of the maxima studied tend to vary from one place to another. This fact justifies the modeling of the location parameter with respect to time and space. The covariates used to estimate spatiotemporal trends were latitude, longitude and time. A hierarchical Bayesian model was fitted to analyze the maxima of earthquakes on the Pacific coast of southern Mexico. At the first level, we modeled the intensities of maxima earthquakes using the generalized extreme value distribution with random censorship; at the second level, we used multivariate smoothing spline functions to model nonstationary spatiotemporal extremes; and, at the third level, we assumed a priori distribution functions for the parameters of the model. Table 2 shows the estimates of the fitted model. In other words, this formulation has the advantage of analyzing maxima through hierarchical models, or stages. The first stage represents the previous information; the second stage represents the trends in the maxima; and the third stage represents the distribution of the data. We constructed the joint probabilities distribution of the parameters and the data, also known as posterior density, using the conditional distributions, mentioned above as stages, formally represented in Equation (2). The procedure used in this work can be easily replicated by: (1) constructing the likelihood function given by Section 2.2.2; and (2) running the Metropolis-Hastings algorithm with the acceptance probabilities given by Equation (4). We fitted our Bayesian hierarchical model with the conditions expressed in Equation (1), setting \(p_{2}=80\). The statistical analysis was performed using the \(R\) 3.3.1 software. The samples were generated using a  
\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline
**Block ID** & **Long (W)** & **Lat (N)** & **Min. ** & **1 st Qu. ** & **Median** & **Mean** & **3 rd Qu. ** & **Max. ** \\ \hline
1 & \(-91.5^{\circ}\) & \(14^{\circ}\) & 6.1 & 6.2 & 6.8 & 6.7 & 7.3 & 7.3 \\
2 & \(-91.5^{\circ}\) & \(15^{\circ}\) & 6.9 & 6.9 & 7 & 7 & 7 &  
Metropolis-Hastings algorithm, with a Gaussian random walk. We generated 300,000 samples using a burn-in period of 50,000. A look at the logarithm of the a posteriori likelihood indicates that the algorithm successfully converged. The spatial smoothing for the years 1995, 2000, 2005, 2010, 2015 and 2018 for the location functions is shown in Figure 4a-f, respectively. Through this, we observed that the model captured the spatiotemporal patterns present in the data. Similarly, we observed variations of the location parameter in the spatial plane.