{
"cells": [
{
"cell_type": "raw",
"id": "056b8bfc",
"metadata": {},
"source": [
"---\n",
"format: html\n",
"title: \"MCMC for a logistic fit from early data: the notebook\"\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "0d5937e8",
"metadata": {},
"source": [
"This is the notebook behind [MCMC for a logistic fit from early data](../../logistic-mcmc.html), which in turn\n",
"follows [Fisher information for a logistic fit from early data](../../logistic.html). The post shows the\n",
"figures and says what they mean. This notebook is the methodology: what is simulated, what is fitted, how the\n",
"sampler is set up and checked, and where the numbers quoted in the post come from. It is meant to be read top to\n",
"bottom.\n",
"\n",
"To re-run it, `pixi run notebook` in `src/link/logistic_mcmc/` of the\n",
"[blog's repository](https://github.com/ickc/blog.kolen.dev). It takes about a minute. Every random number is\n",
"seeded, so a re-run reproduces this page exactly."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "e05e902c",
"metadata": {},
"outputs": [],
"source": [
"import emcee\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pandas as pd\n",
"import seaborn as sns\n",
"from IPython.display import HTML\n",
"from matplotlib.ticker import FuncFormatter\n",
"from scipy.optimize import minimize\n",
"\n",
"%config InlineBackend.figure_formats = ['svg']\n",
"%config InlineBackend.print_figure_kwargs = {\"bbox_inches\": \"tight\", \"metadata\": {\"Date\": None}}"
]
},
{
"cell_type": "markdown",
"id": "9db67089",
"metadata": {},
"source": [
"Figures are SVG in the site's body font, TeX Gyre Schola. matplotlib's default `svg.fonttype = \"path\"` turns\n",
"text into outlines, so a reader does not need the font installed. Anything with thousands of marks is drawn with\n",
"`rasterized=True`, which keeps the axes and text vector but the file small."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "38075d30",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"INK = \"#333333\"\n",
"MUTED = \"#777777\"\n",
"GRID = \"#e6e6e6\"\n",
"# One colour per cutoff in the figures that show three cutoffs: blue, orange, green, in that order everywhere.\n",
"COLOURS = [\"#2a78d6\", \"#eb6834\", \"#1baf7a\"]\n",
"POSTERIOR = \"#2a78d6\"\n",
"\n",
"plt.rcParams.update({\n",
" \"font.family\": \"serif\",\n",
" \"font.serif\": [\"TeX Gyre Schola\", \"DejaVu Serif\"],\n",
" \"mathtext.fontset\": \"stix\",\n",
" \"font.size\": 10,\n",
" \"text.color\": INK,\n",
" \"axes.labelcolor\": INK,\n",
" \"axes.edgecolor\": MUTED,\n",
" \"axes.linewidth\": 0.6,\n",
" \"xtick.color\": MUTED,\n",
" \"ytick.color\": MUTED,\n",
" \"xtick.labelcolor\": INK,\n",
" \"ytick.labelcolor\": INK,\n",
" \"xtick.direction\": \"out\",\n",
" \"ytick.direction\": \"out\",\n",
" \"axes.spines.top\": False,\n",
" \"axes.spines.right\": False,\n",
" \"axes.grid\": True,\n",
" \"grid.color\": GRID,\n",
" \"grid.linewidth\": 0.5,\n",
" \"axes.axisbelow\": True,\n",
" \"legend.frameon\": False,\n",
" \"figure.facecolor\": \"white\",\n",
" \"axes.facecolor\": \"white\",\n",
" \"savefig.facecolor\": \"white\",\n",
" \"svg.fonttype\": \"path\",\n",
" \"svg.hashsalt\": \"logistic-mcmc\", # with Date: None above, a re-run gives byte-identical SVG\n",
"})"
]
},
{
"cell_type": "markdown",
"id": "d902acac",
"metadata": {},
"source": [
"## The model\n",
"\n",
"Same model and notation as the previous post,\n",
"\n",
"$$\n",
"y(t) = \\frac{L}{1 + e^{-k(t - t_0)}} = \\frac{L A e^{kt}}{L + A e^{kt}}, \\qquad A = L e^{-k t_0},\n",
"$$\n",
"\n",
"observed with multiplicative noise: we see $\\log y(t_i) + \\epsilon_i$ with $\\epsilon_i \\sim N(0, \\sigma^2)$.\n",
"\n",
"The parameters are $\\theta = (\\log A, k, \\log L)$ rather than Cook's $(L, k, t_0)$, because they split into\n",
"what early data can pin down, $A$ and $k$ (the exponential $A e^{kt}$ that the left tail looks like), and what\n",
"they cannot, $L$. The inflection point is then derived, $t_0 = (\\log L - \\log A) / k$. Logs are used for $A$\n",
"and $L$ because both are scales, and because the uncertainty on $L$ turns out to span decades.\n",
"\n",
"In these parameters\n",
"\n",
"$$\n",
"\\log y = \\log L - \\log\\left(1 + \\frac{L}{A} e^{-kt}\\right),\n",
"$$\n",
"\n",
"and the second term is a softplus, `np.logaddexp(0, x)`, which stays accurate for any $t$. That matters: as\n",
"$L \\to \\infty$ at fixed $A$ and $k$ the curve becomes exactly the exponential, and the code should be able to\n",
"evaluate that limit without overflowing."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "f91facd5",
"metadata": {},
"outputs": [],
"source": [
"def log_y(t, log_a, k, log_l):\n",
" \"\"\"log y(t) in the (log A, k, log L) parametrization, stable for any t and any L.\"\"\"\n",
" return log_l - np.logaddexp(0.0, log_l - log_a - k * t)"
]
},
{
"cell_type": "markdown",
"id": "d48b3daf",
"metadata": {},
"source": [
"## The simulated data\n",
"\n",
"Units are chosen so that $L = 1$, $k = 1$ and $t_0 = 0$. That is no loss of generality: it is the same as\n",
"measuring $y$ as a fraction of its limit and time in $e$-folding times from the inflection point. The noise is\n",
"10%, $\\sigma = 0.1$, and there are $\\rho = 10$ observations per $e$-folding time, starting at $t = -10$, where\n",
"$y/L \\approx 5 \\times 10^{-5}$ and the curve is exponential to any precision the noise could see.\n",
"\n",
"One noise realization is drawn and then cut off at successive times $T$, as if watching the data come in. The\n",
"cutoffs are named by how far up to the limit the data have reached, $y(T)/L$, which is the quantity the Fisher\n",
"information in the previous post depends on."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "64e6356d",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
T
\n",
"
observations
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
-4.60
\n",
"
55
\n",
"
\n",
"
\n",
"
2%
\n",
"
-3.89
\n",
"
62
\n",
"
\n",
"
\n",
"
5%
\n",
"
-2.94
\n",
"
71
\n",
"
\n",
"
\n",
"
10%
\n",
"
-2.20
\n",
"
79
\n",
"
\n",
"
\n",
"
20%
\n",
"
-1.39
\n",
"
87
\n",
"
\n",
"
\n",
"
50%
\n",
"
0.00
\n",
"
101
\n",
"
\n",
"
\n",
"
80%
\n",
"
1.39
\n",
"
114
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" T observations\n",
"y(T)/L \n",
"1% -4.60 55\n",
"2% -3.89 62\n",
"5% -2.94 71\n",
"10% -2.20 79\n",
"20% -1.39 87\n",
"50% 0.00 101\n",
"80% 1.39 114"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"L_TRUE, K_TRUE, T0_TRUE = 1.0, 1.0, 0.0\n",
"LOG_A_TRUE = np.log(L_TRUE) - K_TRUE * T0_TRUE\n",
"THETA_TRUE = np.array([LOG_A_TRUE, K_TRUE, np.log(L_TRUE)])\n",
"SIGMA, RHO, T_START = 0.1, 10.0, -10.0\n",
"\n",
"REACHED = np.array([0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 0.8])\n",
"T_CUT = T0_TRUE + np.log(REACHED / (1 - REACHED)) / K_TRUE\n",
"SHOWN = [0.01, 0.1, 0.5] # the three cutoffs drawn in the bands and pair figures\n",
"\n",
"rng = np.random.default_rng(20260926)\n",
"t_all = np.arange(T_START, T_CUT.max() + 1e-9, 1 / RHO)\n",
"obs_all = log_y(t_all, *THETA_TRUE) + rng.normal(0.0, SIGMA, t_all.size)\n",
"\n",
"\n",
"def data_up_to(T):\n",
" keep = t_all <= T\n",
" return t_all[keep], obs_all[keep]\n",
"\n",
"\n",
"pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"T\": T_CUT.round(2),\n",
" \"observations\": [(t_all <= T).sum() for T in T_CUT],\n",
"}).set_index(\"y(T)/L\")"
]
},
{
"cell_type": "markdown",
"id": "5c8d38da",
"metadata": {},
"source": [
"## The likelihood\n",
"\n",
"With Gaussian noise in $\\log y$ and $\\sigma$ known, as in the Fisher information calculation, the\n",
"log-likelihood is a sum of squares,\n",
"\n",
"$$\n",
"\\log \\mathcal{L}(\\theta) = -\\frac{1}{2\\sigma^2} \\sum_i \\left( \\log y_i^{\\text{obs}} - \\log y(t_i; \\theta)\n",
"\\right)^2 + \\text{const}.\n",
"$$\n",
"\n",
"Everything below, the maximum-likelihood fit, the profile likelihood and the posterior, is built on this one\n",
"function."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "c692e115",
"metadata": {},
"outputs": [],
"source": [
"def neg_log_like(theta, t, obs):\n",
" log_a, k, log_l = theta\n",
" r = obs - log_y(t, log_a, k, log_l)\n",
" return 0.5 * np.sum(r * r) / SIGMA**2"
]
},
{
"cell_type": "markdown",
"id": "5577b237",
"metadata": {},
"source": [
"## Maximum likelihood and the profile likelihood\n",
"\n",
"The key fact about this likelihood is what happens as $L \\to \\infty$ with $A$ and $k$ held fixed: the logistic\n",
"becomes the exponential $A e^{kt}$, so the likelihood tends to a constant, the likelihood of the best\n",
"exponential fit. It does not go to zero. If the data cannot tell a logistic from an exponential, the\n",
"exponential can even be the best fit, and then the maximum-likelihood estimate (MLE) of $L$ is infinite.\n",
"\n",
"The clean way to see this is the *profile likelihood* of $\\log L$: fix $\\log L$, maximize over the other two\n",
"parameters, and plot the result against $\\log L$. Plotted as\n",
"$-2 \\Delta \\log \\mathcal{L} = 2 \\left( \\min_{A,k} (-\\log \\mathcal{L}) - \\min_\\theta (-\\log \\mathcal{L}) \\right)$,\n",
"its minimum is at the MLE, and by Wilks' theorem the set where it is below $\\chi^2_{1, 0.95} = 3.84$ is an\n",
"approximate 95% confidence interval. No prior is involved anywhere.\n",
"\n",
"The profile is computed on a grid in $\\log L$ out to $e^{20}$, far past anything physical, so that the plateau\n",
"to the exponential is visible. Each grid point is a two-parameter Nelder–Mead fit, warm-started from its\n",
"neighbour; the grid is swept in both directions and the better fit kept, which guards against a warm start\n",
"getting stuck.\n",
"\n",
"The same is done for $t_0$, with $\\log L = \\log A + k t_0$ substituted, since $t_0$ is not one of the sampled\n",
"parameters and its profile has to be computed on its own."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e0122ed3",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
MLE of L
\n",
"
95% interval for L
\n",
"
−2Δ log-likelihood of the exponential
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
∞
\n",
"
[0.0563, ∞]
\n",
"
0.00
\n",
"
\n",
"
\n",
"
2%
\n",
"
0.577
\n",
"
[0.0958, ∞]
\n",
"
0.16
\n",
"
\n",
"
\n",
"
5%
\n",
"
1.42
\n",
"
[0.252, ∞]
\n",
"
0.17
\n",
"
\n",
"
\n",
"
10%
\n",
"
0.577
\n",
"
[0.314, 3.33]
\n",
"
5.62
\n",
"
\n",
"
\n",
"
20%
\n",
"
1.22
\n",
"
[0.681, 5.23]
\n",
"
6.53
\n",
"
\n",
"
\n",
"
50%
\n",
"
1
\n",
"
[0.834, 1.23]
\n",
"
117.27
\n",
"
\n",
"
\n",
"
80%
\n",
"
1.11
\n",
"
[0.987, 1.19]
\n",
"
624.33
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" MLE of L 95% interval for L −2Δ log-likelihood of the exponential\n",
"y(T)/L \n",
"1% ∞ [0.0563, ∞] 0.00\n",
"2% 0.577 [0.0958, ∞] 0.16\n",
"5% 1.42 [0.252, ∞] 0.17\n",
"10% 0.577 [0.314, 3.33] 5.62\n",
"20% 1.22 [0.681, 5.23] 6.53\n",
"50% 1 [0.834, 1.23] 117.27\n",
"80% 1.11 [0.987, 1.19] 624.33"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"LOG_L_GRID = np.arange(-5.0, 20.0 + 1e-9, 0.05)\n",
"T0_GRID = np.arange(-8.0, 25.0 + 1e-9, 0.05)\n",
"NM = {\"xatol\": 1e-9, \"fatol\": 1e-11, \"maxiter\": 5_000}\n",
"CHI2_95 = 3.841\n",
"\n",
"\n",
"def profile(t, obs, grid, of=\"log_l\"):\n",
" \"\"\"min over (log A, k) of -log L, at each value of log L (of=\"log_l\") or t0 (of=\"t0\").\"\"\"\n",
" def theta(p, v):\n",
" return (p[0], p[1], v) if of == \"log_l\" else (p[0], p[1], p[0] + p[1] * v)\n",
"\n",
" best = np.full(grid.size, np.inf)\n",
" for order in (np.arange(grid.size), np.arange(grid.size)[::-1]):\n",
" x = THETA_TRUE[:2]\n",
" for i in order:\n",
" fit = minimize(lambda p: neg_log_like(theta(p, grid[i]), t, obs), x, method=\"Nelder-Mead\", options=NM)\n",
" x = fit.x\n",
" best[i] = min(best[i], fit.fun)\n",
" return best\n",
"\n",
"\n",
"def mle_and_interval(grid, prof):\n",
" \"\"\"MLE and the 95% profile-likelihood interval, with np.inf where the profile never climbs back up.\"\"\"\n",
" d = 2 * (prof - prof.min())\n",
" # On the plateau the profile is flat to optimizer precision, so \"the minimum is at the far end of the grid\"\n",
" # has to allow for noise at the level of the tolerances.\n",
" mle = np.inf if d[-1] < 1e-4 else grid[np.argmin(prof)]\n",
" inside = np.flatnonzero(d <= CHI2_95)\n",
" lo_i, hi_i = inside[0], inside[-1]\n",
" lo = np.interp(CHI2_95, [d[lo_i], d[lo_i - 1]], [grid[lo_i], grid[lo_i - 1]]) if lo_i > 0 else -np.inf\n",
" hi = np.interp(CHI2_95, [d[hi_i], d[hi_i + 1]], [grid[hi_i], grid[hi_i + 1]]) if hi_i < grid.size - 1 else np.inf\n",
" return mle, lo, hi\n",
"\n",
"\n",
"profiles, profiles_t0, mle, mle_t0, profile_int, profile_int_t0 = {}, {}, {}, {}, {}, {}\n",
"for u, T in zip(REACHED, T_CUT):\n",
" t, obs = data_up_to(T)\n",
" profiles[u] = profile(t, obs, LOG_L_GRID)\n",
" profiles_t0[u] = profile(t, obs, T0_GRID, of=\"t0\")\n",
" mle[u], *profile_int[u] = mle_and_interval(LOG_L_GRID, profiles[u])\n",
" mle_t0[u], *profile_int_t0[u] = mle_and_interval(T0_GRID, profiles_t0[u])\n",
"\n",
"\n",
"def fmt(x):\n",
" return \"∞\" if np.isinf(x) else f\"{0 if abs(x) < 5e-4 else x:.3g}\"\n",
"\n",
"\n",
"pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"MLE of L\": [fmt(np.exp(mle[u])) for u in REACHED],\n",
" \"95% interval for L\": [f\"[{fmt(np.exp(profile_int[u][0]))}, {fmt(np.exp(profile_int[u][1]))}]\" for u in REACHED],\n",
" \"−2Δ log-likelihood of the exponential\": [\n",
" f\"{2 * (profiles[u][-1] - profiles[u].min()):.2f}\" for u in REACHED],\n",
"}).set_index(\"y(T)/L\")"
]
},
{
"cell_type": "markdown",
"id": "ccb485d8",
"metadata": {},
"source": [
"At 1% the profile never climbs back up: it is lowest on the plateau, so the best fit is the exponential and\n",
"the MLE of $L$ is infinite. At 2% and 5% the MLE is finite, but the exponential is only 0.2 worse in $-2 \\log \\mathcal{L}$, far\n",
"inside the 3.84 threshold, so the interval is still unbounded above. From 10% the exponential is ruled out, only\n",
"just at 10%, and the interval closes.\n",
"\n",
"The last column is the height of the plateau. It is worth keeping in mind for the posterior below: a plateau\n",
"that is excluded at 95% can still hold a lot of posterior mass if it is long enough."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "08e2bab8",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| label: fig-profile\n",
"#| fig-cap: \"Profile likelihood of $L$. Below the dotted line is the 95% interval, and the curves flatten out to the right at the height of the best exponential fit.\"\n",
"fig, ax = plt.subplots(figsize=(6.5, 3.8), layout=\"constrained\")\n",
"palette = sns.color_palette(\"crest\", len(REACHED))\n",
"for u, colour in zip(REACHED, palette):\n",
" ax.plot(LOG_L_GRID / np.log(10), 2 * (profiles[u] - profiles[u].min()), color=colour, lw=1.6, label=f\"{u:.0%}\")\n",
"ax.axhline(CHI2_95, color=MUTED, lw=0.8, ls=\":\")\n",
"ax.text(8.6, CHI2_95 + 0.3, \"95%\", color=MUTED, fontsize=8, ha=\"right\")\n",
"ax.axvline(0, color=INK, lw=0.6, ls=\"--\")\n",
"ax.set_xlim(-1.6, 8.7)\n",
"ax.set_ylim(0, 15)\n",
"ax.xaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"$10^{{{v:g}}}$\"))\n",
"ax.set_xlabel(r\"$L$ (true value 1)\")\n",
"ax.set_ylabel(r\"$-2\\,\\Delta \\log \\mathcal{L}$, profiled over $A$, $k$\")\n",
"ax.legend(title=\"data up to $y(T)/L$\", fontsize=8, title_fontsize=8, loc=\"upper right\", ncols=2)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "ddd2ce65",
"metadata": {},
"source": [
"## The prior and the posterior\n",
"\n",
"The same plateau is why the posterior needs care. With a flat prior on $\\log L$ over the whole real line, the\n",
"posterior density tends to a positive constant as $\\log L \\to \\infty$ and cannot be normalized. A sampler would\n",
"simply walk off to infinity. So the prior on $\\log L$ has to stop somewhere, and here it is flat between\n",
"$e^{-7}$ and $e^{7}$, about $1/1000$ to $1000$ times the true value. The priors on $\\log A$ and on $k > 0$ are\n",
"flat and wide enough never to matter.\n",
"\n",
"The upper bound on $L$ is not a technicality. Wherever the likelihood is on its plateau, the posterior mass\n",
"there is proportional to the plateau's height *times its length*, and the length is set by the prior. So:\n",
"\n",
"- At 1–5%, where the plateau is as high as the peak, the posterior is spread up to the bound, and its upper\n",
" quantiles are the bound's.\n",
"- At 10–20% the plateau is excluded by the likelihood, but it is 14 units of $\\log L$ long, so it still holds a\n",
" substantial share of the posterior. The credible interval reaches up into it even though the profile\n",
" interval does not.\n",
"\n",
"Move the bound and those upper ends move with it."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1e73a279",
"metadata": {},
"outputs": [],
"source": [
"LOG_L_PRIOR = (-7.0, 7.0)\n",
"\n",
"\n",
"def log_posterior(theta, t, obs):\n",
" log_a, k, log_l = theta\n",
" if not (0.0 < k < 10.0 and LOG_L_PRIOR[0] < log_l < LOG_L_PRIOR[1] and -50.0 < log_a < 50.0):\n",
" return -np.inf\n",
" return -neg_log_like(theta, t, obs)"
]
},