How to estimate b-value and completeness in Python

Problem statement

Ten years of ComCat for a 5° by 3° box of central California, 2015 to 2024, return 96,146 earthquakes of magnitude 1 and above. Fitting log₁₀N = a − bM to them gives the b-value, which says how the events divide between small and large, and the a-value, which sets the rate — the two numbers a hazard model uses to extrapolate from the many small earthquakes a network records to the rare large ones it has not yet seen.

Estimating them is three lines of arithmetic and one difficult decision: the magnitude of completeness, below which the catalogue is missing events. Set it too low and the missing small events flatten the curve and bias b downwards; set it too high and there is nothing left to fit.

Quick answer

Find the completeness magnitude, then use the Aki maximum-likelihood estimator above it:

import numpy as np, pandas as pd, requests

def fetch(bbox, start, end, minmag):
    r = requests.get("https://earthquake.usgs.gov/fdsnws/event/1/query",
                     params={"format": "geojson", "starttime": start, "endtime": end,
                             "minlongitude": bbox[0], "minlatitude": bbox[1],
                             "maxlongitude": bbox[2], "maxlatitude": bbox[3],
                             "minmagnitude": minmag, "orderby": "time-asc", "limit": 20000},
                     timeout=180)
    f = r.json()["features"]
    return pd.DataFrame({"mag": [x["properties"]["mag"] for x in f],
                         "time": pd.to_datetime([x["properties"]["time"] for x in f],
                                                unit="ms", utc=True)})

eq = pd.concat([fetch([-122.5, 35.5, -117.5, 38.5], f"{y}-01-01", f"{y}-12-31T23:59:59", 1.0)
                for y in range(2015, 2025)], ignore_index=True)

def b_value(mags, mc, binwidth=0.1):
    """Aki (1965) maximum likelihood, with the Utsu bin correction."""
    m = mags[mags >= mc - binwidth / 2]
    b = np.log10(np.e) / (m.mean() - (mc - binwidth / 2))
    return b, b / np.sqrt(len(m)), len(m)

counts, edges = np.histogram(eq.mag, bins=np.arange(0.0, 8.0, 0.1))
mc = edges[np.argmax(counts)]                       # maximum-curvature estimate
b, err, n = b_value(eq.mag.values, mc)
print(f"{len(eq):,} events; completeness M{mc:.1f}; b = {b:.2f} ± {err:.2f} from {n:,} events")

96,146 events; completeness M1.1; b = 0.77 ± 0.00 from 88,296 events

The standard error is b divided by the square root of the sample size, so halving the uncertainty needs four times as many events above completeness.

Small events, completeness, a fitted b and its uncertainty

1. Assemble a catalogue with plenty of small events

Fit magnitudes above completeness, so the lower the catalogue reaches, the more events there are to fit.

2. Estimate completeness

The maximum-curvature method — the peak of the non-cumulative histogram — is quick and slightly biased low. A goodness-of-fit test, choosing the lowest magnitude at which a fitted line reproduces the observed counts within a tolerance, is more robust.

3. Fit b above it

Aki's maximum-likelihood estimator is the standard: b = log₁₀(e) / (mean magnitude − completeness), with a half-bin correction.

4. Report the uncertainty

b / √N, where N is the number of events used — and a bootstrap that re-estimates Mc as well (Example 4), because that formula assumes the cut-off is right. A b-value from 50 events is not a measurement.

5. Check the sensitivity to completeness

Recompute b for a range of completeness magnitudes. If it moves a lot, the estimate is being driven by the choice, not by the data.

6. Decide about aftershocks

Aftershock sequences have their own b-values and can dominate a catalogue. Hazard work fits declustered catalogues; report which you used.

Completeness two ways and the cut-off's effect

Example 1 — completeness two ways

def gof_completeness(mags, candidates=np.arange(1.0, 4.0, 0.1), tolerance=5.0):
    """Lowest Mc whose fitted GR line reproduces observed counts within `tolerance` percent."""
    for mc in candidates:
        m = mags[mags >= mc]
        if len(m) < 50:
            break
        b, _, _ = b_value(mags, mc)
        bins = np.arange(mc, m.max() + 0.1, 0.1)
        obs, _ = np.histogram(m, bins=bins)
        pred = len(m) * (10 ** (-b * (bins[:-1] - mc)) - 10 ** (-b * (bins[1:] - mc)))
        resid = 100 * np.abs(obs - pred).sum() / obs.sum()
        if resid <= tolerance:
            return mc, resid
    return np.nan, np.nan

