Is this fish being overfished? Eighteen stocks as savings accounts, the biology we bring, and a cod collapse that catches alone cannot see

import os

for v_ in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
    os.environ[v_] = "1"  # nutpie runs chains as threads: one BLAS pool per chain would fight

import logging
import sys
import textwrap
import time
import warnings

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import nutpie
import pandas as pd
import pymc as pm
import pytensor
import pytensor.tensor as pt
from IPython.display import Markdown, display
from matplotlib.colors import to_rgb
from matplotlib.patches import Ellipse, Polygon, Rectangle

from pymc_challenges import data
from pymc_challenges import public as pub

RANDOM_SEED = 114
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-variat")
logging.getLogger("pymc").setLevel(logging.WARNING)
warnings.filterwarnings("ignore", category=FutureWarning)
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)
sys.modules.pop("jax", None)  # arviz leaves a lazy jax stub that makes matplotlib import all of JAX
pd.set_option("display.width", 160)
T_START = time.time()
TIMES = {}
INK, MUTED, LIGHT, HI, LO = pub.INK, pub.MUTED, pub.LIGHT, pub.HI, pub.LO


def wrap(text, width=90):
    return "\n".join(textwrap.fill(t, width) for t in text.split("\n"))


print(f"PyMC {pm.__version__}, ArviZ {az.__version__}, nutpie {nutpie.__version__}")

PyMC 6.3.2, ArviZ 1.3.0, nutpie 0.16.11

#In plain words

The question. When I buy cod, sole or hake, is the population in the sea healthy, or has it been fished down?

What we did. For 18 fish stocks from North America, Chile, Australia and New Zealand we rebuilt the hidden amount of fish in the sea each year from catches and surveys. Each stock is a savings account: it earns the most interest at half its natural size, and catches are withdrawals. How fast each kind of fish grows is built in from biology.

What we found. In the last year of each stock's data (2001-2015, not today), about 6 of the 18 (roughly 3 to 9) were below the healthy level, the size that gives the largest catch that can be taken year after year. Only one was more likely than not fished too hard. Catches alone missed the northern cod collapse.

HOW_SURE = display(Markdown("**How sure are we.** *Filled in from the fits in part E.*"),
                   display_id=True)
HEADLINE = display(Markdown("*The headline figure is drawn in part G from the fits below and "
                            "appears here once the notebook has run.*"), display_id=True)

#The setting

On 2 July 1992 Canada closed the northern cod fishery off Newfoundland. A stock that had supported fishers for five centuries had fallen, in the official reconstruction, to about a twentieth of the level that gives the largest sustainable catch, and some 30,000 people lost their work. Hutchings & Myers (1994, Can. J. Fish. Aquat. Sci.) asked how the scientists had missed it: catch rates in the commercial fishery held up while the fish shrank, because boats kept finding the last dense schools. What is caught is not what is in the sea, and the gap between them has to be modelled, with its uncertainty stated.

The oldest model for that gap is a surplus-production model (Schaefer 1954; Pella & Tomlinson 1969). A stock of biomass $B$ grows logistically towards a carrying capacity $K$ at intrinsic rate $r$, and the catch $C$ is taken out:

$$B_{t+1} = B_t + r B_t \left(1 - \frac{B_t}{K}\right) - C_t .$$

The growth term is the "interest" of a savings account. It is largest when the stock is at half its unfished size, $B_{MSY} = K/2$, where it equals the maximum sustainable yield $MSY = rK/4$; a fishing mortality above $F_{MSY} = r/2$ drains the account. The ratios $B/B_{MSY}$ and $F/F_{MSY}$ are the coordinates of the Kobe plot that tuna commissions and many national agencies use to report stock status.

We never see $B$. We see the catch, and an index (a research survey, or the catch per hour of the fleet) assumed proportional to $B$ with noise. Meyer & Millar (1999, Can. J. Fish. Aquat. Sci.) wrote this as a Bayesian state-space model: each year's biomass is a latent state with process error (good and bad years for the fish), and the index observes it with observation error. Punt & Hilborn (1997, Reviews in Fish Biology and Fisheries) made the case for Bayesian assessment; JABBA (Winker, Carvalho & Kapur 2018, Fisheries Research) is a modern version. The prior on $r$ is where biology enters: CMSY (Froese et al. 2017, Fish and Fisheries) turns FishBase's resilience category into a range for $r$, and we start from the same ranges.

The data come from the RAM Legacy Stock Assessment Database (Ricard, Minto, Jensen & Baum 2012, Fish and Fisheries), which collects the inputs and outputs of formal assessments around the world, so for 14 of our stocks there is an official answer to compare with. Those are usually age-structured models with far more data: agreement is a check, not proof, and disagreement is a question, not a verdict on either.

part	question	tool
A	What do catches and an index look like?	Pacific cod, Hecate Strait; the surplus curve
B	What can a stock's history tell us?	a toy stock, one-way and two-way trips
C	How does one real stock fit?	centred free states, a prior check, two failures
D	What do 18 stocks share?	family priors for r, one process noise, official assessments
E	Below the healthy level, and still fished hard?	Kobe probabilities, one decisive noise
F	What if we only had the catches?	northern cod, catch-only refits
G	How do we explain it to everyone?	40 possible histories, fish icons, an MSY band, a GIF

#A · Catches, an index and a savings account

#The data

tools/build_e114_ram.py downloads the RAM v4.44 archive, reads its R tables once with Rscript, and keeps one row per stock and year from the first to the last catch year. We kept stocks with at least 25 years of catch and at least 15 years of an index, and the script checks each index against the assessment's own biomass, where RAM has one, to catch copies: the Peruvian anchoveta's "survey" equals its assessed biomass in every year, so it is not an independent index and the anchoveta stays in the file only as a catch-only row. Each stock carries its FishBase resilience and the matching range for $r$. Beyond the copy test the choice was a judgement made while looking at the data: three more candidates (Greenland halibut, Gulf of St Lawrence cod, St Pierre pollock) were left out after we saw how their index compared with the assessed biomass (correlations 0.03 and 0.54, and one with nothing to compare), while Hecate Strait rock sole was kept although its catch rate runs against the assessed biomass (correlation -0.17; its short survey agrees better, 0.40).

data.describe("e114_ram_stocks")
RAM = data.load("e114_ram_stocks")
NICE = {
    "PCODHS": "Pacific cod, Hecate Strait", "PCODWCVANI": "Pacific cod, W. Vancouver Island",
    "PCOD5AB": "Pacific cod, Queen Charlotte Sound", "RSOLEHSTR": "Rock sole, Hecate Strait",
    "RSOLE5AB": "Rock sole, Queen Charlotte Sound", "YELL3LNO": "Yellowtail flounder, Grand Banks",
    "YELLGB": "Yellowtail flounder, Georges Bank",
    "WITFLOUN4RST": "Witch flounder, Gulf of St Lawrence",
    "SPHAKECH": "Hake, Chile", "SHAKE4VWX": "Silver hake, Scotian Shelf",
    "SHAKESGBMATL": "Silver hake, S. Georges Bank-Mid-Atlantic",
    "SHAKEGOMNGB": "Silver hake, Gulf of Maine", "REDFISHSPP3LN": "Redfish, Grand Banks",
    "NZSNAPNZ7": "Snapper, New Zealand (SNA 7)", "TREVALLYTRE7": "Trevally, New Zealand (TRE 7)",
    "NZLINGESE": "Ling, SE Australia", "BLSHARNPAC": "Blue shark, North Pacific",
    "WHAKE4T": "White hake, S. Gulf of St Lawrence", "ANCHOSA": "Anchovy, South Africa",
    "SARDSA": "Sardine, South Africa", "COD2J3KL": "Atlantic cod, northern (2J3KL)",
    "PANCHPERUNC": "Peruvian anchoveta"}
summ = (RAM.groupby("stockid")
        .agg(years=("year", lambda y: f"{y.min()}-{y.max()}"), cpue=("cpue", "count"),
             survey=("survey", "count"), official=("bdivbmsy", "count"), family=("family", "first"),
             resilience=("resilience", "first"), r_range=("r_lo", "first"), r_hi=("r_hi", "first"),
             max_catch_kt=("catch_t", lambda c: round(c.max() / 1e3, 1)), role=("role", "first")))
summ["r_range"] = [f"{a:.2f}-{b:.2f}" for a, b in zip(summ.r_range, summ.r_hi)]
print(summ.drop(columns="r_hi").sort_values(["role", "family"]).to_string())

