Forecasting-Test-Methods

Earthquake forecasting performance evaluation serves as a critical metric for quantitatively assessing and comparing the efficacy and reliability of predictive models.

The predictability of earthquakes and the issue of how to conduct earthquake forecasting have attracted intense attention. In seismology, commonly used models such as ETAS, STEP, and EEPSA have been compared and tested within RELM and CSEP projects, while an increasing number of new forecasting models and evaluation methods have recently been developed.

Example evaluation result

[Now] 2026-08-05 23:19:20

# Final forecast evaluation
# Number-Test:
#   Nobs    = 311
#   delta1  = 0.26
#   delta2  = 0.76
#   status  = consistent
# ---------------------------------------------------------
# Magnitude-Test:
#   dobs    = 0.826
#   gammaM  = 0.61
#   pUpper  = 0.39
#   status  = consistent
# ---------------------------------------------------------
# Spatial-Test:
#   Sobs    = -4.553
#   gammaS  = 0.19
#   status  = consistent
# =========================================================

Results plot

---

Catalog-Based Earthquake Forecast Evaluation

This project demonstrates a simulation-based framework for evaluating catalog-based earthquake forecasts. The forecast is represented by an ensemble of stochastic earthquake catalogs, and its predictive performance is assessed by comparing the observed catalog with the forecast ensemble.

Three complementary consistency tests are implemented:

1. Number-Test — evaluates the total number of earthquakes;
2. Magnitude-Test — evaluates the magnitude-frequency distribution;
3. Spatial-Test — evaluates the spatial distribution of earthquakes.

Together, these tests isolate three important components of an earthquake forecast:

Test	Forecast component	Main question
Number-Test	Rate	Is the observed number of earthquakes consistent with the forecast?
Magnitude-Test	Magnitude distribution	Are the observed magnitudes consistent with the forecast magnitude-frequency distribution?
Spatial-Test	Spatial distribution	Are the observed earthquake locations consistent with the forecast spatial distribution?

The demonstration script generates >500 stochastic forecast catalogs and one synthetic observed catalog. The three tests are then applied independently to diagnose possible disagreements between the forecast and observation.

---

1. Forecast and observed catalogs

1.1 Forecast catalog ensemble

Let

$$ \Lambda_j,\qquad j=1,\ldots,J $$

denote the simulated forecast catalogs, where $J$ is the number of catalogs. In the demonstration,

The number of earthquakes in each forecast catalog simply follows a Poisson distribution:

where the expected number of earthquakes is

$$ \bar{N}=300. $$

Each simulated earthquake has:

• a spatial cell;
• a magnitude-bin index;
• an associated magnitude.

The forecast ensemble is therefore represented by:

• ForecastCounts(j) is the number of earthquakes in catalog $j$;
• ForecastSpatialCounts(k,j) is the number of earthquakes in spatial cell $k$ of catalog $j$;
• ForecastMag{j} contains the magnitudes in catalog $j$.

---

1.2 Spatial forecast

The testing region is divided into a regular $14\times14$ spatial grid:

Nx=14;
Ny=14;
Xv=linspace(0,100,Nx);
Yv=linspace(0,100,Ny);

The spatial forecast contains three Gaussian-like earthquake clusters. For a cluster centered at $(X_c,Y_c)$ with spatial scale $S$, the unnormalized spatial field is

$$ f(x,y)=\exp\left[ -\frac{(x-X_c)^2}{2S^2} -\frac{(y-Y_c)^2}{2S^2} \right]. $$

The complete forecast field is constructed as

$$ F(x,y)=1.2f_1(x,y)+ 0.4f_2(x,y)+ 0.6f_3(x,y)+ \varepsilon, $$

where $\varepsilon=10^{-4}$ is a small background value used to prevent zero-probability cells.

The spatial probability in cell $k$ is

$$ p_k=\frac{F_k}{\sum_{l=1}^{K}F_l}, $$

where $K=N_xN_y$ is the number of spatial cells.

The forecast map shows the largest earthquake probability near the lower-left cluster, with weaker clusters near the center and upper-right portion of the testing region.

---