mc_max, = [edges[np.argmax(counts)]]
mc_gof, resid = gof_completeness(eq.mag.values)
print(f"maximum curvature: M{mc_max:.1f}")
print(f"goodness of fit:   M{mc_gof:.1f} (residual {resid:.1f}%)")
for mc in (mc_max, mc_gof):
    b, err, n = b_value(eq.mag.values, mc)
    print(f"  Mc {mc:.1f}: b = {b:.2f} ± {err:.2f} from {n:,} events")

maximum curvature: M1.1
goodness of fit:   M1.4 (residual 4.3%)
  Mc 1.1: b = 0.77 ± 0.00 from 88,296 events
  Mc 1.4: b = 0.82 ± 0.00 from 54,145 events

Example 2 — sensitivity to the cut-off

for mc in np.arange(1.2, 3.2, 0.2):
    b, err, n = b_value(eq.mag.values, mc)
    print(f"Mc {mc:.1f}: b = {b:.2f} ± {err:.2f}  ({n:,} events)")

Mc 1.2: b = 0.79 ± 0.00  (75,324 events)
Mc 1.4: b = 0.83 ± 0.00  (54,757 events)
Mc 1.6: b = 0.86 ± 0.00  (38,656 events)
Mc 1.8: b = 0.88 ± 0.01  (26,588 events)
Mc 2.0: b = 0.89 ± 0.01  (17,915 events)
Mc 2.2: b = 0.91 ± 0.01  (12,193 events)
Mc 2.4: b = 0.93 ± 0.01  (8,152 events)
Mc 2.6: b = 0.96 ± 0.01  (5,456 events)
Mc 2.8: b = 0.99 ± 0.02  (3,604 events)
Mc 3.0: b = 0.97 ± 0.02  (2,272 events)

Example 3 — what b buys you: extrapolating the rate

years = 10
for mc in (mc_max, mc_gof, 2.6):
    b, _, n = b_value(eq.mag.values, mc)
    rate_mc = n / years
    print(f"from Mc {mc:.1f} (b = {b:.2f}):")
    for target in (4.0, 5.0, 6.0, 7.0):
        rate = rate_mc * 10 ** (-b * (target - mc))
        observed = (eq.mag >= target).sum() / years
        print(f"  M{target}+: predicted {rate:6.2f} per year, observed {observed:5.2f} "
              f"(one every {1 / rate:5.1f} years predicted)")

from Mc 1.1 (b = 0.77):
  M4.0+: predicted  51.47 per year, observed 19.60 (one every   0.0 years predicted)
  M5.0+: predicted   8.73 per year, observed  1.60 (one every   0.1 years predicted)
  M6.0+: predicted   1.48 per year, observed  0.30 (one every   0.7 years predicted)
  M7.0+: predicted   0.25 per year, observed  0.10 (one every   4.0 years predicted)
from Mc 1.4 (b = 0.82):
  M4.0+: predicted  40.69 per year, observed 19.60 (one every   0.0 years predicted)
  M5.0+: predicted   6.20 per year, observed  1.60 (one every   0.2 years predicted)
  M6.0+: predicted   0.95 per year, observed  0.30 (one every   1.1 years predicted)
  M7.0+: predicted   0.14 per year, observed  0.10 (one every   6.9 years predicted)
from Mc 2.6 (b = 0.94):
  M4.0+: predicted  25.88 per year, observed 19.60 (one every   0.0 years predicted)
  M5.0+: predicted   2.97 per year, observed  1.60 (one every   0.3 years predicted)
  M6.0+: predicted   0.34 per year, observed  0.30 (one every   2.9 years predicted)
  M7.0+: predicted   0.04 per year, observed  0.10 (one every  25.4 years predicted)

The two automatic estimates put completeness at M1.1 and M1.4, and both overpredict M4 and above by a factor of two or three. The sensitivity table is the tell: b keeps climbing until about M2.6 and then settles near 0.96, which is where the catalogue is genuinely complete. Fitted from there (b = 0.94), M4+ comes out at 25.9 a year against 19.6 observed and M6+ at 0.34 against 0.30 — much closer, and M7 now looks rarer than the decade's single event, which is what one event in ten years should look like. Automatic completeness estimators are a starting point, not an answer — and this catalogue also contains an aftershock sequence that no single b-value describes well.