e114_ram_stocks: One row per stock and year from the first to the last catch year (1,278 rows, 22 stocks): stockid, year, catch_t (total catch TCbest, tonnes), cpue and survey (RAM's preferred abundance indices, units in cpue_units / survey_units), bdivbmsy and udivumsy (the assessment's B/Bmsy and U/Umsy), tb_t (assessment total biomass, tonnes); stock columns repeated per row: stocklong, commonname, scientificname, family, region, assessorid, assessyear, role ('index': 20 finfish stocks with an index of >= 15 years; 'catch_only': northern cod COD2J3KL and Peruvian anchoveta PANCHPERUNC, whose RAM survey series equals the assessment's own biomass and is dropped), resilience (FishBase), r_lo, r_hi (Froese et al. 2017 prior range of r).
  source: RAM Legacy Stock Assessment Database v4.44 (Ricard, Minto, Jensen & Baum 2012, Fish and Fisheries 13:380-398), Zenodo DOI 10.5281/zenodo.2542919, licence CC BY 4.0; retrieved 2026-11-05. Resilience categories from FishBase species summaries (read 2026-11-06), with the CMSY r ranges of Froese et al. (2017). DERIVED by tools/build_e114_ram.py (needs R once, to read DBdata.RData); the url is the 55 MB zip, not this table - keep the cached file in data/. A newer v4.66 (2024, DOI 10.5281/zenodo.14043031) exists.
  url:    https://zenodo.org/records/2542919/files/RLSADB%20v4.44.zip?download=1
                   years  cpue  survey  official          family    resilience    r_range  max_catch_kt        role
stockid                                                                                                            
PANCHPERUNC    1959-2016     0       0        53     Engraulidae          High  0.60-1.50       10900.0  catch_only
COD2J3KL       1959-2014     0       0        50         Gadidae        Medium  0.20-0.80         810.0  catch_only
TREVALLYTRE7   1944-2014    23       0        71      Carangidae           Low  0.05-0.50           2.3       index
BLSHARNPAC     1971-2012    18       0        40  Carcharhinidae           Low  0.05-0.50         115.0       index
SARDSA         1949-2015     0      32         0       Clupeidae        Medium  0.20-0.80         410.0       index
ANCHOSA        1949-2015     0      32         0     Engraulidae        Medium  0.20-0.80         600.0       index
PCOD5AB        1945-2012    57       0         0         Gadidae           Low  0.05-0.50           3.2       index
PCODHS         1944-2012    57      11        57         Gadidae           Low  0.05-0.50           9.5       index
PCODWCVANI     1956-2001    46       0        46         Gadidae           Low  0.05-0.50           5.8       index
SHAKE4VWX      1977-2014     0      38        22    Merlucciidae        Medium  0.20-0.80          87.7       index
SHAKEGOMNGB    1963-2013    47       0         0    Merlucciidae        Medium  0.20-0.80          94.5       index
SHAKESGBMATL   1963-2013    47       0         0    Merlucciidae        Medium  0.20-0.80         593.0       index
SPHAKECH       1968-2015    31       0        48    Merlucciidae        Medium  0.20-0.80         130.0       index
NZLINGESE      1968-2007    21       0        40      Ophidiidae           Low  0.05-0.50           1.3       index
WHAKE4T        1960-2010     0      40         0        Phycidae           Low  0.05-0.50          14.0       index
RSOLE5AB       1945-2013    47       7        69  Pleuronectidae           Low  0.05-0.50           1.1       index
RSOLEHSTR      1945-2005    52      12        61  Pleuronectidae           Low  0.05-0.50           2.6       index
WITFLOUN4RST   1960-2011     0      41        51  Pleuronectidae        Medium  0.20-0.80           6.9       index
YELL3LNO       1960-2014    45       0        50  Pleuronectidae           Low  0.05-0.50          39.3       index
YELLGB         1935-2014     0      28        35  Pleuronectidae           Low  0.05-0.50          21.4       index
REDFISHSPP3LN  1959-2013    36       0        55      Sebastidae  Very low-Low  0.01-0.50          79.0       index
NZSNAPNZ7      1931-2013    25       0        83        Sparidae           Low  0.05-0.50           3.3       index

Twenty stocks have an index, from 18 years (blue shark) to 57 (Pacific cod in Hecate Strait and Queen Charlotte Sound), and fourteen also have the assessment's B/B_MSY. They span eleven families; cod (Gadidae), flatfish (Pleuronectidae) and hake (Merlucciidae) have three to five stocks each, the rest one. FishBase calls most of them "low" resilience, with $r$ between 0.05 and 0.5 a year, and the four Merluccius hakes, witch flounder and the two South African small pelagics "medium" (0.2-0.8). These two swing with ocean conditions more than with fishing, a story a surplus-production model tells badly; they caused most of the divergences in a 20-stock fit, so the main model uses the other 18 and leaves them for "Try it yourself".

#One stock: Pacific cod in Hecate Strait

Hecate Strait lies between the Haida Gwaii islands and the mainland of British Columbia. Its Pacific cod have been fished since the 1940s, with a commercial catch rate (kg per hour of trawling) from 1956 and a research survey for 1984-1994.

def stock_rows(sid):
    return RAM[RAM.stockid == sid].sort_values("year").reset_index(drop=True)


g = stock_rows("PCODHS")
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(12, 4.4), width_ratios=[1.6, 1])
ax.bar(g.year, g.catch_t / 1e3, color=LIGHT, width=0.8)
ax.set_ylabel("catch (thousand tonnes)")
axr = ax.twinx()
for col, c, lab in (("cpue", HI, "trawl catch rate"), ("survey", LO, "research survey")):
    ok = g[col].notna()
    axr.plot(g.year[ok], g[col][ok] / g[col][ok].mean(), "o-", color=c, ms=3, lw=1, label=lab)
axr.set_ylabel("index / its own mean")
axr.grid(False)
axr.legend(loc="upper left", fontsize=10)
ax.set_title("Pacific cod, Hecate Strait: catch (bars), two indices", loc="left", fontsize=12)
P_ = np.linspace(0, 1, 200)
for rr, ls in ((0.13, "-"), (0.4, "--")):
    ax2.plot(P_, rr * P_ * (1 - P_), color=INK, ls=ls, label=f"r = {rr}")
ax2.axvline(0.5, color=HI, lw=1)
ax2.axvline(0.25, color=LO, lw=1, ls="--")
ax2.text(0.52, 0.103, "$B_{MSY}$: the most\ninterest", color=HI, fontsize=10, va="top")
ax2.text(0.24, 0.103, "overfished:\nbelow half\nof $B_{MSY}$", color=LO, fontsize=10, ha="right",
         va="top")
ax2.set_xlabel("stock as a share of its unfished size, $B/K$")
ax2.set_ylabel("surplus per year (share of $K$)")
ax2.legend(loc="lower center", fontsize=10)
ax2.set_title("The savings account", loc="left", fontsize=12);

The catch peaked twice, at about 9,000 t in the mid-1960s and 9,500 t in 1987, and fell to a few hundred tonnes after 2000. The catch rate peaks in the 1970s and is lowest around 2000, with year-to-year jumps of a factor of two: observation error. The eleven-year survey does not obviously agree with it, and for 1944-1955 there is no index at all.

The right panel is the savings account: the surplus is a parabola in $B/K$ that peaks at half the unfished size whatever $r$ is; $r$ only sets the height. The dashed line is the usual "overfished" threshold, half of $B_{MSY}$.

#B · A toy stock: what a one-way trip cannot tell

#Two histories with a known truth

Before a real stock, a simulated one with $r = 0.4$, $K = 1{,}000$ t and so $MSY = 100$ t, an index with catchability $q = 0.001$ and 15% observation noise, and a little process noise. In the one-way trip the catch rises steadily from 10% to 150% of MSY over 40 years, and the stock declines from $K$ to about 40% of it. In the two-way trip the catch rises to 170% of MSY over 20 years and is then cut to 40%, so the stock falls and recovers. Hilborn & Walters (1992, Quantitative Fisheries Stock Assessment) made the point this toy illustrates: data from a fishery that has only ever gone down tell you little about how it would come back.

R_TRUE, K_TRUE, Q_TRUE, T_TOY = 0.4, 1000.0, 0.001, 40
MSY_TRUE = R_TRUE * K_TRUE / 4
tt = np.arange(T_TOY)
TRIPS = {"one-way trip": MSY_TRUE * (0.1 + 1.4 * tt / (T_TOY - 1)),
         "two-way trip": np.where(tt < 20, MSY_TRUE * (0.1 + 1.6 * tt / 19), 0.4 * MSY_TRUE)}
toy_rng = np.random.default_rng(1140)


def simulate(C, gen, sp=0.05, so=0.15):
    P = np.empty(len(C))
    P[0] = 1.0
    for k in range(1, len(C)):
        f = P[k - 1] + R_TRUE * P[k - 1] * (1 - P[k - 1]) - C[k - 1] / K_TRUE
        P[k] = max(f, 0.01) * np.exp(sp * gen.normal())
    return P, Q_TRUE * K_TRUE * P * np.exp(so * gen.normal(size=len(C)))


TOY = {k: (C,) + simulate(C, toy_rng) for k, C in TRIPS.items()}
for k, (C, P, I) in TOY.items():
    print(f"{k}: B/K from {P[0]:.2f} to a low of {P.min():.2f}, ending at {P[-1]:.2f}")

one-way trip: B/K from 1.00 to a low of 0.40, ending at 0.40
two-way trip: B/K from 1.00 to a low of 0.32, ending at 0.86

#A fit with observation error only

The simplest fit, and the one fisheries scientists used for decades, ignores process error: given $r$, $K$ and a start at $K$, the dynamics are deterministic, and only the index is noisy, $\log I_t \sim N(\log q + \log B_t, \sigma_o)$. The recursion runs in pytensor.scan, and a smooth floor $\epsilon \,\mathrm{softplus}(f/\epsilon)$ keeps the stock positive when a draw of $(r, K)$ makes the catch larger than the fish. The priors are deliberately vague: $\log r \sim N(\log 0.4, 1)$ and $K$ around five times the largest catch with a log sd of 1.

def toy_fit(C, I, seed=1):
    cmax = C.max()
    with pm.Model() as m:
        log_r = pm.Normal("log_r", np.log(0.4), 1.0)
        log_K = pm.Normal("log_K", np.log(5 * cmax), 1.0)
        r, K = pt.exp(log_r), pt.exp(log_K)

        def step(c, p, r, K):
            f = p + r * p * (1 - p) - c / K
            return 0.01 * pt.softplus(f / 0.01)

        Ps = pytensor.scan(step, sequences=[pt.as_tensor(C[:-1])],
                           outputs_info=[pt.constant(1.0, dtype="float64")],
                           non_sequences=[r, K], return_updates=False)
        P_all = pt.concatenate([pt.ones(1), Ps])
        lq = pm.Normal("log_qK", np.log(I).mean() - np.log(0.5), 2.0)
        so = pm.LogNormal("sigma_o", np.log(0.2), 0.5)
        pm.Normal("I", lq + pt.log(P_all), so, observed=np.log(I))
        return pm.sample(random_seed=seed, progressbar=False, target_accept=0.9)


t0 = time.time()
TOY_POST = {}
for k, (C, P, I) in TOY.items():
    idt = toy_fit(C, I)
    lr = idt.posterior["log_r"].values.ravel()
    lk = idt.posterior["log_K"].values.ravel()
    sm = az.summary(idt, var_names=["log_r", "log_K"], round_to=4)
    q_r, q_k, q_m = (np.quantile(v, [0.055, 0.5, 0.945])
                     for v in (np.exp(lr), np.exp(lk), np.exp(lr + lk) / 4))
    TOY_POST[k] = (lr, lk)
    print(f"{k}: divergences {int(idt.sample_stats['diverging'].sum())}, "
          f"max r_hat {sm.r_hat.max():.3f}, min ESS {sm.ess_bulk.min():.0f}, "
          f"corr(log r, log K) {np.corrcoef(lr, lk)[0, 1]:.2f}\n"
          f"   r {q_r[1]:.2f} ({q_r[0]:.2f}-{q_r[2]:.2f}), K {q_k[1]:,.0f} ({q_k[0]:,.0f}-"
          f"{q_k[2]:,.0f}), MSY {q_m[1]:.0f} ({q_m[0]:.0f}-{q_m[2]:.0f});  truth r 0.4, K 1,000, "
          f"MSY 100")
TIMES["toy"] = time.time() - t0

one-way trip: divergences 160, max r_hat 1.012, min ESS 186, corr(log r, log K) -0.99
   r 0.35 (0.10-1.08), K 1,213 (467-2,930), MSY 107 (69-128);  truth r 0.4, K 1,000, MSY 100

two-way trip: divergences 0, max r_hat 1.006, min ESS 727, corr(log r, log K) -0.98
   r 0.40 (0.33-0.49), K 971 (837-1,130), MSY 97 (91-103);  truth r 0.4, K 1,000, MSY 100

fig, axs = plt.subplots(2, 2, figsize=(11, 7.4), height_ratios=[1, 1.15])
for j, (k, (C, P, I)) in enumerate(TOY.items()):
    a = axs[0, j]
    a.bar(tt, C / MSY_TRUE, color=LIGHT, width=0.8)
    a.axhline(1, color=MUTED, lw=0.8, ls=":")
    a.text(0.5, 1.03, "catch = MSY", color=MUTED, fontsize=9.5)
    a.plot(tt, 2 * P, color=INK, lw=1.5, label="true B/B$_{MSY}$")
    a.plot(tt, 2 * I / (Q_TRUE * K_TRUE), "o", color=HI, ms=3.5, label="index (scaled)")
    a.set_ylim(0, 2.3)
    a.set_title(f"{k}; bars: catch / MSY", loc="left", fontsize=11.5)
    a.legend(loc="upper right", fontsize=9.5)
    lr, lk = TOY_POST[k]
    b = axs[1, j]
    b.scatter(np.exp(lr[::2]), np.exp(lk[::2]), s=3, alpha=0.25, color=HI)
    b.plot(R_TRUE, K_TRUE, "*", color=LO, ms=15, label="truth")
    rg = np.geomspace(0.05, 3, 50)
    b.plot(rg, 4 * MSY_TRUE / rg, color=MUTED, lw=0.8, ls="--", label="r K / 4 = 100 (true MSY)")
    b.set_xscale("log")
    b.set_yscale("log")
    b.set_xlim(0.08, 1.5)
    b.set_ylim(300, 6000)
    b.set_xlabel("r")
    b.set_ylabel("K (t)")
    b.legend(loc="upper right", fontsize=9.5)
    b.set_title(f"posterior draws, corr(log r, log K) {np.corrcoef(lr, lk)[0, 1]:.2f}",
                loc="left", fontsize=11.5)
fig.suptitle("A one-way trip leaves r and K on a thin ridge; a recovery pins them down",
             x=0.01, ha="left", fontsize=13);

On the one-way trip the posterior is a long, thin banana: $\log r$ and $\log K$ are correlated at $-0.99$, and $r$ could be anything from 0.1 to 1.1 a year, with $K$ moving the other way from about 470 to 2,900 t. The data fix one combination, roughly the surplus the stock produced while it declined (MSY 69-128 t), and say almost nothing about whether that is a fast growth rate on a small stock or a slow one on a big stock. The 160 divergences come from the ends of the banana, where it bends and narrows.

On the two-way trip the recovery after the cut is the missing experiment: how fast the stock came back measures $r$, and the ridge shrinks to a short segment around the truth ($r$ 0.33-0.49, MSY 91-103 t, no divergences). Real fisheries are mostly one-way trips, so information about $r$ has to come from biology and from other stocks of the same kind. That is what part D does.

#C · One real stock in PyMC

#The state-space model

For stock $s$ we work with the depletion $P_{st} = B_{st}/K_s$ and with catches divided by the stock's largest catch $\bar C_s$, so that $K_s$ is in units of the largest catch and every stock has the same scale. The model, after Meyer & Millar (1999), is

$$f_{st} = P_{s,t-1} + r_s P_{s,t-1}(1 - P_{s,t-1}) - C_{s,t-1}/K_s, \qquad \log P_{st} \sim N\big(\log \tilde f_{st},\ \sigma_p\big),$$

$$\log I_{it} \sim N\big(\log (qK)_i + \log P_{s(i),t},\ \sigma_{o,i}\big),$$

where $\tilde f = \epsilon\,\mathrm{softplus}(f/\epsilon)$ is a smooth floor at about $\epsilon = 1\%$ of $K$, and $i$ runs over the index series (Pacific cod in Hecate Strait has two). A stock cannot lose more fish than it has, so a penalty $-\sum \mathrm{softplus}\big((0.02 - f_{st})/0.005\big)$ removes draws in which the catch would take the stock below 2% of $K$; without it the floor alone would let a tiny $K$ explain any catch. The priors:

• $\log P_{s1} \sim N(\log 0.8, 0.3)$: lightly fished when the record starts, an assumption that flatters records starting after heavy fishing (Scotian Shelf silver hake from 1977, after the distant-water fleets; Chilean hake's catches here from 1968);
• $\log K_s \sim N(\log 5, 0.8)$ in units of the largest catch: about 1.4 to 18 times it (89%);
• $\log r_s \sim N(m_s, w_s)$ with $m_s$ the log midpoint of the FishBase-CMSY range and $w_s$ a quarter of its log width, so the range is a 95% interval (part D keeps the centre and changes the width);
• $\sigma_p \sim \mathrm{LogNormal}(\log 0.1, 0.5)$, $\sigma_{o,i} \sim \mathrm{LogNormal}(\log 0.25, 0.5)$, and $\log(qK)_i$ vague around the index's mean.

The states $\log P_{st}$ are free parameters (pm.Flat), and the process and initial densities enter as pm.Potential terms: the centred form. The non-centred form, usually recommended for hierarchical models, pushes standardised innovations through a scan; we try it too. One helper flattens any list of stocks into stock-year arrays, and one builds the model.

def prep(stocks, source=RAM, use_index=True, multi_only=True):
    """Flatten stock-years for several stocks; catches in units of each stock's largest catch."""
    st, years, C, first, prev, cmax, m_r, s_r, fam = [], [], [], [], [], [], [], [], []
    obs_state, obs_ser, logI, ser_name, ser_m, ser_stock = [], [], [], [], [], []
    off = 0
    for s, sid in enumerate(stocks):
        g = source[source.stockid == sid].sort_values("year").reset_index(drop=True)
        T = len(g)
        cmax.append(g.catch_t.max())
        st += [s] * T
        years += list(g.year)
        C += list(g.catch_t / g.catch_t.max())
        first.append(off)
        prev += list(range(off, off + T - 1))
        lo, hi = g.r_lo.iloc[0], g.r_hi.iloc[0]
        m_r.append(0.5 * (np.log(lo) + np.log(hi)))
        s_r.append((np.log(hi) - np.log(lo)) / 4)
        fam.append(g.family.iloc[0])
        for kind in ("cpue", "survey") if use_index else ():
            ok = (g[kind].notna() & (g[kind] > 0)).to_numpy()
            if ok.sum() >= 5:
                li = np.log(g[kind][ok].to_numpy())
                obs_state += list(off + np.where(ok)[0])
                obs_ser += [len(ser_name)] * int(ok.sum())
                logI += list(li)
                ser_name.append(f"{sid}:{kind}")
                ser_m.append(li.mean())
                ser_stock.append(s)
        off += T
    if multi_only:  # family effects only for families with two or more stocks
        fam = [f if fam.count(f) >= 2 else "(single)" for f in fam]
    fams = sorted(set(fam))
    st = np.array(st)
    first = np.array(first)
    last = np.r_[first[1:] - 1, off - 1]
    return dict(stocks=list(stocks), S=len(stocks), N=off, st=st, years=np.array(years),
                C=np.array(C), first=first, last=last, rest=np.setdiff1d(np.arange(off), first),
                prev=np.array(prev), cmax=np.array(cmax), obs_state=np.array(obs_state, int),
                obs_ser=np.array(obs_ser, int), logI=np.array(logI), ser_name=ser_name,
                ser_m=np.array(ser_m), ser_stock=np.array(ser_stock, int), m_r=np.array(m_r),
                s_r=np.array(s_r), fam_names=fams, fam=np.array([fams.index(f) for f in fam]))


def ss_model(d, *, hier=False, sigma_p="stock", sp_mu=np.log(0.1), sp_sd=0.5, k_mu=np.log(5.0),
             k_sd=0.8, eps=0.01, pen=0.005, floor=0.02, end_prior=None):
    """Schaefer state-space model with centred free states (Meyer & Millar 1999).

    sigma_p: "stock" (one per stock), "pooled" (one for all) or a number (fixed).
    end_prior: (lo, hi) range for the final depletion B/K, imposed as a steep wall (CMSY's filter).
    """
    coords = {"stock": d["stocks"], "state": np.arange(d["N"]), "series": d["ser_name"],
              "family": d["fam_names"]}
    st, prev, rest, first = d["st"], d["prev"], d["rest"], d["first"]
    with pm.Model(coords=coords) as m:
        if hier:
            sigma_fam = pm.HalfNormal("sigma_fam", 0.4)
            tau_r = pm.HalfNormal("tau_r", 0.4)
            w = pm.Normal("w_fam", 0, 1, dims="family")
            z = pm.Normal("z_r", 0, 1, dims="stock")
            fe = sigma_fam * w
            if "(single)" in d["fam_names"]:
                fe = pt.set_subtensor(fe[d["fam_names"].index("(single)")], 0.0)
            log_r = pm.Deterministic("log_r", d["m_r"] + fe[d["fam"]] + tau_r * z, dims="stock")
        else:
            log_r = pm.Normal("log_r", d["m_r"], d["s_r"], dims="stock")
        r = pm.Deterministic("r", pt.exp(log_r), dims="stock")
        log_K = pm.Normal("log_K", k_mu, k_sd, dims="stock")
        K = pt.exp(log_K)
        if sigma_p == "pooled":
            sp = pm.LogNormal("sigma_p", sp_mu, sp_sd) * pt.ones(d["S"])
        elif sigma_p == "stock":
            sp = pm.LogNormal("sigma_p", sp_mu, sp_sd, dims="stock")
        else:
            sp = pt.as_tensor(np.full(d["S"], float(sigma_p)))
        logP = pm.Flat("logP", dims="state")
        Pp = pt.exp(logP)[prev]
        f = Pp + r[st[rest]] * Pp * (1 - Pp) - d["C"][prev] / K[st[rest]]
        f_pos = eps * pt.softplus(f / eps)
        pm.Potential("start", pm.logp(pm.Normal.dist(np.log(0.8), 0.3), logP[first]).sum())
        pm.Potential("process", pm.logp(pm.Normal.dist(pt.log(f_pos), sp[st[rest]]),
                                        logP[rest]).sum())
        pm.Potential("feasible", -pt.softplus((floor - f) / pen).sum())
        if end_prior is not None:
            lo, hi = np.log(end_prior[0]), np.log(end_prior[1])
            lpl = logP[d["last"]]  # log-density falls by 20 per unit of log B/K outside the range
            pm.Potential("end_prior", (-pt.softplus((lpl - hi) / 0.05)
                                       - pt.softplus((lo - lpl) / 0.05)).sum())
        if d["ser_name"]:
            lq = pm.Normal("log_qK", d["ser_m"] - np.log(0.5), 2.0, dims="series")
            sigma_o = pm.LogNormal("sigma_o", np.log(0.25), 0.5, dims="series")
            pm.Normal("logI", lq[d["obs_ser"]] + logP[d["obs_state"]], sigma_o[d["obs_ser"]],
                      observed=d["logI"])
    return m


def sample(m, d, seed, ta=0.9, **kw):
    t0 = time.time()
    idata = pm.sample(model=m, random_seed=seed, target_accept=ta, progressbar=False,
                      initvals={"logP": np.full(d["N"], np.log(0.5))}, **kw)
    return idata, time.time() - t0


def draws(idata, v):
    """Posterior draws as an array (draw, ...), chains stacked."""
    return idata.posterior[v].stack(s=("chain", "draw")).transpose("s", ...).values


def status(idata, d):
    """B/B_MSY = 2P and F/F_MSY = (C/B)/(r/2) for every stock-year and draw."""
    P = np.exp(draws(idata, "logP"))
    r, K = draws(idata, "r")[:, d["st"]], np.exp(draws(idata, "log_K"))[:, d["st"]]
    return 2 * P, d["C"][None, :] / (K * P) / (r / 2)


def report(idata, d, names, label, secs):
    div = int(idata.sample_stats["diverging"].sum())
    sm = pd.concat([az.summary(idata, var_names=[v], round_to=4) for v in names])
    print(f"{label}: {secs:.0f} s, {div} divergences of {idata.posterior.sizes['draw'] * 4:,}, "
          f"max r_hat {sm.r_hat.max():.3f}, min bulk ESS {sm.ess_bulk.min():.0f}, min tail ESS "
          f"{sm.ess_tail.min():.0f} (warm-up {idata.posterior.attrs.get('tuning_steps')})")
    return div


D1 = prep(["PCODHS"], multi_only=False)
print(f"{D1['N']} states, {len(D1['logI'])} index values in {len(D1['ser_name'])} series; "
      f"prior r 95% range {np.exp(D1['m_r'][0] - 2 * D1['s_r'][0]):.2f}-"
      f"{np.exp(D1['m_r'][0] + 2 * D1['s_r'][0]):.2f}")

69 states, 68 index values in 2 series; prior r 95% range 0.05-0.50

#The prior, checked by simulation

Because the dynamics are potentials, PyMC cannot draw a prior predictive for us, so we simulate it in NumPy: draw $r$, $K$, $\sigma_p$ and the start from their priors, run Pacific cod's real catches through the dynamics, and look at the biomass paths. A sensible prior should allow stocks from barely touched to heavily fished, and should not make most stocks collapse under the catches that were actually taken.

def prior_paths(d, n, gen, eps=0.01, k_mu=np.log(5.0), k_sd=0.8):
    r = np.exp(gen.normal(d["m_r"][0], d["s_r"][0], n))
    K = np.exp(gen.normal(k_mu, k_sd, n))
    sp = np.exp(gen.normal(np.log(0.1), 0.5, n))
    P = np.empty((n, d["N"]))
    P[:, 0] = np.exp(gen.normal(np.log(0.8), 0.3, n))
    fmin = np.full(n, np.inf)
    for t in range(1, d["N"]):
        f = P[:, t - 1] + r * P[:, t - 1] * (1 - P[:, t - 1]) - d["C"][t - 1] / K
        fmin = np.minimum(fmin, f)
        P[:, t] = eps * np.logaddexp(0, f / eps) * np.exp(sp * gen.normal(size=n))
    return 2 * P, K, fmin


BB_PRIOR, K_PRIOR, FMIN = prior_paths(D1, 4000, np.random.default_rng(7))
crash = FMIN < 0.02
print(f"K in tonnes, 89% prior range: {np.quantile(K_PRIOR, 0.055) * D1['cmax'][0]:,.0f}-"
      f"{np.quantile(K_PRIOR, 0.945) * D1['cmax'][0]:,.0f} (largest catch "
      f"{D1['cmax'][0]:,.0f} t)")
print(f"prior paths in which the catch would empty the stock: {crash.mean():.0%} "
      f"(the feasibility penalty removes them)")
print(f"surviving paths: final B/B_MSY median {np.median(BB_PRIOR[~crash, -1]):.2f}, "
      f"89% {np.quantile(BB_PRIOR[~crash, -1], 0.055):.2f}-"
      f"{np.quantile(BB_PRIOR[~crash, -1], 0.945):.2f}")
fig, ax = plt.subplots(figsize=(10, 4))
yrs1 = D1["years"]
for i in np.where(~crash)[0][:40]:
    ax.plot(yrs1, BB_PRIOR[i], color=HI, alpha=0.25, lw=0.9)
for i in np.where(crash)[0][:10]:
    ax.plot(yrs1, BB_PRIOR[i], color=LO, alpha=0.5, lw=0.9)
ax.axhline(1, color=INK, lw=1)
ax.set_ylim(0, 4)
ax.set_ylabel("B / B$_{MSY}$")
ax.set_title(wrap(f"Prior paths for Pacific cod under its real catches: 40 that survive (blue) "
                  f"and 10 of the {crash.mean():.0%} in which the catch empties the stock "
                  f"(vermilion)", 95), loc="left");

K in tonnes, 89% prior range: 13,138-162,694 (largest catch 9,540 t)
prior paths in which the catch would empty the stock: 71% (the feasibility penalty removes them)
surviving paths: final B/B_MSY median 1.83, 89% 0.93-2.48

The prior allows $K$ from about 1.4 to 18 times the largest catch (13,000 to 163,000 t in this simulation), yet about 70% of the prior draws cannot support the catches that were taken: with a growth rate near 0.16, a stock five times the largest catch has an MSY of only a fifth of it, and Pacific cod was fished near its largest catch for years.

The prior that reaches the data is the surviving 30%, and it has an opinion: those paths mostly end above the healthy level (final B/B_MSY median about 1.8, 89% range 0.9-2.5), because a stock big enough to survive the catches was never heavily depleted by them. An index can overrule that opinion. Without one, in part F, it is the answer.

#Failure 1: non-centred states

The usual advice for a hierarchical model with a poorly identified scale is to non-centre: sample standardised innovations $e_t \sim N(0, 1)$ and build the states as $\log P_t = \log \tilde f_t + \sigma_p e_t$ in a scan. We try that first, with a short run.

def nc_model(d, eps=0.01):
    """Non-centred single-stock version: innovations through a scan recursion."""
    T = d["N"]
    with pm.Model() as m:
        log_r = pm.Normal("log_r", d["m_r"][0], d["s_r"][0])
        r = pm.Deterministic("r", pt.exp(log_r))
        log_K = pm.Normal("log_K", np.log(5.0), 0.8)
        K = pt.exp(log_K)
        sp = pm.LogNormal("sigma_p", np.log(0.1), 0.5)
        e0 = pm.Normal("e0", 0, 1)
        e = pm.Normal("e", 0, 1, shape=T - 1)

        def step(e_t, c_t, lp, r, K, sp):
            p = pt.exp(lp)
            f = p + r * p * (1 - p) - c_t / K
            return pt.log(eps * pt.softplus(f / eps)) + sp * e_t, f

        lps, fs = pytensor.scan(step, sequences=[e, pt.as_tensor(d["C"][:-1])],
                                outputs_info=[np.log(0.8) + 0.3 * e0, None],
                                non_sequences=[r, K, sp], return_updates=False)
        logP = pm.Deterministic("logP", pt.concatenate([(np.log(0.8) + 0.3 * e0)[None], lps]))
        pm.Potential("feasible", -pt.softplus((0.02 - fs) / 0.005).sum())
        lq = pm.Normal("log_qK", d["ser_m"] - np.log(0.5), 2.0, shape=len(d["ser_name"]))
        so = pm.LogNormal("sigma_o", np.log(0.25), 0.5, shape=len(d["ser_name"]))
        pm.Normal("logI", lq[d["obs_ser"]] + logP[d["obs_state"]], so[d["obs_ser"]],
                  observed=d["logI"])
    return m


t0 = time.time()
id_nc = pm.sample(model=nc_model(D1), random_seed=RANDOM_SEED, target_accept=0.95, tune=400,
                  draws=500, progressbar=False)
DIV_NC = report(id_nc, D1, ["log_r", "log_K", "sigma_p"], "non-centred, Pacific cod",
                time.time() - t0)
TIMES["non-centred"] = time.time() - t0
del id_nc

non-centred, Pacific cod: 44 s, 194 divergences of 2,000, max r_hat 1.556, min bulk ESS 7, min tail ESS 7 (warm-up 400)

About one transition in ten diverges, r_hat reaches 1.6 and the worst bulk ESS is 7: the chains have not mixed. The reason is the one E46 met with a chaotic map. Every state depends on every earlier innovation through the nonlinear dynamics, so a small change in $r$ or an early $e_t$ moves the whole path, while the index pins each $\log P_t$ tightly: a curved, narrow tube in the innovations that no single step size can follow. Non-centring helps when the data are weak about each state; here they are strong.

#The centred fit

id1, secs = sample(ss_model(D1), D1, RANDOM_SEED, ta=0.95)
DIV_C1 = report(id1, D1, ["log_r", "log_K", "sigma_p", "sigma_o", "logP"],
                "centred, Pacific cod", secs)
TIMES["centred"] = secs
dv = id1.sample_stats["diverging"].stack(s=("chain", "draw")).values
sp1 = draws(id1, "sigma_p").ravel()
print(f"sigma_p: median {np.median(sp1):.3f}; in divergent draws "
      f"{np.median(sp1[dv]) if dv.any() else float('nan'):.3f}")
print(f"r: prior median {np.exp(D1['m_r'][0]):.2f}, posterior median "
      f"{np.median(draws(id1, 'r')):.2f}")

centred, Pacific cod: 3 s, 28 divergences of 4,000, max r_hat 1.071, min bulk ESS 40, min tail ESS 55 (warm-up 400)
sigma_p: median 0.400; in divergent draws 0.523
r: prior median 0.16, posterior median 0.15

With free states the same stock samples in a few seconds, and the worst r_hat drops from 1.6 to about 1.07. That is better, not clean: a few dozen divergences remain, the lowest ESS is about 40, and the divergent draws have a larger process noise than the rest (median about 0.52 against 0.40). This is not the familiar funnel at a scale near zero. When $\sigma_p$ is large, the 12 years before the first index (1944-1955) are free to wander, and the posterior gets a thin, wide-spreading region that the sampler cannot enter smoothly. We will come back to it.

The fit should first reproduce what it was fitted to. The posterior predictive for each index is the fitted state times the catchability with fresh observation noise.

lp1, lq1, so1 = draws(id1, "logP"), draws(id1, "log_qK"), draws(id1, "sigma_o")
pp_rng = np.random.default_rng(3)
fig, axs = plt.subplots(1, 2, figsize=(12, 3.8))
for i, (a, nm) in enumerate(zip(axs, D1["ser_name"])):
    sel = D1["obs_ser"] == i
    stt = D1["obs_state"][sel]
    rep = lq1[:, [i]] + lp1[:, stt] + so1[:, [i]] * pp_rng.normal(size=(len(lq1), sel.sum()))
    lo_, hi_ = np.quantile(rep, [0.055, 0.945], axis=0)
    a.fill_between(D1["years"][stt], np.exp(lo_), np.exp(hi_), color=HI, alpha=0.25, lw=0,
                   label="89% predictive")
    a.plot(D1["years"][stt], np.exp(D1["logI"][sel]), "o", color=INK, ms=3.5, label="observed")
    inside = ((D1["logI"][sel] >= lo_) & (D1["logI"][sel] <= hi_)).mean()
    a.set_yscale("log")
    a.legend(fontsize=9.5)
    a.set_title(f"{nm.split(':')[1]}: {inside:.0%} of {sel.sum()} values inside the 89% band",
                loc="left")
    print(f"{nm}: {inside:.0%} inside; sigma_o median {np.median(so1[:, i]):.2f}")

PCODHS:cpue: 100% inside; sigma_o median 0.40
PCODHS:survey: 91% inside; sigma_o median 0.72

All 57 catch-rate values and 10 of the 11 survey values sit inside their 89% predictive bands, if anything too good. The catch rate gets an observation sd of about 0.4 and the short survey about 0.7: the survey disagrees with the catch rate in the late 1980s and the model trusts the longer series. This check cannot fail very badly in a state-space model, because the states can follow the index wherever the process noise allows.

The data barely move $r$ (posterior median 0.15, prior 0.16): like the toy's one-way trip, the index says how the stock went down and partly back up, not how fast it grows.

#Failure 2: a long record with no index

Georges Bank yellowtail flounder has 80 years of catch (from 1935) but a survey only from 1987, so for half a century the states are tied only to the catches and the dynamics. Fitted alone with its own process noise:

DG = prep(["YELLGB"], multi_only=False)
idg, secs = sample(ss_model(DG), DG, RANDOM_SEED, ta=0.95)
DIV_G = report(idg, DG, ["log_r", "log_K", "sigma_p", "logP"], "Georges Bank yellowtail alone",
               secs)
dv = idg.sample_stats["diverging"].stack(s=("chain", "draw")).values
spg = draws(idg, "sigma_p").ravel()
print(f"sigma_p: median {np.median(spg):.3f}, 89% {np.quantile(spg, 0.055):.3f}-"
      f"{np.quantile(spg, 0.945):.3f}; in divergent draws {np.median(spg[dv]):.3f}")
TIMES["yellowtail alone"] = secs
del idg

Georges Bank yellowtail alone: 5 s, 298 divergences of 4,000, max r_hat 1.024, min bulk ESS 270, min tail ESS 313 (warm-up 400)
sigma_p: median 0.540, 89% 0.377-0.737; in divergent draws 0.682

About 300 divergences. The process sd comes out at about 0.54 (89% 0.38-0.74), more than for Pacific cod, and the divergent draws again sit at its upper end (median about 0.68), where five decades of index-free states can roam. One stock cannot tell how big its own process noise is when most of its states are unobserved; part D lets the stocks share it.

#D · Eighteen stocks, shared biology

#What the stocks share

Two things are pooled. First, growth. The prior for each stock's $\log r$ is centred at the midpoint of its resilience range, and stocks of the same family share a family effect:

$$\log r_s = m_s + \sigma_{fam}\, w_{f(s)} + \tau_r\, z_s, \qquad w, z \sim N(0, 1),\quad \sigma_{fam}, \tau_r \sim \mathrm{HalfNormal}(0.4).$$

Only cod, flatfish and hake have more than one stock here. A family effect for a family with a single stock cannot be told apart from that stock's own effect: the data see only their sum, and in our trial fits that unidentified pair added divergences. So the six single-stock families get no family term.

Second, process noise. Good and bad years differ between stocks, but their typical size may be similar enough for one shared $\sigma_p$, which lets the stocks with long index records inform those with long gaps (part C showed that a stock's own noise is poorly learned there). We fit the 18 stocks both ways, first with one $\sigma_p$ per stock (a short run of 250 draws per chain), then with one shared value. The model has 1,030 latent states, 21 index series with 728 values, and about 80 other parameters.

MAIN = ["PCODHS", "PCODWCVANI", "PCOD5AB", "RSOLE5AB", "RSOLEHSTR", "YELL3LNO", "YELLGB",
        "WITFLOUN4RST", "SPHAKECH", "SHAKE4VWX", "SHAKESGBMATL", "SHAKEGOMNGB", "REDFISHSPP3LN",
        "NZSNAPNZ7", "TREVALLYTRE7", "NZLINGESE", "BLSHARNPAC", "WHAKE4T"]
D = prep(MAIN)
print(f"{D['S']} stocks, {D['N']:,} states, {len(D['ser_name'])} series, "
      f"{len(D['logI'])} index values; families with a shared effect: "
      f"{[f for f in D['fam_names'] if f != '(single)']}")
MAIN_VARS = ["r", "log_K", "sigma_p", "sigma_o", "log_qK", "sigma_fam", "tau_r", "w_fam", "z_r",
             "logP"]
id_ps, secs = sample(ss_model(D, hier=True, sigma_p="stock"), D, RANDOM_SEED, ta=0.9, draws=250)
TIMES["18 stocks, per-stock noise"] = secs
DIV_PS = report(id_ps, D, MAIN_VARS, "18 stocks, per-stock process noise", secs)
dep_ps = id_ps.sample_stats["depth"].values
dv_ps = id_ps.sample_stats["diverging"].stack(s=("chain", "draw")).values
sp_ps = draws(id_ps, "sigma_p")
jg = MAIN.index("YELLGB")
print(f"tree depth 10 in {(dep_ps >= 10).mean():.0%} of draws; Georges Bank yellowtail sigma_p "
      f"median {np.median(sp_ps[:, jg]):.2f}, in divergent draws "
      f"{np.median(sp_ps[dv_ps, jg]):.2f}; per-stock medians {np.median(sp_ps, 0).min():.2f}-"
      f"{np.median(sp_ps, 0).max():.2f}")
del id_ps

18 stocks, 1,030 states, 21 series, 728 index values; families with a shared effect: ['Gadidae', 'Merlucciidae', 'Pleuronectidae']

18 stocks, per-stock process noise: 53 s, 114 divergences of 1,000, max r_hat 1.131, min bulk ESS 30, min tail ESS 40 (warm-up 400)
tree depth 10 in 55% of draws; Georges Bank yellowtail sigma_p median 0.54, in divergent draws 0.62; per-stock medians 0.08-0.54

Per-stock noise fails: more than one draw in ten diverges, r_hat reaches about 1.13, and over half the transitions hit nutpie's maximum tree depth of 10. The stocks' own noise estimates range from about 0.1 to 0.55, and Georges Bank yellowtail again sits at the top, with its divergent draws higher still. That spread suggests the noise really does differ between stocks, so a shared value trades some bias for a model that samples; its prior, centred at 0.1, also sits well below where the stocks will put it. Now the shared value:

M_MAIN = ss_model(D, hier=True, sigma_p="pooled")
id_main, secs = sample(M_MAIN, D, RANDOM_SEED, ta=0.9)
TIMES["18 stocks"] = secs
DIV_MAIN = report(id_main, D, MAIN_VARS, "18 stocks, shared process noise", secs)
SP = draws(id_main, "sigma_p").ravel()
print(f"shared process sd {np.median(SP):.3f} (89% {np.quantile(SP, 0.055):.3f}-"
      f"{np.quantile(SP, 0.945):.3f}); sigma_fam {np.median(draws(id_main, 'sigma_fam')):.2f}, "
      f"tau_r {np.median(draws(id_main, 'tau_r')):.2f}")
dep = id_main.sample_stats["depth"].values
print(f"tree depth: mean {dep.mean():.1f}, max {dep.max()}, at depth 10 in "
      f"{int((dep >= 10).sum())} draws ({(dep >= 10).mean():.1%})")
dv_m = id_main.sample_stats["diverging"].stack(s=("chain", "draw")).values
r_m = draws(id_main, "r")
top10 = (r_m[dv_m] > np.quantile(r_m, 0.9, axis=0)).mean(axis=0)  # per stock, among divergent draws
for s in np.argsort(top10)[::-1][:3]:
    print(f"{NICE[MAIN[s]]}: {top10[s]:.0%} of the divergent draws in the top tenth of its r "
          f"(median r {np.median(r_m[dv_m, s]):.2f} in them, {np.median(r_m[:, s]):.2f} overall)")

18 stocks, shared process noise: 48 s, 11 divergences of 4,000, max r_hat 1.016, min bulk ESS 308, min tail ESS 312 (warm-up 400)
shared process sd 0.236 (89% 0.211-0.264); sigma_fam 0.43, tau_r 0.30
tree depth: mean 7.2, max 10, at depth 10 in 1 draws (0.0%)
Trevally, New Zealand (TRE 7): 64% of the divergent draws in the top tenth of its r (median r 0.36 in them, 0.17 overall)
Ling, SE Australia: 45% of the divergent draws in the top tenth of its r (median r 0.23 in them, 0.16 overall)
Rock sole, Queen Charlotte Sound: 36% of the divergent draws in the top tenth of its r (median r 0.18 in them, 0.17 overall)

About ten divergences out of 4,000 draws remain, against 114 in 1,000 with per-stock noise. r_hat is at most about 1.02, every parameter and state has at least about 300 effective draws in the bulk and in the tails, a single draw reached the maximum tree depth of 10 (more than half did with per-stock noise), and the fit takes under a minute. The divergent draws sit in the upper tail of trevally's growth rate: about two thirds of them are in its top tenth, with a median $r$ of about 0.36 against 0.17 overall. That is few enough to report and move on, not few enough to ignore: the posterior has a corner the sampler visits reluctantly.

The shared process sd comes out at about 0.24: a typical stock departs from the Schaefer prediction by about a quarter each year. Part of that is real (recruitment varies by factors of several), and part is everything a logistic curve leaves out: age structure, changing selectivity, a changing environment.

#Does the 18-stock fit reproduce its indices?

The same predictive check as for Pacific cod, now for all 21 index series: the share of each series' values inside its 89% posterior predictive band, and the worst series drawn out.

LP_M, LQ_M, SO_M = draws(id_main, "logP"), draws(id_main, "log_qK"), draws(id_main, "sigma_o")
ppc_rng = np.random.default_rng(5)
sub = ppc_rng.choice(LP_M.shape[0], 1000, replace=False)
os_, oi_ = D["obs_ser"], D["obs_state"]
REP = (LQ_M[sub][:, os_] + LP_M[sub][:, oi_]
       + SO_M[sub][:, os_] * ppc_rng.normal(size=(len(sub), len(os_))))
lo_m, hi_m = np.quantile(REP, [0.055, 0.945], axis=0)
inside_m = (D["logI"] >= lo_m) & (D["logI"] <= hi_m)
PPC = pd.DataFrame({"series": [f"{NICE[n.split(':')[0]]}: {n.split(':')[1]}"
                               for n in D["ser_name"]],
                    "values": np.bincount(os_),
                    "inside_89": np.bincount(os_, inside_m) / np.bincount(os_),
                    "sigma_o": np.median(SO_M, 0)}).set_index("series").sort_values("inside_89")
print(PPC.round(2).to_string())
print(f"\nall series: {inside_m.mean():.0%} of {len(inside_m)} values inside their 89% band; "
      f"series below 80%: {int((PPC.inside_89 < 0.8).sum())}")
worst = int(np.argmin(np.bincount(os_, inside_m) / np.bincount(os_)))
w_sel = os_ == worst
fig, ax = plt.subplots(figsize=(10, 3.6))
ax.fill_between(D["years"][oi_[w_sel]], np.exp(lo_m[w_sel]), np.exp(hi_m[w_sel]), color=HI,
                alpha=0.25, lw=0, label="89% predictive")
ax.plot(D["years"][oi_[w_sel]], np.exp(np.median(REP[:, w_sel], 0)), color=HI, lw=1.2,
        label="predictive median")
ax.plot(D["years"][oi_[w_sel]], np.exp(D["logI"][w_sel]), "o", color=INK, ms=3.5, label="observed")
ax.set_yscale("log")
ax.legend(fontsize=9.5)
ax.set_title(wrap(f"Worst-covered series, {PPC.index[0]}: {PPC.inside_89.iloc[0]:.0%} of "
                  f"{PPC['values'].iloc[0]} values inside the 89% band", 70), loc="left");
del LP_M, REP

                                                 values  inside_89  sigma_o
series                                                                     
Yellowtail flounder, Georges Bank: survey            28       0.86     0.83
Pacific cod, Hecate Strait: survey                   11       0.91     0.63
Pacific cod, Hecate Strait: cpue                     57       0.93     0.55
White hake, S. Gulf of St Lawrence: survey           39       0.95     0.39
Silver hake, Gulf of Maine: cpue                     47       0.96     0.36
Pacific cod, W. Vancouver Island: cpue               46       0.98     0.33
Pacific cod, Queen Charlotte Sound: cpue             57       0.98     0.33
Rock sole, Queen Charlotte Sound: cpue               47       1.00     0.10
Rock sole, Hecate Strait: survey                     12       1.00     0.35
Rock sole, Hecate Strait: cpue                       52       1.00     0.13
Witch flounder, Gulf of St Lawrence: survey          41       1.00     0.38
Yellowtail flounder, Grand Banks: cpue               45       1.00     0.08
Hake, Chile: cpue                                    31       1.00     0.16
Silver hake, Scotian Shelf: survey                   38       1.00     0.38
Silver hake, S. Georges Bank-Mid-Atlantic: cpue      47       1.00     0.29
Rock sole, Queen Charlotte Sound: survey              7       1.00     0.31
Redfish, Grand Banks: cpue                           36       1.00     0.11
Snapper, New Zealand (SNA 7): cpue                   25       1.00     0.31
Ling, SE Australia: cpue                             21       1.00     0.10
Trevally, New Zealand (TRE 7): cpue                  23       1.00     0.12
Blue shark, North Pacific: cpue                      18       1.00     0.10

all series: 98% of 728 values inside their 89% band; series below 80%: 0

Of the 728 index values, 98% fall inside their 89% bands, more than the nominal 89% because the states are free to follow each index, the same leniency as for Pacific cod. No series falls below 80%. The worst, Georges Bank yellowtail's survey (86%), has the largest observation sd (about 0.8). Two isolated spikes, 107 and 115 kg per tow in 2008-09 against 21 in 2007 and 15 in 2010, sit above the band; they look like outliers, and they inflate that sd. The last two values, 2013-14, about a hundredth of the spikes, sit well below it: so noisy a survey lets the model treat the final collapse partly as noise. Part E tests both. The check cannot say whether the dynamics are right; the official assessments test that below.

#How much did the data teach us about growth?

For each stock we compare the prior and posterior of $\log r$. The prior here is the hierarchical one, simulated in NumPy from the same half-normal scales. The contraction is $1 - \mathrm{var}_{post}/\mathrm{var}_{prior}$: 0 means the data taught nothing, 1 that they settled it.

DRAW = {v: draws(id_main, v) for v in ("r", "log_K")}
BB, UU = status(id_main, D)
pr_rng = np.random.default_rng(11)
n_pr = 20000
sf_pr, tr_pr = np.abs(pr_rng.normal(0, 0.4, n_pr)), np.abs(pr_rng.normal(0, 0.4, n_pr))
LOGR_PRIOR = {}
rows = []
for s, sid in enumerate(MAIN):
    single = D["fam_names"][D["fam"][s]] == "(single)"
    pl = (D["m_r"][s] + (0 if single else sf_pr * pr_rng.normal(size=n_pr))
          + tr_pr * pr_rng.normal(size=n_pr))
    LOGR_PRIOR[sid] = pl
    po = np.log(DRAW["r"][:, s])
    msy = np.exp(po + DRAW["log_K"][:, s]) / 4 * D["cmax"][s]
    rows.append(dict(stock=NICE[sid], resilience=RAM.loc[RAM.stockid == sid, "resilience"].iloc[0],
                     cmsy_range=f"{np.exp(D['m_r'][s] - 2 * D['s_r'][s]):.2g}-"
                                f"{np.exp(D['m_r'][s] + 2 * D['s_r'][s]):.2g}",
                     r_prior=np.exp(np.median(pl)),
                     prior_95=f"{np.exp(np.quantile(pl, 0.025)):.2f}-"
                              f"{np.exp(np.quantile(pl, 0.975)):.2f}",
                     r_post=np.exp(np.median(po)),
                     r_post_89=f"{np.exp(np.quantile(po, 0.055)):.2f}-"
                               f"{np.exp(np.quantile(po, 0.945)):.2f}",
                     contraction=1 - po.var() / pl.var(),
                     MSY_t=f"{np.median(msy):,.0f} ({np.quantile(msy, 0.055):,.0f}-"
                           f"{np.quantile(msy, 0.945):,.0f})"))
RTAB = pd.DataFrame(rows).set_index("stock")
print(RTAB.round(2).to_string())
print(f"contraction: median {RTAB.contraction.median():.0%}, range {RTAB.contraction.min():.0%} "
      f"to {RTAB.contraction.max():.0%}")

                                             resilience cmsy_range  r_prior   prior_95  r_post  r_post_89  contraction                     MSY_t
stock                                                                                                                                           
Pacific cod, Hecate Strait                          Low   0.05-0.5     0.16  0.05-0.53    0.14  0.06-0.23         0.46       2,735 (1,271-5,732)
Pacific cod, W. Vancouver Island                    Low   0.05-0.5     0.16  0.05-0.55    0.13  0.06-0.23         0.43         1,089 (513-2,072)
Pacific cod, Queen Charlotte Sound                  Low   0.05-0.5     0.16  0.05-0.51    0.13  0.06-0.22         0.47           979 (429-2,492)
Rock sole, Queen Charlotte Sound                    Low   0.05-0.5     0.16  0.05-0.51    0.17  0.09-0.31         0.51           623 (319-1,512)
Rock sole, Hecate Strait                            Low   0.05-0.5     0.16  0.05-0.51    0.16  0.09-0.27         0.58         1,210 (652-2,139)
Yellowtail flounder, Grand Banks                    Low   0.05-0.5     0.16  0.05-0.52    0.17  0.09-0.30         0.51     17,701 (9,134-34,142)
Yellowtail flounder, Georges Bank                   Low   0.05-0.5     0.16  0.05-0.53    0.15  0.08-0.27         0.51     11,835 (5,796-24,289)
Witch flounder, Gulf of St Lawrence              Medium    0.2-0.8     0.40  0.12-1.33    0.32  0.15-0.54         0.48      5,661 (2,103-18,392)
Hake, Chile                                      Medium    0.2-0.8     0.40  0.12-1.34    0.18  0.09-0.31         0.47   68,827 (35,441-142,292)
Silver hake, Scotian Shelf                       Medium    0.2-0.8     0.40  0.12-1.36    0.21  0.10-0.39         0.39   54,658 (22,249-153,613)
Silver hake, S. Georges Bank-Mid-Atlantic        Medium    0.2-0.8     0.40  0.12-1.29    0.16  0.07-0.29         0.37  159,357 (66,753-529,143)
Silver hake, Gulf of Maine                       Medium    0.2-0.8     0.40  0.12-1.32    0.24  0.11-0.48         0.32    40,340 (19,641-95,729)
Redfish, Grand Banks                       Very low-Low  0.015-0.5     0.09  0.04-0.21    0.09  0.06-0.19         0.06     18,468 (8,366-51,056)
Snapper, New Zealand (SNA 7)                        Low   0.05-0.5     0.16  0.06-0.38    0.15  0.09-0.21         0.48         1,193 (603-3,011)
Trevally, New Zealand (TRE 7)                       Low   0.05-0.5     0.16  0.06-0.38    0.17  0.12-0.43        -0.00         1,745 (920-4,043)
Ling, SE Australia                                  Low   0.05-0.5     0.16  0.07-0.38    0.16  0.10-0.33         0.10           693 (323-1,583)
Blue shark, North Pacific                           Low   0.05-0.5     0.16  0.07-0.38    0.16  0.10-0.33         0.14   67,862 (33,411-150,838)
White hake, S. Gulf of St Lawrence                  Low   0.05-0.5     0.16  0.07-0.38    0.16  0.09-0.25         0.33      6,989 (2,877-18,605)
contraction: median 45%, range -0% to 58%

In the hierarchy the CMSY range sets only the centre of each stock's prior; the spread comes from the two half-normal scales, so the prior 95% ranges are narrower than CMSY's for stocks alone in their family, wider for the hakes, and much narrower for redfish (0.04-0.21 against 0.015-0.5). The redfish's small contraction therefore means the data added little to a range the hierarchy had already narrowed.

For most stocks the data narrowed $r$ by a third to a half (median about 45%) and left its centre close to the prior's; for trevally, ling and blue shark, each alone in its family, by about 15% or less. The hakes show the family effect working: their "medium" prior is centred at 0.4, the four stocks together say about 0.2, and the family effect carries that lesson to all four. So the growth rate is largely prior, adjusted a little by the data, the state of knowledge for one-way trips. MSY's 89% ranges span factors of about three to nine.

#A check against the official assessments

Fourteen of the eighteen stocks have an official B/B_MSY series in RAM. For each, we compare the shape of the two histories (the correlation of the log posterior median with the log official value, across years), the share of years in which the official value lies inside our 89% interval, and whether the two agree on the status in the official series' last year (our P(B < B_MSY) above or below one half, against the official value below or above 1).

rows = []
for s, sid in enumerate(MAIN):
    sel = D["st"] == s
    g = stock_rows(sid)
    off = g.bdivbmsy.to_numpy()
    ok = ~np.isnan(off)
    if ok.sum() < 6:
        continue
    b = BB[:, sel]
    med = np.median(b, 0)
    lo_, hi_ = np.quantile(b[:, ok], [0.055, 0.945], axis=0)
    j = np.where(ok)[0][-1]
    p_below = (b[:, j] < 1).mean()
    rows.append(dict(stock=NICE[sid], years=int(ok.sum()), last=int(g.year[j]),
                     corr_log=np.corrcoef(np.log(med[ok]), np.log(off[ok]))[0, 1],
                     inside_89=((off[ok] >= lo_) & (off[ok] <= hi_)).mean(),
                     ours_median=np.median(b[:, j]), P_below=p_below,
                     P_half=(b[:, j] < 0.5).mean(), official=off[j],
                     agree=(p_below > 0.5) == (off[j] < 1)))
AGREE = pd.DataFrame(rows).set_index("stock").sort_values("corr_log", ascending=False)
print(AGREE.round(2).assign(corr_log=AGREE.corr_log.round(3)).to_string())
print(f"\n{len(AGREE)} stocks: median correlation {AGREE.corr_log.median():.2f}; "
      f"{int((AGREE.corr_log >= 0.8).sum())} at 0.8 or above; status agrees for "
      f"{int(AGREE.agree.sum())} of {len(AGREE)}; official inside our 89% band in "
      f"{AGREE.inside_89.mean():.0%} of stock-years on average")
OFF_LOW, OURS_LOW = AGREE.official < 1, AGREE.P_below > 0.5
print(pd.crosstab(OURS_LOW.rename("ours: P(B<Bmsy) > 0.5"),
                  OFF_LOW.rename("official B/Bmsy < 1")).to_string())
N_OPT, N_PES = int((OFF_LOW & ~OURS_LOW).sum()), int((~OFF_LOW & OURS_LOW).sum())
print(f"officially below B_MSY: {int(OFF_LOW.sum())}, of which we call "
      f"{int((OFF_LOW & OURS_LOW).sum())}; "
      f"disagreements: {N_OPT} with our model more optimistic, {N_PES} more pessimistic")
OVERF = AGREE[AGREE.official < 0.5]
print(f"officially overfished (below half B_MSY) in their last assessed year: {len(OVERF)}; "
      f"ours P(B < Bmsy/2) there: "
      + ", ".join(f"{nm} {int(r['last'])} {r['P_half']:.2f}" for nm, r in OVERF.iterrows()))