
{
"cell_type": "markdown",
"id": "97902807",
"metadata": {},
"source": [
"## Sampling with emcee\n",
"\n",
"The posterior is sampled with [emcee](https://emcee.readthedocs.io/), an ensemble sampler: 48 walkers move\n",
"together, and each proposal for one walker is built from the positions of the others, so the sampler adapts to\n",
"the scale and correlations of the posterior without tuning.\n",
"\n",
"Choices, and why:\n",
"\n",
"- **Start.** The walkers start in a tiny ball around the true parameters. That is the most favourable start\n",
" possible, and it is deliberate: the question here is what the data allow, not whether a sampler can find the\n",
" answer from a bad start.\n",
"- **Burn-in.** The first 2,000 of 12,000 steps are dropped. The walkers start bunched together, and those steps\n",
" reflect the start, not the posterior. How long that takes is checked below with the autocorrelation time.\n",
"- **Moves.** emcee's default is the stretch move. Here it is replaced by differential-evolution moves\n",
" (`DEMove`, with some `DESnookerMove`), because the posterior in $\\log L$ is a long flat plateau for the early\n",
" cutoffs, which the stretch move explores slowly. That claim is checked below rather than taken on trust.\n",
"- **Thinning.** The chain is thinned by half the autocorrelation time, which keeps the samples nearly\n",
" independent and the arrays small."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "5261b4a7",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
autocorrelation time, steps
\n",
"
burn-in / autocorrelation time
\n",
"
effective samples
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
16.0
\n",
"
125.0
\n",
"
29900.0
\n",
"
\n",
"
\n",
"
2%
\n",
"
19.0
\n",
"
108.0
\n",
"
25800.0
\n",
"
\n",
"
\n",
"
5%
\n",
"
16.0
\n",
"
126.0
\n",
"
30200.0
\n",
"
\n",
"
\n",
"
10%
\n",
"
43.0
\n",
"
47.0
\n",
"
11200.0
\n",
"
\n",
"
\n",
"
20%
\n",
"
44.0
\n",
"
46.0
\n",
"
10900.0
\n",
"
\n",
"
\n",
"
50%
\n",
"
13.0
\n",
"
157.0
\n",
"
37600.0
\n",
"
\n",
"
\n",
"
80%
\n",
"
13.0
\n",
"
153.0
\n",
"
36600.0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" autocorrelation time, steps burn-in / autocorrelation time \\\n",
"y(T)/L \n",
"1% 16.0 125.0 \n",
"2% 19.0 108.0 \n",
"5% 16.0 126.0 \n",
"10% 43.0 47.0 \n",
"20% 44.0 46.0 \n",
"50% 13.0 157.0 \n",
"80% 13.0 153.0 \n",
"\n",
" effective samples \n",
"y(T)/L \n",
"1% 29900.0 \n",
"2% 25800.0 \n",
"5% 30200.0 \n",
"10% 11200.0 \n",
"20% 10900.0 \n",
"50% 37600.0 \n",
"80% 36600.0 "
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N_WALKERS, N_STEPS, N_BURN = 48, 12_000, 2_000\n",
"DE_MOVES = [(emcee.moves.DEMove(), 0.8), (emcee.moves.DESnookerMove(), 0.2)]\n",
"\n",
"\n",
"def run_sampler(T, seed, moves=DE_MOVES):\n",
" t, obs = data_up_to(T)\n",
" p0 = THETA_TRUE + 1e-3 * np.random.default_rng(seed).normal(size=(N_WALKERS, 3))\n",
" sampler = emcee.EnsembleSampler(N_WALKERS, 3, log_posterior, args=(t, obs), moves=moves)\n",
" sampler.random_state = np.random.RandomState(seed).get_state()\n",
" sampler.run_mcmc(p0, N_STEPS)\n",
" return sampler\n",
"\n",
"\n",
"samples, taus = {}, {}\n",
"for i, (u, T) in enumerate(zip(REACHED, T_CUT)):\n",
" sampler = run_sampler(T, seed=i)\n",
" taus[u] = sampler.get_autocorr_time(discard=N_BURN, quiet=True)\n",
" log_a, k, log_l = sampler.get_chain(discard=N_BURN, thin=max(1, int(taus[u].max() / 2)), flat=True).T\n",
" samples[u] = pd.DataFrame({\"log_a\": log_a, \"k\": k, \"log_l\": log_l, \"t0\": (log_l - log_a) / k})\n",
"\n",
"pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"autocorrelation time, steps\": [taus[u].max().round(0) for u in REACHED],\n",
" \"burn-in / autocorrelation time\": [(N_BURN / taus[u].max()).round(0) for u in REACHED],\n",
" \"effective samples\": [(N_WALKERS * (N_STEPS - N_BURN) / taus[u].max()).round(-2) for u in REACHED],\n",
"}).set_index(\"y(T)/L\")"
]
},
{
"cell_type": "markdown",
"id": "29a49caf",
"metadata": {},
"source": [
"The autocorrelation time is the number of steps a walker takes to forget where it was. The usual rule of thumb\n",
"is that the burn-in should be several of them and the chain at least 50; here the burn-in is at least 45 of\n",
"them and the chain over 200, and every cutoff has more than ten thousand effective samples.\n",
"\n",
"The comparison with the default stretch move, on the two cutoffs where the plateau matters most:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "9a1fc488",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
\n",
"
autocorrelation time of log L
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
moves
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
stretch (default)
\n",
"
52.0
\n",
"
\n",
"
\n",
"
differential evolution
\n",
"
14.0
\n",
"
\n",
"
\n",
"
10%
\n",
"
stretch (default)
\n",
"
121.0
\n",
"
\n",
"
\n",
"
differential evolution
\n",
"
41.0
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" autocorrelation time of log L\n",
"y(T)/L moves \n",
"1% stretch (default) 52.0\n",
" differential evolution 14.0\n",
"10% stretch (default) 121.0\n",
" differential evolution 41.0"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"move_check = []\n",
"for u in (0.01, 0.1):\n",
" T = T_CUT[REACHED == u][0]\n",
" for name, moves in ((\"stretch (default)\", None), (\"differential evolution\", DE_MOVES)):\n",
" tau = run_sampler(T, seed=0, moves=moves).get_autocorr_time(discard=N_BURN, quiet=True)\n",
" move_check.append({\"y(T)/L\": f\"{u:.0%}\", \"moves\": name, \"autocorrelation time of log L\": tau[2].round(0)})\n",
"pd.DataFrame(move_check).set_index([\"y(T)/L\", \"moves\"])"
]
},
{
"cell_type": "markdown",
"id": "b0210fe0",
"metadata": {},
"source": [
"The stretch move would work, but it needs about three times as many steps for the same number of effective\n",
"samples."
]
},
{
"cell_type": "markdown",
"id": "db75fc68",
"metadata": {},
"source": [
"## Which point estimate\n",
"\n",
"Three candidates for \"the\" value of $L$:\n",
"\n",
"- **Posterior median.** Invariant under monotone reparametrization, but where the data are weak it is simply\n",
" the middle of whatever prior range is left, and moves when the prior bound moves.\n",
"- **MAP**, the maximum of the posterior density. With flat priors in $(\\log A, k, \\log L)$ it is the MLE clipped\n",
" to the prior box, so at 1% it lands on the bound. It also depends on the parametrization: a density in $L$ is\n",
" the density in $\\log L$ divided by $L$, so the peak in one is not the peak in the other.\n",
"- **MLE.** Needs no prior and is invariant under reparametrization: the MLE of $L$ is $e$ to the MLE of\n",
" $\\log L$, and $t_0$ at the MLE is the MLE of $t_0$.\n",
"\n",
"There is no real prior information in this exercise, and a flat prior is not \"no information\" either (flat in\n",
"$\\log L$, flat in $L$ and flat in $y(T)/L$ give different answers precisely where the data are weak). So the MLE\n",
"is used as the point estimate, with the profile-likelihood interval that goes with it. The posterior is used\n",
"for what it is good at, showing the whole shape.\n",
"\n",
"For comparison:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "1b5ec551",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
MLE
\n",
"
MAP (MLE within the prior)
\n",
"
posterior median
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
∞
\n",
"
1.1e+03
\n",
"
11.4
\n",
"
\n",
"
\n",
"
2%
\n",
"
0.577
\n",
"
0.577
\n",
"
11
\n",
"
\n",
"
\n",
"
5%
\n",
"
1.42
\n",
"
1.42
\n",
"
17.5
\n",
"
\n",
"
\n",
"
10%
\n",
"
0.577
\n",
"
0.577
\n",
"
0.838
\n",
"
\n",
"
\n",
"
20%
\n",
"
1.22
\n",
"
1.22
\n",
"
1.54
\n",
"
\n",
"
\n",
"
50%
\n",
"
1
\n",
"
1
\n",
"
1
\n",
"
\n",
"
\n",
"
80%
\n",
"
1.11
\n",
"
1.11
\n",
"
1.08
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" MLE MAP (MLE within the prior) posterior median\n",
"y(T)/L \n",
"1% ∞ 1.1e+03 11.4\n",
"2% 0.577 0.577 11\n",
"5% 1.42 1.42 17.5\n",
"10% 0.577 0.577 0.838\n",
"20% 1.22 1.22 1.54\n",
"50% 1 1 1\n",
"80% 1.11 1.11 1.08"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"MLE\": [fmt(np.exp(mle[u])) for u in REACHED],\n",
" \"MAP (MLE within the prior)\": [fmt(np.exp(min(mle[u], LOG_L_PRIOR[1]))) for u in REACHED],\n",
" \"posterior median\": [f\"{np.exp(samples[u].log_l.median()):.3g}\" for u in REACHED],\n",
"}).set_index(\"y(T)/L\")"
]
},
{
"cell_type": "markdown",
"id": "0b6cf368",
"metadata": {},
"source": [
"At 1–5% the three disagree by orders of magnitude, which is itself the answer: no point estimate means much\n",
"before the data reach about 10% of the limit.\n",
"\n",
"For drawing, the MLE curve needs $A$ and $k$ too, from a full three-parameter fit started at the best grid\n",
"point. At 1% the MLE is the exponential, which the parametrization handles by setting $\\log L$ very large."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "bf86bc90",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"mle_theta = {}\n",
"for u, T in zip(REACHED, T_CUT):\n",
" t, obs = data_up_to(T)\n",
" i = np.argmin(profiles[u])\n",
" start = np.r_[THETA_TRUE[:2], LOG_L_GRID[i]]\n",
" if np.isinf(mle[u]):\n",
" fit = minimize(lambda p: neg_log_like((p[0], p[1], 50.0), t, obs), THETA_TRUE[:2], method=\"Nelder-Mead\", options=NM)\n",
" mle_theta[u] = np.r_[fit.x, 50.0]\n",
" else:\n",
" mle_theta[u] = minimize(neg_log_like, start, args=(t, obs), method=\"Nelder-Mead\", options=NM).x"
]
},
{
"cell_type": "markdown",
"id": "5c6f1337",
"metadata": {},
"source": [
"## What the data allow\n",
"\n",
"The MLE curve and the 50% and 95% posterior bands of $y(t)$, for three cutoffs. The dots are the data, and\n",
"everything right of the last one is extrapolation."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "6aaad913",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| label: fig-bands\n",
"#| fig-cap: \"The maximum-likelihood fit, and the 50% and 95% posterior bands of $y(t)$, with data up to 1%, 10% and 50% of the limit. The vertical line is the last observation.\"\n",
"def bands(df, t):\n",
" curves = np.exp(log_y(t[None, :], df.log_a.values[:, None], df.k.values[:, None], df.log_l.values[:, None]))\n",
" return np.percentile(curves, [2.5, 25, 75, 97.5], axis=0)\n",
"\n",
"\n",
"t_plot = np.linspace(-8, 8, 400)\n",
"truth = np.exp(log_y(t_plot, *THETA_TRUE))\n",
"\n",
"fig, axes = plt.subplots(1, 3, figsize=(9, 3.4), sharey=True, layout=\"constrained\")\n",
"for ax, u, colour in zip(axes, SHOWN, COLOURS):\n",
" T = T_CUT[REACHED == u][0]\n",
" lo95, lo50, hi50, hi95 = bands(samples[u].sample(4000, random_state=0), t_plot)\n",
" ax.fill_between(t_plot, lo95, hi95, color=colour, alpha=0.15, lw=0, label=\"95%\")\n",
" ax.fill_between(t_plot, lo50, hi50, color=colour, alpha=0.35, lw=0, label=\"50%\")\n",
" ax.plot(t_plot, np.exp(log_y(t_plot, *mle_theta[u])), color=colour, lw=1.5, label=\"MLE\")\n",
" ax.plot(t_plot, truth, color=INK, lw=0.8, ls=\"--\", label=\"truth\")\n",
" t, obs = data_up_to(T)\n",
" keep = t >= t_plot[0]\n",
" ax.plot(t[keep], np.exp(obs[keep]), \"o\", ms=2.2, color=INK, alpha=0.7, mec=\"none\", label=\"data\")\n",
" ax.axvline(T, color=MUTED, lw=0.6)\n",
" ax.axhline(np.exp(LOG_L_PRIOR[1]), color=MUTED, lw=0.8, ls=\":\", label=\"prior bound on $L$\")\n",
" ax.set_yscale(\"log\")\n",
" ax.set_ylim(3e-4, 3e3)\n",
" ax.set_xlim(t_plot[0], t_plot[-1])\n",
" ax.set_title(f\"data up to $y(T)/L$ = {u:.0%}\", fontsize=10, loc=\"left\")\n",
" ax.set_xlabel(r\"$k(t - t_0)$\")\n",
"axes[0].set_ylabel(r\"$y / L$\")\n",
"axes[0].yaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"{v:g}\"))\n",
"axes[-1].legend(loc=\"lower right\", fontsize=8, handlelength=1.5)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "d9be188d",
"metadata": {},
"source": [
"At 1% the MLE is the exponential itself and runs off the top. The posterior bands at 1% and 10% reach the prior\n",
"bound, and at 50% they have closed around the truth.\n",
"\n",
"## The joint distribution\n",
"\n",
"$A$, $k$, $t_0$ and $L$ against each other, for the same three cutoffs, with $A$ and $L$ on log axes. The dashed\n",
"lines are the true values."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "9805b840",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| label: fig-pairs\n",
"#| fig-cap: \"Posterior of $A$, $k$, $t_0$ and $L$ at the same three cutoffs. Dashed lines are the true values, and each histogram is scaled to its own peak.\"\n",
"LABELS = {\"log_a\": r\"$\\log_{10} A$\", \"k\": r\"$k$\", \"t0\": r\"$t_0$\", \"log_l\": r\"$\\log_{10} L$\"}\n",
"TRUE = {\"log_a\": LOG_A_TRUE / np.log(10), \"k\": K_TRUE, \"t0\": T0_TRUE, \"log_l\": np.log10(L_TRUE)}\n",
"\n",
"pair = pd.concat(\n",
" [samples[u].sample(3000, random_state=0).assign(reached=f\"{u:.0%}\") for u in SHOWN],\n",
" ignore_index=True,\n",
")\n",
"pair[\"log_a\"] /= np.log(10)\n",
"pair[\"log_l\"] /= np.log(10)\n",
"pair = pair[list(LABELS) + [\"reached\"]]\n",
"\n",
"g = sns.pairplot(\n",
" pair, hue=\"reached\", palette=COLOURS, corner=True, height=1.9,\n",
" plot_kws={\"s\": 3, \"alpha\": 0.25, \"linewidth\": 0, \"rasterized\": True},\n",
" diag_kind=\"hist\",\n",
")\n",
"# Each cutoff's histogram is redrawn scaled to its own peak: on a common density scale the 50% one is a spike\n",
"# that flattens the other two into the axis.\n",
"for ax, var in zip(g.diag_axes, LABELS):\n",
" for artist in [*ax.patches, *ax.collections]:\n",
" artist.remove()\n",
" bins = np.linspace(pair[var].min(), pair[var].max(), 60)\n",
" for (_, grp), colour in zip(pair.groupby(\"reached\", sort=False), COLOURS):\n",
" counts, edges = np.histogram(grp[var], bins=bins)\n",
" ax.stairs(counts / counts.max(), edges, color=colour, fill=True, alpha=0.25, lw=0)\n",
" ax.stairs(counts / counts.max(), edges, color=colour, lw=1.0)\n",
" ax.set_ylim(0, 1.1)\n",
" # The corner layout leaves the diagonal without a y label, so name each histogram in place.\n",
" ax.set_title(LABELS[var], fontsize=10, loc=\"left\", pad=2)\n",
"for i, row in enumerate(LABELS):\n",
" for j, col in enumerate(LABELS):\n",
" if j > i:\n",
" continue\n",
" ax = g.axes[i, j]\n",
" ax.axvline(TRUE[col], color=INK, lw=0.6, ls=\"--\")\n",
" if i != j:\n",
" ax.axhline(TRUE[row], color=INK, lw=0.6, ls=\"--\")\n",
" g.axes[-1, i].set_xlabel(LABELS[row])\n",
" g.axes[i, 0].set_ylabel(LABELS[row])\n",
"g.axes[0, 0].set_ylabel(\"\")\n",
"sns.move_legend(g, \"upper right\", bbox_to_anchor=(0.9, 0.9), title=\"data up to $y(T)/L$\", frameon=False,\n",
" markerscale=4)\n",
"for handle in g.legend.legend_handles:\n",
" handle.set_alpha(1)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "194cad9c",
"metadata": {},
"source": [
"$A$ and $k$ are pinned down to a few percent even at 1%; they are strongly correlated with each other, which is\n",
"the usual trade-off between the amplitude and the rate of an exponential fitted over a finite window. $t_0$ and\n",
"$\\log L$ lie on the line $\\log L = \\log A + k t_0$, with a width set by how well $A$ and $k$ are known, so they\n",
"carry the same information. At 1% the histogram of $\\log L$ is flat up to the prior bound, which is the prior\n",
"showing through. At 10% the peak is near the truth but the tail runs to the bound, the plateau of the previous\n",
"section.\n",
"\n",
"## The intervals, as the data come in\n",
"\n",
"One figure for all seven cutoffs, for $L$ and for how far ahead of the last observation the turning point is,\n",
"$t_0 - T$. The violins are the posterior, cut at its extremes. The black bars are the 95% profile-likelihood\n",
"intervals and the dots the MLE; an arrow means the interval does not close, and an MLE at infinity is marked\n",
"with ∞ at the top.\n",
"\n",
"$L$ is on a log axis, which is wide enough for the whole prior range, so nothing needs cropping. $t_0 - T$ is\n",
"on a linear axis, cropped at 12 $e$-folding times, and the few unbounded intervals are drawn as arrows."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "7cab258c",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| label: fig-intervals\n",
"#| fig-cap: \"Posterior (violins), maximum-likelihood estimate (dots) and 95% profile-likelihood interval (bars) of the limit $L$ and of how far ahead of the last observation the turning point is, as the data come in. Arrows are intervals that do not close; the dashed lines are the truth.\"\n",
"labels = [f\"{u:.0%}\" for u in REACHED]\n",
"x = np.arange(REACHED.size)\n",
"T0_TOP = 12.0\n",
"\n",
"fig, (ax_l, ax_t) = plt.subplots(1, 2, figsize=(7.6, 3.6), layout=\"constrained\")\n",
"\n",
"violin = {\"inner\": None, \"cut\": 0, \"density_norm\": \"width\", \"width\": 0.8, \"linewidth\": 0}\n",
"sns.violinplot(\n",
" data=pd.concat([pd.DataFrame({\"reached\": lab, \"v\": samples[u].log_l / np.log(10)}) for u, lab in zip(REACHED, labels)]),\n",
" x=\"reached\", y=\"v\", ax=ax_l, color=POSTERIOR, alpha=0.3, **violin,\n",
")\n",
"sns.violinplot(\n",
" data=pd.concat([pd.DataFrame({\"reached\": lab, \"v\": (samples[u].t0 - T).clip(upper=T0_TOP)})\n",
" for u, lab, T in zip(REACHED, labels, T_CUT)]),\n",
" x=\"reached\", y=\"v\", ax=ax_t, color=POSTERIOR, alpha=0.3, **violin,\n",
")\n",
"\n",
"L_TOP = np.log10(np.exp(LOG_L_PRIOR[1])) + 0.45\n",
"for xi, u, T in zip(x, REACHED, T_CUT):\n",
" # L, in log10 units\n",
" lo, hi = np.array(profile_int[u]) / np.log(10)\n",
" top = min(hi, L_TOP - 0.25)\n",
" ax_l.plot([xi, xi], [lo, top], color=INK, lw=2, solid_capstyle=\"butt\")\n",
" if np.isinf(hi):\n",
" ax_l.annotate(\"\", xy=(xi, L_TOP - 0.05), xytext=(xi, top - 0.01),\n",
" arrowprops={\"arrowstyle\": \"-|>\", \"color\": INK, \"lw\": 2, \"shrinkA\": 0, \"shrinkB\": 0})\n",
" if np.isinf(mle[u]):\n",
" ax_l.text(xi, L_TOP + 0.05, \"∞\", ha=\"center\", va=\"bottom\", fontsize=11, color=INK)\n",
" else:\n",
" ax_l.plot(xi, mle[u] / np.log(10), \"o\", ms=6, color=INK, mec=\"white\", mew=1.2, zorder=3)\n",
" # t0 - T\n",
" lo, hi = np.array(profile_int_t0[u]) - T\n",
" top = min(hi, T0_TOP - 0.6)\n",
" ax_t.plot([xi, xi], [lo, top], color=INK, lw=2, solid_capstyle=\"butt\")\n",
" if hi > T0_TOP - 0.6:\n",
" ax_t.annotate(\"\", xy=(xi, T0_TOP), xytext=(xi, top - 0.01),\n",
" arrowprops={\"arrowstyle\": \"-|>\", \"color\": INK, \"lw\": 2, \"shrinkA\": 0, \"shrinkB\": 0})\n",
" if np.isinf(mle_t0[u]) or mle_t0[u] - T > T0_TOP:\n",
" ax_t.text(xi, T0_TOP + 0.15, \"∞\" if np.isinf(mle_t0[u]) else f\"{mle_t0[u] - T:.0f}\",\n",
" ha=\"center\", va=\"bottom\", fontsize=11 if np.isinf(mle_t0[u]) else 8, color=INK)\n",
" else:\n",
" ax_t.plot(xi, mle_t0[u] - T, \"o\", ms=6, color=INK, mec=\"white\", mew=1.2, zorder=3)\n",
" ax_t.plot([xi - 0.4, xi + 0.4], [T0_TRUE - T] * 2, color=INK, lw=0.8, ls=\"--\")\n",
"\n",
"ax_l.axhline(0, color=INK, lw=0.8, ls=\"--\")\n",
"ax_l.axhline(LOG_L_PRIOR[1] / np.log(10), color=MUTED, lw=0.8, ls=\":\")\n",
"ax_l.text(REACHED.size - 0.45, LOG_L_PRIOR[1] / np.log(10) + 0.06, \"prior bound\", color=MUTED, fontsize=8,\n",
" ha=\"right\", va=\"bottom\")\n",
"ax_l.set_ylim(-1.6, L_TOP)\n",
"ax_l.yaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"{10 ** v:g}\"))\n",
"ax_l.set_ylabel(\"limit $L$ (true value 1)\")\n",
"ax_t.set_ylim(-2.5, T0_TOP)\n",
"ax_t.set_ylabel(r\"turning point ahead of the data, $k(t_0 - T)$\")\n",
"for ax in (ax_l, ax_t):\n",
" ax.set_xlabel(\"data up to $y(T)/L$\")\n",
" ax.set_xlim(-0.6, REACHED.size - 0.4)\n",
" ax.grid(axis=\"x\", visible=False)\n",
" for spine in (\"top\",):\n",
" ax.spines[spine].set_visible(False)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "7bc28e34",
"metadata": {},
"source": [
"Reading left to right, as the data come in:\n",
"\n",
"- Up to 5% the profile interval only has a lower end, and so does anything honest one can say about $L$: the\n",
" data rule out a limit close to where they are, because that would already show as bending, and nothing more.\n",
" The MLE jumps around, and at 1% it is infinite.\n",
"- At 10% and 20% the profile interval closes, a few times the truth at the top. The posterior still stretches up\n",
" to the prior bound, which is the plateau effect described above.\n",
"- From the inflection point on, everything agrees.\n",
"\n",
"The right panel is the same thing in terms of the turning point, since $t_0$ and $\\log L$ carry the same\n",
"information once $A$ and $k$ are known."
]
},