Example 4 — an uncertainty that includes the choice of Mc

b / √N assumes the cut-off is right. Resampling the catalogue and re-running both the completeness estimate and the fit gives an uncertainty that includes the cut-off:

rng = np.random.default_rng(42)
mags = eq.mag.values
boot = []
for _ in range(200):
    s = rng.choice(mags, size=len(mags), replace=True)
    mc, _ = gof_completeness(s)
    boot.append((round(float(mc), 1), b_value(s, mc)[0]))
boot = pd.DataFrame(boot, columns=["mc", "b"])
print("bootstrap Mc:", boot.mc.value_counts().sort_index().to_dict())
print(f"bootstrap b: mean {boot.b.mean():.3f}, sd {boot.b.std():.3f}, "
      f"95% range {boot.b.quantile(0.025):.3f} to {boot.b.quantile(0.975):.3f}")
print(f"formula at Mc {mc_gof:.1f}: ± {b_value(mags, mc_gof)[1]:.3f}")

bootstrap Mc: {1.4: 166, 1.5: 24, 1.8: 1, 2.0: 3}
bootstrap b: mean 0.820, sd 0.011, 95% range 0.811 to 0.843
formula at Mc 1.4: ± 0.004

The formula's ±0.004 is the uncertainty if M1.4 is right. Resampled, the goodness-of-fit estimator lands on M1.4 in 166 of 200 runs and anywhere up to M2.0 in the rest, and the spread of b nearly triples to 0.011. On a catalogue this size the cut-off, not the sample, is where the uncertainty lives — and Example 2 shows the true completeness is higher still.

Why maximum likelihood, not regression

Why completeness bends the curve

Below the completeness magnitude a network records only some of the events — the closer ones, the better-recorded ones. The counts therefore fall away from the straight line, and a fit that includes them is pulled towards a shallower slope. That underestimates b, which in hazard terms overestimates the relative frequency of large events.

Why the maximum-likelihood estimator, not a regression

Fitting a line to log-counts by least squares weights every magnitude bin equally, although the highest bins contain a handful of events and the lowest contain thousands. The Aki estimator uses the mean magnitude above completeness directly, which is the maximum-likelihood solution for an exponential distribution and has a known standard error.

Why b near 1 matters

A b-value of 1 means ten times fewer events per magnitude unit. Below 1, large events are relatively more common than that; above 1, less. Changes in b are used as a stress proxy in volcano and induced-seismicity monitoring, which is why the uncertainty — and the sensitivity to the completeness choice — needs reporting with it.

Binning, mixed magnitudes and small samples

• Magnitude binning matters: the half-bin correction assumes magnitudes are rounded to 0.1.
• Mixed magnitude types widen the distribution and bias b; filter to one type where possible.
• Completeness varies in time; a long catalogue may need a different Mc per period.
• Aftershocks raise the apparent rate and can change b; state whether the catalogue was declustered.
• Maximum magnitude is not part of the relation; hazard models truncate it separately.
• b from fewer than about 100 events above completeness is not worth quoting.

Internal links

• Earthquake catalogues explained — completeness, magnitude types, declustering
• How to download and map an earthquake catalogue in Python — getting the data
• Earthquake counts rise over time in the catalogue — the detection artefact
• Return periods and annual exceedance probability explained — rates as return periods
• Hazard, exposure and vulnerability explained — where rates feed into risk
• How to fit return periods to annual maxima in Python — the flood analogue

FAQ

What is the b-value?

The slope of the frequency–magnitude relation: how quickly earthquake numbers fall as magnitude rises. It is usually near 1.

What is the magnitude of completeness?

The lowest magnitude above which the catalogue records essentially every event.

How do I estimate b?

With the Aki maximum-likelihood estimator above the completeness magnitude, not by regressing log counts.

How many events do I need?

At least a hundred above completeness for a usable estimate; the standard error is b/√N.

Should I decluster first?

For hazard rates, yes. Say which catalogue the b-value came from.

Can I use the fit to predict big earthquakes?

It predicts long-run rates, not individual events, and the extrapolation is only as good as the completeness and the assumption of stationarity.