1.3 Magnitude forecast

Forecast magnitudes are generated from a truncated Gutenberg–Richter distribution:

$$ P(M\geq m) \propto 10^{-b(m-M_{\min})}. $$

For the magnitude interval $[m_k,m_{k+1})$, the forecast probability is

$$ p_k^{(m)}=10^{-b(m_k-M_{\min})}-10^{-b(m_{k+1}-M_{\min})}. $$

The probabilities are normalized so that

$$ \sum_k p_k^{(m)}=1. $$

---

2. Number-Test

2.1 Purpose

The Number-Test evaluates whether the total number of observed earthquakes is consistent with the number distribution generated by the forecast catalogs.

For a catalog-based forecast, the test distribution is obtained directly from the simulated catalogs:

$$ N_j=|\Lambda_j|, \qquad j=1,\ldots,J. $$

The observed statistic is simply

$$ N_{\mathrm{obs}}=|\Omega|, $$

where $\Omega$ denotes the observed catalog.

The Number-Test isolates the rate component because it ignores the locations and magnitudes of individual earthquakes.

In grid-based forecast evaluation, this test is commonly called the N-Test.
When an ensemble of complete earthquake catalogs is evaluated directly, the name Number-Test is more explicit. The figure title may still use the shorter label “N-Test”.

---

2.2 Tail probabilities

Two empirical tail probabilities are calculated:

$$ \delta_1=P(N_j\geq N_{\mathrm{obs}}), $$

and

$$ \delta_2=P(N_j\leq N_{\mathrm{obs}}). $$

They are estimated from the forecast ensemble as

$$ \widehat{\delta}_1=\frac{ \sum_{j=1}^{J} I(N_j\geq N_{\mathrm{obs}}) }{J}, $$

$$ \widehat{\delta}_2=\frac{ \sum_{j=1}^{J} I(N_j\leq N_{\mathrm{obs}}) }{J}, $$

where $I(\cdot)$ is an indicator function.

Because both probabilities include equality, their sum may be slightly larger than one when some simulations have exactly $N_{\mathrm{obs}}$ events.

---

2.3 Interpretation

For a two-sided test with significance level

$$ \alpha=0.05, $$

each tail is compared with

$$ \frac{\alpha}{2}=0.025. $$

The interpretation is:

• if $\delta_1<0.025$, the observed count is unusually large and the forecast underpredicts the number of earthquakes;
• if $\delta_2<0.025$, the observed count is unusually small and the forecast overpredicts the number of earthquakes;
• otherwise, the observed count is considered consistent with the forecast.

The forecast predictive interval is obtained from the corresponding lower and upper empirical quantiles of ForecastCounts.

---

3. Magnitude-Test

3.1 Purpose

The Magnitude-Test evaluates whether the magnitude-frequency distribution of the observed catalog is consistent with the simulated forecast catalogs.

The catalog-based test uses:

1. the union of all simulated catalogs;
2. every individual simulated catalog;
3. the observed catalog.

The union forecast catalog is

$$ \Lambda_U=\Lambda_1\cup\Lambda_2\cup\cdots\cup\Lambda_J, $$

and its total number of earthquakes is

$$ N_U=\sum_{j=1}^{J}N_j. $$

Let:

• $\Lambda_U^{(m)}(k)$ be the union-catalog count in magnitude bin $k$;
• $\Lambda_j^{(m)}(k)$ be the count in bin $k$ of simulated catalog $j$;
• $\Omega^{(m)}(k)$ be the observed count in magnitude bin $k$.

---

3.2 Normalization

The union and simulated magnitude histograms are normalized to the observed catalog size. The normalized union count in bin $k$ is

$$ \widetilde{\Lambda}_U^{(m)}(k)=\frac{N_{\mathrm{obs}}}{N_U} \Lambda_U^{(m)}(k). $$

For simulated catalog $j$,

$$ \widetilde{\Lambda}_j^{(m)}(k)=\frac{N_{\mathrm{obs}}}{N_j} \Lambda_j^{(m)}(k). $$

This normalization removes differences in total event count, allowing the test to focus on the shape of the magnitude distribution.

Empty simulated catalogs cannot be normalized and are therefore excluded from the Magnitude-Test statistic distribution.

---

3.3 Observed statistic

The observed discrepancy is calculated from the squared logarithmic residuals:

$$ d_{\mathrm{obs}}= \sum_{k=1}^{K} \left( \log\left[ \frac{N_{\mathrm{obs}}}{N_U} \Lambda_U^{(m)}(k)+1 \right]-\log\left[ \Omega^{(m)}(k)+1 \right] \right)^2. $$

The value one is added to every bin to avoid the singularity associated with $\log(0)$.

Using logarithmic counts reduces the dominance of low-magnitude bins. Because small earthquakes are much more numerous than large earthquakes, a statistic based directly on count differences can be controlled almost entirely by the smallest magnitude bins. Logarithmic residuals make discrepancies in less frequent, relatively large earthquakes more visible.

---

3.4 Simulated statistics

For every valid simulated catalog, the corresponding statistic is

$$ D_j=\sum_{k=1}^{K} \left( \log\left[ \frac{N_{\mathrm{obs}}}{N_U} \Lambda_U^{(m)}(k)+1 \right]-\log\left[ \frac{N_{\mathrm{obs}}}{N_j} \Lambda_j^{(m)}(k)+1 \right] \right)^2, \qquad j = 1,\ldots,J. $$

The values

$$ D_1,D_2,\ldots,D_J $$

form the forecast test distribution.

Because $D$ is a discrepancy statistic, larger values indicate poorer agreement.

The empirical quantile is

$$ \gamma_M=P(D_j\leq d_{\mathrm{obs}}), $$

and the upper-tail probability is

$$ p_{\mathrm{upper}}=P(D_j\geq d_{\mathrm{obs}}). $$

Without a finite-sample correction,

$$ p_{\mathrm{upper}}\approx1-\gamma_M. $$

The forecast is considered inconsistent with the observed magnitude distribution when

$$ p_{\mathrm{upper}}<\alpha. $$

---

4. Spatial-Test

4.1 Purpose

The Spatial-Test evaluates whether the observed earthquake locations are consistent with the spatial probability distribution generated by the forecast catalogs.

The total number of earthquakes is removed from the comparison. Therefore, the test focuses on the relative spatial distribution rather than the total forecast rate.

Let $n_{kj}$ be the number of earthquakes in spatial cell $k$ of simulated catalog $j$. The ensemble spatial probability can be estimated as

$$ p_k=\frac{ \sum_{j=1}^{J}n_{kj} }{ \sum_{l=1}^{K}\sum_{j=1}^{J}n_{lj} }. $$

Thus,

$$ \sum_{k=1}^{K}p_k=1. $$

In the demonstration, a small positive spatial background is included during simulation. This prevents impossible zero-probability cells, so

SpatialOptions.ProbabilityFloor=0;

can be used safely.

---

4.2 Observed spatial score

Let $n_k^{\mathrm{obs}}$ denote the observed count in cell $k$. The normalized observed log score is

$$ S_{\mathrm{obs}}=\frac{1}{N_{\mathrm{obs}}} \sum_{k=1}^{K} n_k^{\mathrm{obs}}\log(p_k). $$

This is equivalent to averaging the logarithmic forecast probability over all observed earthquake locations.

A larger, or less negative, value indicates that the observed earthquakes occurred preferentially in cells assigned relatively high forecast probabilities.

A very small value indicates that many observed earthquakes occurred in cells assigned low probabilities by the forecast.

---

4.3 Simulated spatial scores

For each simulated catalog $j$, the corresponding normalized spatial score is

$$ S_j=\frac{1}{N_j} \sum_{k=1}^{K} n_{kj}\log(p_k). $$

The simulated scores form the Spatial-Test distribution:

$$ S_1,S_2,\ldots,S_J. $$

The lower-tail quantile is

$$ \gamma_S=P(S_j\leq S_{\mathrm{obs}}). $$

Because poor spatial forecasts produce unusually small log scores, the forecast is considered spatially inconsistent when

$$ \gamma_S<\alpha. $$