
{
"cell_type": "markdown",
"id": "7eb80754",
"metadata": {},
"source": [
"## Against the Fisher information\n",
"\n",
"The previous post gave the Cramér–Rao bound on $\\log L$ with $k$ and $A$ known,\n",
"\n",
"$$\n",
"\\operatorname{sd}(\\log \\hat{L}) \\gtrsim \\left[ \\frac{\\rho}{k \\sigma^2} \\left( -\\log\\left(1 - \\frac{y(T)}{L}\\right)\n",
"- \\frac{y(T)}{L} \\right) \\right]^{-1/2} \\approx \\sigma \\sqrt{\\frac{2k}{\\rho}} \\, \\frac{L}{y(T)}.\n",
"$$\n",
"\n",
"With all three fitted the bound comes from the full $3 \\times 3$ Fisher matrix. With $u(t) = y(t)/L$ the\n",
"derivatives of $\\log y$ are\n",
"\n",
"$$\n",
"\\frac{\\partial \\log y}{\\partial \\log A} = 1 - u, \\qquad\n",
"\\frac{\\partial \\log y}{\\partial k} = t (1 - u), \\qquad\n",
"\\frac{\\partial \\log y}{\\partial \\log L} = u,\n",
"$$\n",
"\n",
"so $I = J^\\top J / \\sigma^2$ with those as the columns of $J$, and the bound is $\\sqrt{(I^{-1})_{\\log L, \\log L}}$.\n",
"Both are evaluated at the true parameters. For the curves, the observation times end exactly at each $T$, so the\n",
"curve is not stepped by the sampling grid; for the table, they are the times actually used.\n",
"\n",
"Two things from the simulation go against it: the posterior standard deviation of $\\log L$, and the\n",
"half-width of the profile-likelihood interval divided by 1.96, which is what a standard deviation would be if the\n",
"likelihood were Gaussian."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "d998b06d",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"\n",
"
\n",
" \n",
"
\n",
"
\n",
"
CR, A and k known
\n",
"
CR, closed form
\n",
"
CR, joint
\n",
"
profile half-width / 1.96
\n",
"
posterior sd
\n",
"
\n",
"
\n",
"
y(T)/L
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
4.265
\n",
"
4.457
\n",
"
9.782
\n",
"
inf
\n",
"
2.692
\n",
"
\n",
"
\n",
"
2%
\n",
"
2.133
\n",
"
2.221
\n",
"
4.423
\n",
"
inf
\n",
"
2.616
\n",
"
\n",
"
\n",
"
5%
\n",
"
0.885
\n",
"
0.879
\n",
"
1.653
\n",
"
inf
\n",
"
2.321
\n",
"
\n",
"
\n",
"
10%
\n",
"
0.414
\n",
"
0.432
\n",
"
0.714
\n",
"
0.603
\n",
"
1.831
\n",
"
\n",
"
\n",
"
20%
\n",
"
0.202
\n",
"
0.208
\n",
"
0.323
\n",
"
0.520
\n",
"
1.289
\n",
"
\n",
"
\n",
"
50%
\n",
"
0.070
\n",
"
0.072
\n",
"
0.097
\n",
"
0.098
\n",
"
0.097
\n",
"
\n",
"
\n",
"
80%
\n",
"
0.036
\n",
"
0.035
\n",
"
0.044
\n",
"
0.047
\n",
"
0.046
\n",
"
\n",
" \n",
"
\n",
"
"
],
"text/plain": [
" CR, A and k known CR, closed form CR, joint \\\n",
"y(T)/L \n",
"1% 4.265 4.457 9.782 \n",
"2% 2.133 2.221 4.423 \n",
"5% 0.885 0.879 1.653 \n",
"10% 0.414 0.432 0.714 \n",
"20% 0.202 0.208 0.323 \n",
"50% 0.070 0.072 0.097 \n",
"80% 0.036 0.035 0.044 \n",
"\n",
" profile half-width / 1.96 posterior sd \n",
"y(T)/L \n",
"1% inf 2.692 \n",
"2% inf 2.616 \n",
"5% inf 2.321 \n",
"10% 0.603 1.831 \n",
"20% 0.520 1.289 \n",
"50% 0.098 0.097 \n",
"80% 0.047 0.046 "
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def fisher(t):\n",
" \"\"\"Fisher information for (log A, k, log L) at the true parameters, on observation times t.\"\"\"\n",
" u = np.exp(log_y(t, *THETA_TRUE)) / L_TRUE\n",
" J = np.column_stack([1 - u, t * (1 - u), u])\n",
" return J.T @ J / SIGMA**2\n",
"\n",
"\n",
"def cr_joint(t):\n",
" return np.sqrt(np.linalg.inv(fisher(t))[2, 2])\n",
"\n",
"\n",
"def cr_known(t):\n",
" return 1 / np.sqrt(fisher(t)[2, 2])\n",
"\n",
"\n",
"def times_up_to(T):\n",
" return np.arange(T, T_START, -1 / RHO)[::-1]\n",
"\n",
"\n",
"def cr_known_closed(u):\n",
" return 1 / np.sqrt(RHO / (K_TRUE * SIGMA**2) * (-np.log1p(-u) - u))\n",
"\n",
"\n",
"profile_sd = {u: (profile_int[u][1] - profile_int[u][0]) / (2 * 1.96) for u in REACHED}\n",
"\n",
"summary = pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"CR, A and k known\": [cr_known(data_up_to(T)[0]) for T in T_CUT],\n",
" \"CR, closed form\": cr_known_closed(REACHED),\n",
" \"CR, joint\": [cr_joint(data_up_to(T)[0]) for T in T_CUT],\n",
" \"profile half-width / 1.96\": [profile_sd[u] for u in REACHED],\n",
" \"posterior sd\": [samples[u].log_l.std() for u in REACHED],\n",
"}).set_index(\"y(T)/L\")\n",
"summary.round(3)"
]
},