The Tonga–Fiji region is one of the most seismically active places on Earth. Earthquakes there happen constantly, but most are small.

Plan: Question → Statistical idea → Mathematics → R → Visualization → Insight

Data: R’s built-in quakes data set: 1,000 earthquakes of body-wave magnitude 4.0 or more near Fiji, with latitude, longitude, depth (km), magnitude, and the number of seismic stations that reported each one.

Press Play or drag the slider to go 100 km deeper each step; hover for details. Above 300 km, quakes spread across two plate boundaries; below 300 km they collapse into one narrow band, the Pacific plate sinking under Tonga to almost 700 km.

Small earthquakes should be the most common, yet counts rise until M 4.5 and only then fall. Is nature strange here, or is something wrong with how we see small quakes? Hold that thought.

Seismologists Gutenberg and Richter (1944) found that the number of earthquakes of magnitude at least \(M\) follows

\[ \log_{10} N(\geq M) = a - bM \]

• Each +1 in magnitude makes quakes about \(10^{b}\) times rarer. Worldwide, \(b \approx 1\), so ten times rarer.
• \(a\) measures overall activity; \(b\) measures the balance of small to large quakes.

Statistically, this is the same as saying that above a completeness magnitude \(M_c\), magnitudes follow an exponential distribution:

\[ f(m) = \beta\, e^{-\beta (m - M_c)}, \quad m \geq M_c, \qquad \beta = b \ln 10 \]

So “estimate the law” becomes point estimation: estimate \(b\) from sampled magnitudes, then attach an interval and test it.

For \(n\) independent magnitudes \(m_1, \dots, m_n \geq M_c\), the likelihood and log-likelihood are

\[ L(\beta) = \prod_{i=1}^{n} \beta\, e^{-\beta (m_i - M_c)} = \beta^{\,n} e^{-\beta \sum_{i}(m_i - M_c)}, \qquad \ell(\beta) = n \ln \beta - \beta \sum_{i=1}^{n} (m_i - M_c) \]

Setting \(\ell'(\beta) = \frac{n}{\beta} - \sum_i (m_i - M_c) = 0\) gives Aki’s (1965) estimator:

\[ \hat\beta = \frac{1}{\bar m - M_c} \quad\Longrightarrow\quad \hat b = \frac{\log_{10} e}{\bar m - M_c} \]

Magnitudes are rounded to \(\Delta = 0.1\), so we use Utsu’s correction: replace \(M_c\) by \(M_c - \Delta/2\).

Exact confidence interval. If \(M_i - M_c \sim \text{Exp}(\beta)\), then \(2\beta \sum_i (M_i - M_c) \sim \chi^2_{2n}\), so a 95% CI is

\[ \frac{\chi^2_{0.025,\,2n}}{2 S \ln 10} \;\leq\; b \;\leq\; \frac{\chi^2_{0.975,\,2n}}{2 S \ln 10}, \qquad S = \sum_{i=1}^{n} (m_i - M_c) \]

Standard error and test of the global value \(H_0: b = 1\) vs. \(H_a: b \neq 1\):

\[ SE(\hat b) \approx \frac{\hat b}{\sqrt n}, \qquad z = \frac{\hat b - 1}{SE(\hat b)} \;\overset{H_0}{\approx}\; N(0,1) \]

Comparing two groups (e.g., depths) with \(H_0: b_1 = b_2\):

\[ z = \frac{\hat b_1 - \hat b_2}{\sqrt{SE(\hat b_1)^2 + SE(\hat b_2)^2}} \]

library(ggplot2); library(plotly)
delta <- 0.1                                  # magnitude rounding

# MLE of b (Aki 1965 + Utsu correction) with exact chi-square 95% CI
b_mle <- function(m, Mc) {
  m <- m[m >= Mc]; n <- length(m)
  S <- sum(m - (Mc - delta / 2))
  b <- n / S / log(10)
  ci <- qchisq(c(0.025, 0.975), df = 2 * n) / (2 * S) / log(10)
  data.frame(Mc, n, b, lower = ci[1], upper = ci[2], se = b / sqrt(n))
}

Mc  <- as.numeric(names(which.max(table(quakes$mag))))   # max curvature
res <- b_mle(quakes$mag, Mc)

z <- (res$b - 1) / res$se;  p <- 2 * pnorm(-abs(z))    # H0: b = 1
lm(stations ~ mag, data = quakes)                       # detection model

plot_ly(quakes, x = ~long, y = ~lat, color = ~mag,
        frame = ~cut(depth, seq(0, 700, 100))) %>%     # depth slider map
  add_markers()

The regression predicts only 5 stations for an M 4.0 quake, yet every event in the catalog was reported by at least 10. Many small quakes simply never made it in. Below the peak, the catalog is incomplete, so we estimate \(b\) only from \(M \geq M_c =\) 4.5.

Using every quake (cutoff 4.0) gives \(\hat b =\) 0.65, badly biased by the missing small events. The estimate keeps rising for high cutoffs too, a hint that something also happens at the top end (coming up).

• Test of \(H_0: b = 1\): \(z =\) 1.84, p-value \(=\) 0.066. At \(\alpha = 0.05\) we fail to reject: Fiji’s quakes are consistent with the global ratio, and the CI just includes 1.
• Two estimators, two answers. Least squares weights the last few sparse, heavily correlated tail points as much as the hundreds of small events, so it overshoots. The MLE uses each quake once.

A straight line on this log scale means the law holds. It does from 4.5 to about 5.5, then the observed counts drop below the line. This matches a known limitation of body-wave magnitude: it saturates for large earthquakes, squeezing the biggest events into lower values.

If the law holds, then out of all quakes with \(M \geq 4.5\):

\[ P(M \geq 6 \mid M \geq 4.5) = 10^{-\hat b\,(6 - 4.5)} = 10^{-1.079 \times 1.5} \approx 0.024 \]

Seismologists call earthquakes below 300 km deep-focus. We split the catalog there and compare the shallower and deep groups.

Testing \(H_0: b_{\text{shallower}} = b_{\text{deep}}\): \(z =\) 2.24, p-value \(=\) 0.025. At \(\alpha = 0.05\) we reject: deep quakes have a significantly higher \(b\), meaning large events are relatively rarer far down in the sinking slab.

• R Core Team. quakes data set, package datasets (source: Harvard PRIM-H project, Dr. John Woodhouse). See ?quakes.
• Gutenberg, B. and Richter, C. F. (1944). Frequency of earthquakes in California. Bulletin of the Seismological Society of America, 34, 185–188.
• Aki, K. (1965). Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits. Bulletin of the Earthquake Research Institute, 43, 237–239.
• Utsu, T. (1966). A statistical significance test of the difference in b-value between two earthquake groups. Journal of Physics of the Earth, 14, 37–40.
• Wiemer, S. and Wyss, M. (2000). Minimum magnitude of completeness in earthquake catalogs. Bulletin of the Seismological Society of America, 90, 859–869.
• Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis. Springer. Sievert, C. (2020). Interactive Web-Based Data Visualization with R, plotly, and shiny. CRC Press.

Built with R Markdown, ggplot2, and plotly.