## 4 Extreme Value Probability Density Maps  
Let \(U\) be any value from the matrix (7), the estimates of which were obtained in a moving time window. For each grid node \((i,j)\) and for each time window labeled \(t\) at the right end of the window, we find the 10 closest seismic working stations, which gives 10 values of \(U\). Let us take their median value \(U_{ij}^{(t)}\). The values \(U_{ij}^{(t)}\) form an \"elementary\" map corresponding to a time window of 365 days. Consider the values of the parameter as a function of two-dimensional vectors \(z_{ij}=(x_{i},y_{j})\) of longitudes and latitudes of nodes \((i,j)\) in an explicit form: \(U_{ij}^{(t)}\equiv U^{(t)}(z_{ij})\). For each \"elementary map\" with a discrete time index \(t\), we find the coordinates of the nodes \(z_{mn}^{(t)}=(x_{m}^{(t)},y_{n}^{(t)})\) at which the extreme value \(U\) is reached with respect to all other nodes of the regular grid. As noted above, the choice of the minimum for entropies and the maximum for spectral slopes is due to the considerations that the regions of active tremor should be characterized by the values of statistics that reflect the maximum deviation from the properties of white noise. For passive tremor, the opposite is true: maxima for entropies and minima for spectral slopes. A cloud of 2D vectors \(z_{mn}^{(t)}\) considered within a certain time interval \(t\in[t_{0},t_{1}]\) composes some random set. Let us estimate two-dimensional probability density function for each node of the regular grid. For this, the Parzen-Rosenblatt estimate with the Gaussian kernel function [23] will be applied:  
\[p(z_{ij}|t_{0},t_{1})=\frac{1}{2\pi h^{2}(t_{1}-t_{0}+1)}\sum_{t=t_{0}}^{t_{1 }}\sum_{mn}\exp\left(-\frac{\left|z_{ij}-z_{mn}^{(t)}\right|^{2}}{2h^{2}}\right) \tag{8}\]Here, \(h\) is the averaging radius, integer indices \(t_{0},t_{1}\) enumerate the maps of each time window. The \((t_{1}-t_{0}+1)\) equals the number of time windows in the considered time interval. The width of the smoothing band \(h=0.25^{\circ}\) was used, which is approximately equal to 28 km. In our analysis, time windows of two lengths are used: a \"short\" window of 365 days in length, which is used to obtain estimates of tremor parameters, and a \"long\" window of 1095 days (3 years) in length, in which a sequence of estimates of the distribution densities of extreme values of statistics is the obtained tremor. Estimates of characteristics from matrix (7) are calculated using \"short\" time windows 365 days long. Next, the probability densities of extreme values are calculated from the values that fall within the current \"long\" time window of 1095 days, which are also shifted by 7 days. Each \"long\" window includes 105 \"short\" time windows. As a result of such assessments, 24 maps of the distribution density of extreme values of tremor parameters were constructed, presented in Figure 3, on which 12 maps for active tremor are in the upper half, and 12 maps for passive tremor are in the lower half. Figure 3: Average probability densities of distribution of extreme values of 4 statistics of “active” and “passive” tremor of the Earth’s surface for 3 components of the displacement vector. Due to the fact that each property value is estimated over a 365-day window offset by 7 days, the average probability density estimate is obtained by averaging 548 probability density maps, each derived from 105 tremor property value maps. Purely visually, a correlation of probability density maps of the distribution of extreme values of tremor statistics is noticeable. The principal component method [29] is used to extract the most common spatial characteristics. In each \"long\" 3-year time window, a correlation matrix of the size 12 \(\times\) 12 between the probability densities of the distribution of extreme values of tremor statistics is calculated.