
{
"cell_type": "markdown",
"id": "3bca748f",
"metadata": {},
"source": [
"- The closed form with $A$ and $k$ known is within 5% of the sum over the actual observation times.\n",
"- Fitting $A$ and $k$ jointly roughly doubles the bound, but it still goes as $L / y(T)$, the divergence the\n",
" previous post said remains.\n",
"- From the inflection point on, the posterior sd and the profile interval both agree with the joint bound.\n",
"- Before it they do not, and the reason is that neither the likelihood nor the posterior is anywhere near\n",
" Gaussian there, while the bound is a local, Gaussian statement. At 10–20% the profile interval is lopsided,\n",
" short below and long above, and the posterior adds the plateau on top, so its sd is 2.5–4 times the bound. At\n",
" 1–5% the profile interval is unbounded, and the posterior sd *drops below* the bound, which an unbiased\n",
" estimator could never do. That is the prior cutting the tail off: a flat prior of width 14 in $\\log L$ has an\n",
" sd of $14/\\sqrt{12} \\approx 4.0$, the most the posterior sd can be."
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "78bd80f6",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"
"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#| label: fig-fisher\n",
"#| fig-cap: \"Spread of $\\\\log L$ from the simulation against the Cramér–Rao bounds of the last post.\"\n",
"u_grid = np.logspace(-2.3, np.log10(0.9), 200)\n",
"T_grid = T0_TRUE + np.log(u_grid / (1 - u_grid)) / K_TRUE\n",
"PRIOR_SD = (LOG_L_PRIOR[1] - LOG_L_PRIOR[0]) / np.sqrt(12)\n",
"\n",
"fig, ax = plt.subplots(figsize=(6, 4), layout=\"constrained\")\n",
"ax.plot(u_grid, [cr_joint(times_up_to(T)) for T in T_grid], color=COLOURS[0], lw=1.8,\n",
" label=r\"Cramér–Rao, $A$, $k$, $L$ fitted\")\n",
"ax.plot(u_grid, [cr_known(times_up_to(T)) for T in T_grid], color=COLOURS[1], lw=1.8,\n",
" label=r\"Cramér–Rao, $A$, $k$ known\")\n",
"ax.plot(u_grid, SIGMA * np.sqrt(2 * K_TRUE / RHO) / u_grid, color=COLOURS[1], lw=0.9, ls=\":\",\n",
" label=r\"$\\sigma \\sqrt{2k/\\rho}\\; L / y(T)$\")\n",
"ax.axhline(PRIOR_SD, color=MUTED, lw=0.8, ls=\"--\")\n",
"ax.text(0.3, PRIOR_SD * 1.15, \"flat prior on $\\\\log L$\", color=MUTED, fontsize=8, ha=\"center\")\n",
"ax.plot(REACHED, summary[\"posterior sd\"], \"o\", ms=6, color=INK, mec=\"white\", mew=1.2,\n",
" label=\"MCMC posterior sd\", zorder=3)\n",
"finite = np.isfinite(summary[\"profile half-width / 1.96\"])\n",
"ax.plot(REACHED[finite], summary[\"profile half-width / 1.96\"][finite], \"D\", ms=5, mfc=\"white\", mec=INK, mew=1.2,\n",
" label=\"profile interval half-width / 1.96\", zorder=3)\n",
"for u in REACHED[~finite]:\n",
" ax.annotate(\"\", xy=(u, 28), xytext=(u, 14),\n",
" arrowprops={\"arrowstyle\": \"-|>\", \"color\": INK, \"lw\": 1.2, \"shrinkA\": 0, \"shrinkB\": 0})\n",
"ax.set_xscale(\"log\")\n",
"ax.set_yscale(\"log\")\n",
"ax.set_ylim(top=30)\n",
"ax.set_xlabel(r\"how far up the data have reached, $y(T)/L$\")\n",
"ax.set_ylabel(r\"spread of $\\log L$\")\n",
"ax.xaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"{v:.0%}\"))\n",
"ax.yaxis.set_major_formatter(FuncFormatter(lambda v, _: f\"{v:g}\"))\n",
"ax.set_xlim(u_grid[0], u_grid[-1])\n",
"ax.legend(loc=\"lower left\", fontsize=8)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "9106cd50",
"metadata": {},
"source": [
"The arrows are the cutoffs where the profile interval is unbounded, so its \"standard deviation\" is infinite.\n",
"\n",
"## All the numbers\n",
"\n",
"For reference, the intervals behind the interval figure, as MLE with 95% profile-likelihood interval, and\n",
"posterior median with 95% equal-tailed credible interval. The truth is $L = 1$ and $t_0 = 0$. The last column\n",
"is the posterior probability that $L$ is above the top of the profile interval, which is the plateau effect in\n",
"one number."
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "b8327b1d",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"
\n",
" \n",
"
\n",
"
y(T)/L
\n",
"
L, profile
\n",
"
L, posterior
\n",
"
t₀ − T, profile
\n",
"
t₀ − T, posterior
\n",
"
P(L > profile top)
\n",
"
\n",
" \n",
" \n",
"
\n",
"
1%
\n",
"
∞ [0.0563, ∞]
\n",
"
11.4 [0.103, 865]
\n",
"
∞ [1.61, ∞]
\n",
"
7.09 [2.27, 11.5]
\n",
"
—
\n",
"
\n",
"
\n",
"
2%
\n",
"
0.577 [0.0958, ∞]
\n",
"
11 [0.153, 866]
\n",
"
3.34 [1.44, ∞]
\n",
"
6.36 [1.96, 10.8]
\n",
"
—
\n",
"
\n",
"
\n",
"
5%
\n",
"
1.42 [0.252, ∞]
\n",
"
17.5 [0.401, 879]
\n",
"
3.34 [1.49, ∞]
\n",
"
5.88 [1.99, 9.83]
\n",
"
—
\n",
"
\n",
"
\n",
"
10%
\n",
"
0.577 [0.314, 3.33]
\n",
"
0.838 [0.355, 423]
\n",
"
1.65 [0.958, 3.5]
\n",
"
2.05 [1.1, 8.46]
\n",
"
20%
\n",
"
\n",
"
\n",
"
20%
\n",
"
1.22 [0.681, 5.23]
\n",
"
1.54 [0.749, 202]
\n",
"
1.64 [0.967, 3.17]
\n",
"
1.87 [1.08, 6.93]
\n",
"
14%
\n",
"
\n",
"
\n",
"
50%
\n",
"
1 [0.834, 1.23]
\n",
"
1 [0.841, 1.23]
\n",
"
0 [-0.216, 0.273]
\n",
"
0.0177 [-0.204, 0.277]
\n",
"
3%
\n",
"
\n",
"
\n",
"
80%
\n",
"
1.11 [0.987, 1.19]
\n",
"
1.08 [0.99, 1.19]
\n",
"
-1.29 [-1.41, -1.14]
\n",
"
-1.28 [-1.41, -1.14]
\n",
"
2%
\n",
"
\n",
" \n",
"
\n"
],
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def ci(x):\n",
" lo, mid, hi = np.percentile(x, [2.5, 50, 97.5])\n",
" return f\"{mid:.3g} [{lo:.3g}, {hi:.3g}]\"\n",
"\n",
"\n",
"def pl(m, lo, hi, f=lambda v: v):\n",
" return f\"{fmt(f(m))} [{fmt(f(lo))}, {fmt(f(hi))}]\"\n",
"\n",
"\n",
"numbers = pd.DataFrame({\n",
" \"y(T)/L\": [f\"{u:.0%}\" for u in REACHED],\n",
" \"L, profile\": [pl(mle[u], *profile_int[u], f=np.exp) for u in REACHED],\n",
" \"L, posterior\": [ci(np.exp(samples[u].log_l)) for u in REACHED],\n",
" \"t₀ − T, profile\": [pl(mle_t0[u] - T, profile_int_t0[u][0] - T, profile_int_t0[u][1] - T)\n",
" for u, T in zip(REACHED, T_CUT)],\n",
" \"t₀ − T, posterior\": [ci(samples[u].t0 - T) for u, T in zip(REACHED, T_CUT)],\n",
" \"P(L > profile top)\": [\n",
" \"—\" if np.isinf(profile_int[u][1]) else f\"{(samples[u].log_l > profile_int[u][1]).mean():.0%}\"\n",
" for u in REACHED],\n",
"}).style.hide(axis=\"index\").set_uuid(\"numbers\")\n",
"HTML(numbers.to_html())"
]
},
{
"cell_type": "markdown",
"id": "72df6140",
"metadata": {},
"source": [
"## Caveats\n",
"\n",
"- **One realization.** Everything here is one draw of the noise, cut at seven places. Another draw would move\n",
" the MLEs and the interval ends around, especially at 1–5%, though not the overall picture.\n",
"- **Coverage.** The profile intervals rely on Wilks' theorem, which is asymptotic, and this is the regime where\n",
" it is least trustworthy. Checking their coverage would need many simulated datasets with a fit to each.\n",
"- **$\\sigma$ known.** Fitting the noise level as well would widen everything a little and change nothing\n",
" qualitatively.\n",
"- **The prior.** Every statement about the posterior's upper tail depends on the bound at $e^{7}$. The profile\n",
" likelihood does not."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.15"
}
},
"nbformat": 4,
"nbformat_minor": 5
}