{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# KPT audit V8 - numerical evidence for the V8 protocol\n", "\n", "This notebook is the evidence base for `TCIU_SK_Appendix03_KPT_DS_fMRI_V8.Rmd`. The two files are a\n", "pair: every row of the changelog table in the Rmd cites a section number here, and every finding here\n", "names the R function that implements the fix. Nothing in the Rmd is asserted without a corresponding\n", "experiment below, run in NumPy and therefore independent of the R implementation.\n", "\n", "The audit targets `TCIU_SK_Appendix03_KPT_DS_fMRI_V6.Rmd` and the manuscript\n", "`IOP_IP_SKA_KPT_V6_Phy.tex` (model, Lemma \"Convolution structure\", Theorem \"Identifiability of phase\n", "distributions\", Theorem \"Asymptotic null for de-biased or undersmoothed Rayleigh inference\").\n", "\n", "**Findings**\n", "\n", "| # | Finding | Section | Implemented in the V8 Rmd by |\n", "|---|---|---|---|\n", "| 1 | The V6 anchor rotates with the wrong sign: a density with mode $\\mu$ lands at $2\\mu$, not $0$ | \u00a72 | `anchor_rows()`, self-test 3 |\n", "| 2 | Both V6 simulators use surfaces even in $\\theta$, so $\\varphi$ is identifiable only up to reflection and the first-harmonic information matrix is exactly singular | \u00a73 | `reflection_diagnostic()`, `break_reflection=` |\n", "| 3 | Identifiability condition (I2), $\\inf_\\theta\\|\\mathcal S\\|>0$, fails in the V6 fMRI simulator | \u00a73 | `spectral_diagnostic()`, `well_posed=` |\n", "| 4 | Applying the $\\chi^2_2$ Rayleigh rule to a degenerate surface rejects at rate 1.000 under the null; a rank-aware statistic restores calibration | \u00a74 | `rayleigh_uniformity()`, `validate_rayleigh()` |\n", "| 5 | The unguarded alternating KPT iteration diverges; a safeguarded damped step is monotone and stable | \u00a75 | `kpt_fit(safeguard=TRUE)`, self-test 7 |\n", "| 6 | KPT beats the uniform-phase null, but the V6 anchor sign error hides most of the advantage, and the gain is concentrated where the truth is non-uniform | \u00a76 | `rep_uniform()` and the skill section |\n", "| 7 | The V6 regularizer is a pointwise weight in $\\theta$, not a Sobolev seminorm, but the difference is immaterial at the paper's $\\lambda$ | \u00a77 | documented in `kpt_fit()` |\n", "| 8 | Averaging two partial sums of the alternating ILT tail gives $6\\times10^{-3}$ error on $e^{-t}$; Euler weights give $6\\times10^{-9}$ | \u00a78 | `ilt_valsa()`, self-test 9 |\n", "\n", "Run time is about two minutes with `FAST = True`. The companion R notebook is\n", "`TCIU_SK_Appendix03_KPT_DS_fMRI_V8.Rmd`; its self-tests re-check findings 1, 2, 5 and 8 at render time." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from scipy import stats\n", "import matplotlib.pyplot as plt\n", "\n", "FAST = True # False reproduces the larger Monte Carlo runs quoted in the text\n", "B_SIZE = 600 if FAST else 2000 # replications for the null-calibration study\n", "B_POWER = 300 if FAST else 1000\n", "\n", "rng_global = np.random.default_rng(20260912)\n", "np.set_printoptions(precision=4, suppress=True)\n", "print(\"numpy\", np.__version__)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Shared helpers\n", "\n", "The model is the time-specific latent-phase model of the manuscript,\n", "\n", "$$Y_{j,k}=\\mathcal S(t_k,\\Theta_j(t_k))+\\varepsilon_{j,k},\\qquad \\Theta_j(t_k)\\sim\\varphi_{t_k},\n", "\\qquad \\varepsilon_{j,k}\\sim\\mathcal N(0,\\sigma^2),$$\n", "\n", "with all angular integrals taken against the Haar probability measure $d\\theta/(2\\pi)$, so every phase\n", "density has **mean one** on the grid and the uniform law is the vector of ones.\n", "\n", "The helpers below mirror the R routines exactly, including R's FFT sign conventions, so that a\n", "discrepancy found here is a discrepancy in the R code and not in a translation." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "L = 256\n", "theta = 2*np.pi*np.arange(L)/L\n", "\n", "def von_mises(th, mu, kappa):\n", " '''von Mises density normalised to mean one (Haar convention).'''\n", " d = np.exp(kappa*np.cos(th-mu))/(2*np.pi*np.i0(kappa))\n", " return d/d.mean()\n", "\n", "def freq_index(n):\n", " '''R's fft frequency ordering.'''\n", " if n % 2 == 0:\n", " return np.concatenate([np.arange(0, n//2+1), np.arange(-(n//2-1), 0)])\n", " return np.concatenate([np.arange(0, (n-1)//2+1), np.arange(-((n-1)//2), 0)])\n", "\n", "def shift_periodic(x, delta):\n", " '''R's shift_periodic_row: returns x(theta - delta), computed spectrally.'''\n", " n = freq_index(len(x))\n", " return np.real(np.fft.ifft(np.fft.fft(x)*np.exp(-1j*n*delta)))\n", "\n", "def simplex_project(v, z=None):\n", " '''Exact Euclidean projection onto {x >= 0, mean(x) = 1}.'''\n", " n = len(v); z = n if z is None else z\n", " u = np.sort(v)[::-1]\n", " css = np.cumsum(u) - z\n", " rho = np.max(np.where(u - css/np.arange(1, n+1) > 0)[0])\n", " return np.maximum(v - css[rho]/(rho+1), 0)\n", "\n", "def phase_coeff(row, th, n):\n", " '''phi-hat(n) = mean(phi * exp(-i n theta)).'''\n", " return np.mean(row*np.exp(-1j*n*th))\n", "\n", "def jsd_rows(P, Q):\n", " '''Jensen-Shannon divergence in bits, row by row.'''\n", " P = P/P.sum(1, keepdims=True); Q = Q/Q.sum(1, keepdims=True)\n", " M = 0.5*(P+Q)\n", " f = lambda A, B: np.sum(A*np.log2(np.maximum(A, 1e-15)/np.maximum(B, 1e-15)), 1)\n", " return 0.5*f(P, M) + 0.5*f(Q, M)\n", "\n", "# harmonic truncation and the coefficient <-> omega maps used by the KPT filter\n", "J = 12\n", "nv = np.arange(-J, J+1); Nn = len(nv)\n", "E_neg = np.exp(-1j*np.outer(theta, nv)) # L x Nn\n", "E_pos = np.exp( 1j*np.outer(theta, nv))\n", "omega = 2*np.pi*np.arange(Nn)/Nn\n", "Lam = (2*np.sin(omega/2))**2 # the V6 regularisation symbol, p = 1\n", "\n", "def to_omega(c):\n", " z = np.zeros(Nn, complex)\n", " z[(nv + Nn) % Nn] = c\n", " return np.fft.fft(z)\n", "\n", "def from_omega(Z):\n", " y = np.fft.ifft(Z)\n", " return y[(nv + Nn) % Nn]\n", "\n", "def coeffs(x): return (x @ E_neg)/L\n", "def synth(c): return np.real(E_pos @ c)\n", "\n", "print(\"grid L =\", L, \" harmonics |n| <=\", J)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Finding 1 \u2014 the anchoring rotation has the wrong sign\n", "\n", "The rotational gauge is a genuine model symmetry: $(\\mathcal S(\\cdot-\\alpha),\\varphi(\\cdot-\\alpha))$ gives\n", "the same law of $Y$. The manuscript fixes it by choosing $\\alpha$ so that\n", "$e^{-in_\\star\\alpha}\\widehat\\varphi(n_\\star)\\in(0,\\infty)$, i.e. $\\alpha=\\arg\\widehat\\varphi(n_\\star)/n_\\star$.\n", "\n", "V6 instead computes `m = mean(phi * exp(+1i*theta))` and rotates by `+Arg(m)`. Since\n", "$m=\\widehat\\varphi(-1)=\\overline{\\widehat\\varphi(1)}$, that is a rotation by $-\\arg\\widehat\\varphi(1)$ \u2014 the\n", "wrong sign. The consequence is that an anchored density has its mode at $2\\mu$ rather than at $0$, so two\n", "rows are mapped to a **common frame only if their first moments already agree**." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def anchor_v6(row, th):\n", " '''V6 rule: delta = +Arg(mean(phi e^{+i theta})).'''\n", " return shift_periodic(row, np.angle(np.mean(row*np.exp(1j*th))))\n", "\n", "def anchor_fixed(row, th, n_star=1):\n", " '''Manuscript rule: alpha = Arg(phi-hat(n_star)) / n_star.'''\n", " return shift_periodic(row, np.angle(phase_coeff(row, th, n_star))/n_star)\n", "\n", "mode_of = lambda g: theta[np.argmax(g)]\n", "print(f\"{'true mode mu':>13s} {'V6 anchored':>13s} {'corrected':>11s}\")\n", "for mu in (0.4, 1.0, 2.0, 4.5):\n", " f = von_mises(theta, mu, 4.0)\n", " print(f\"{mu:13.3f} {mode_of(anchor_v6(f, theta)):13.3f} {mode_of(anchor_fixed(f, theta)):11.3f}\")\n", "print(\"\\nThe corrected rule sends every mode to 0; the V6 rule sends mode mu to 2*mu (mod 2*pi).\")\n", "\n", "# the practical consequence: two rows land in a common frame only under the corrected rule\n", "a, b = von_mises(theta, 0.4, 4.0), von_mises(theta, 2.4, 4.0)\n", "d_v6 = 0.5*np.abs(anchor_v6(a, theta)/L - anchor_v6(b, theta)/L).sum()\n", "d_fix = 0.5*np.abs(anchor_fixed(a, theta)/L - anchor_fixed(b, theta)/L).sum()\n", "print(f\"\\nTV distance between two anchored von Mises rows that differ only by rotation:\")\n", "print(f\" V6 rule : {d_v6:.4f} (should be 0 - the rows are the same shape)\")\n", "print(f\" corrected rule : {d_fix:.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Findings 2 and 3 \u2014 the V6 surfaces are not identifiable\n", "\n", "Because $F(e^{i\\theta},t)=\\mathcal S(t,\\theta)$ and $\\Phi(e^{i\\theta},t)=\\varphi_t(\\theta)$, the\n", "factorization $M=F\\Phi$ of the manuscript's convolution lemma is *pointwise multiplication in $\\theta$*.\n", "Every identifiability property of the inverse problem is therefore a property of the surface:\n", "\n", "* **Reflection symmetry.** If $\\mathcal S(t,\\alpha+u)=\\mathcal S(t,\\alpha-u)$ then $\\varphi$ and its\n", " reflection about $\\alpha$ induce the same law of $Y$. Anchoring does not remove this.\n", "* **Spectral lower bound (I2).** The identifiability theorem needs $\\inf_\\theta|\\mathcal S|\\ge c_F>0$.\n", " Where the surface crosses zero, the data carry no information about $\\varphi$ and the regularizer\n", " decides the answer.\n", "\n", "Both V6 simulators use $\\mathcal S=\\text{baseline}+a\\cos\\theta+b\\cos2\\theta$, which is even **and**\n", "crosses zero." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def reflection_asymmetry(S, th, n_axis=180):\n", " '''Smallest normalised asymmetry over all reflection axes; ~0 means reflection-degenerate.'''\n", " axes = np.linspace(0, np.pi, n_axis)\n", " v = S.var()\n", " r = [np.mean((np.interp((a+th) % (2*np.pi), th, S, period=2*np.pi) -\n", " np.interp((a-th) % (2*np.pi), th, S, period=2*np.pi))**2)/max(4*v, 1e-15)\n", " for a in axes]\n", " return float(np.min(r)), float(axes[int(np.argmin(r))])\n", "\n", "S_v6 = 0.05 + 1.0*np.cos(theta) + 0.25*np.cos(2*theta) # V6 fMRI surface\n", "S_v8 = 1.60 + 1.0*np.cos(theta) + 0.25*np.cos(2*theta) + 0.35*np.sin(2*theta) # V8 default\n", "\n", "for name, S in ((\"V6 (even, zero-crossing)\", S_v6), (\"V8 (generic, DC offset)\", S_v8)):\n", " asym, axis = reflection_asymmetry(S, theta)\n", " print(f\"{name:26s} min|S| = {np.abs(S).min():.4f} reflection asymmetry = {asym:.2e} \"\n", " f\"(axis {axis:.2f}) (I2) {'FAILS' if np.abs(S).min() < 0.05*np.abs(S).max() else 'holds'}\")\n", "\n", "fig, ax = plt.subplots(1, 2, figsize=(10, 3.2))\n", "for a, S, ttl in ((ax[0], S_v6, \"V6: even and touches zero\"), (ax[1], S_v8, \"V8: generic, bounded away from 0\")):\n", " a.plot(theta, S, lw=1.6); a.axhline(0, color=\"grey\", lw=0.7)\n", " a.set_xticks([0, np.pi, 2*np.pi]); a.set_xticklabels([\"0\", r\"$\\pi$\", r\"$2\\pi$\"])\n", " a.set_title(ttl, fontsize=10); a.set_xlabel(r\"$\\theta$\"); a.set_ylabel(r\"$S(t_k,\\theta)$\")\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Finding 4 \u2014 the uniformity test, and how to calibrate it\n", "\n", "The manuscript's Rayleigh-type statistic de-biases the regularized first harmonic and forms\n", "$R_N=N\\,\\tilde u^\\top\\widehat\\Sigma^{-1}\\tilde u$ with $\\tilde u=(\\Re\\tilde\\varphi(1),\\Im\\tilde\\varphi(1))$,\n", "with a $\\chi^2_2$ limit under $H_0:\\varphi_t\\equiv1$.\n", "\n", "The per-replicate contribution to $\\tilde\\varphi(1)$ is a linear functional of\n", "$\\xi_{j,n}=Y_j\\,\\mathbb E[e^{-in\\Theta}\\mid Y_j]$, so $\\widehat\\Sigma$ can be estimated by the empirical\n", "covariance of those contributions. **For a reflection-symmetric surface every contribution is real**, so\n", "$\\widehat\\Sigma$ is singular and the $\\chi^2_2$ reference is invalid.\n", "\n", "The fix is to threshold the eigenvalues of $\\widehat\\Sigma$ and refer the statistic to $\\chi^2_r$ with\n", "$r$ the effective rank." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def rayleigh(y, S_row, sigma, tau=1e-4, harmonic=1, rank_tol=1e-8, force_df=None):\n", " '''De-biased first-harmonic statistic with a rank-aware quadratic form.'''\n", " N = len(y)\n", " ll = -(y[:, None] - S_row[None, :])**2/(2*sigma**2) # uniform prior under H0\n", " ll -= ll.max(1, keepdims=True)\n", " W = np.exp(ll); W /= W.sum(1, keepdims=True)\n", " xi = y[:, None]*(W @ E_neg) # N x Nn\n", " F = to_omega(coeffs(S_row))\n", " Hinv = np.conj(F)/(np.abs(F)**2 + tau) # de-biasing filter\n", " Xi = np.zeros((Nn, N), complex); Xi[(nv + Nn) % Nn, :] = xi.T\n", " Cc = np.fft.ifft(np.fft.fft(Xi, axis=0)*Hinv[:, None], axis=0)\n", " c = Cc[(harmonic + Nn) % Nn, :]\n", " u = np.column_stack([c.real, c.imag])\n", " ubar = u.mean(0); Sig = np.cov(u, rowvar=False)\n", " w, V = np.linalg.eigh(Sig)\n", " keep = w > rank_tol*max(w.max(), 1e-300)\n", " if force_df is not None: # the unguarded chi2_2 rule\n", " z = V.T @ ubar\n", " return float(N*np.sum(z**2/np.maximum(w, 1e-300))), force_df\n", " z = V[:, keep].T @ ubar\n", " return float(N*np.sum(z**2/w[keep])), int(keep.sum())\n", "\n", "sigma, N = 0.2, 200\n", "print(\"Effective rank of the first-harmonic information matrix:\")\n", "for name, S in ((\"V6 (even)\", S_v6), (\"V8 (generic)\", S_v8)):\n", " th_j = rng_global.uniform(0, 2*np.pi, N)\n", " y = np.interp(th_j, theta, S, period=2*np.pi) + rng_global.normal(0, sigma, N)\n", " _, r = rayleigh(y, S, sigma)\n", " print(f\" {name:14s} effective rank = {r}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def null_sizes(S, sigma=0.2, N=200, B=B_SIZE, force_df=None, seed=1):\n", " rng = np.random.default_rng(seed)\n", " st, dfs = [], []\n", " for _ in range(B):\n", " idx = rng.integers(0, L, N)\n", " y = S[idx] + rng.normal(0, sigma, N)\n", " s, r = rayleigh(y, S, sigma, force_df=force_df)\n", " st.append(s); dfs.append(r)\n", " st = np.array(st); r = int(np.median(dfs))\n", " rows = []\n", " for a in (0.10, 0.05, 0.01):\n", " rows.append((a, float((st > stats.chi2.ppf(1-a, r)).mean()), float(np.sqrt(a*(1-a)/B))))\n", " ks = stats.kstest(st, lambda x: stats.chi2.cdf(x, r)).pvalue\n", " return rows, r, ks\n", "\n", "print(f\"Empirical type I error, B = {B_SIZE} replications\\n\")\n", "print(f\"{'surface':26s} {'rule':>12s} {'df':>3s} a=0.10 a=0.05 a=0.01 KS p\")\n", "for name, S, force in ((\"V6 (even)\", S_v6, 2),\n", " (\"V6 (even)\", S_v6, None),\n", " (\"V8 (generic)\", S_v8, None)):\n", " rows, r, ks = null_sizes(S, force_df=force)\n", " rule = \"chi2_2 (V6)\" if force else \"rank-aware\"\n", " print(f\"{name:26s} {rule:>12s} {r:3d} \" + \" \".join(f\"{v:7.3f}\" for _, v, _ in rows) + f\" {ks:5.3f}\")\n", "print(f\"\\nMonte Carlo standard errors: {np.sqrt(0.10*0.90/B_SIZE):.3f}, \"\n", " f\"{np.sqrt(0.05*0.95/B_SIZE):.3f}, {np.sqrt(0.01*0.99/B_SIZE):.3f}\")\n", "print(\"The unguarded chi2_2 rule on a degenerate surface rejects essentially always under the null.\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Power of the rank-aware test on the generic surface\n", "kappas = [0.0, 0.25, 0.5, 1.0, 2.0]\n", "rng = np.random.default_rng(4)\n", "power = []\n", "for kap in kappas:\n", " rej = []\n", " for _ in range(B_POWER):\n", " if kap == 0:\n", " idx = rng.integers(0, L, N)\n", " else:\n", " p = von_mises(theta, 1.0, kap); p = p/p.sum()\n", " idx = rng.choice(L, size=N, p=p)\n", " y = S_v8[idx] + rng.normal(0, sigma, N)\n", " s, r = rayleigh(y, S_v8, sigma)\n", " rej.append(s > stats.chi2.ppf(0.95, r))\n", " power.append(np.mean(rej))\n", " print(f\" kappa = {kap:4.2f} rejection rate = {power[-1]:.3f}\")\n", "\n", "plt.figure(figsize=(5, 3))\n", "plt.axhline(0.05, ls=\"--\", c=\"grey\", lw=0.8)\n", "plt.plot(kappas, power, \"o-\")\n", "plt.ylim(0, 1.05); plt.xlabel(r\"$\\kappa$ (phase concentration)\"); plt.ylabel(\"rejection rate\")\n", "plt.title(f\"Power at the 5% level (N = {N})\", fontsize=10)\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Finding 5 \u2014 the unguarded iteration diverges\n", "\n", "The manuscript is explicit that the projected practical routines are not exact EM and carries no\n", "monotonicity claim for them. V6 nevertheless runs a fixed number of unguarded updates. Below, the plain\n", "alternating scheme is compared with a safeguarded version that damps the proposed step and accepts it only\n", "if the observed-data log-likelihood does not decrease." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def simulate(K=40, N=150, sigma=0.2, seed=3):\n", " '''Generative model with a non-even surface and a time-varying phase law.'''\n", " rng = np.random.default_rng(seed)\n", " t = np.linspace(0, 1, K)\n", " bold = 0.3 + 0.7*(np.sin(2*np.pi*t*3) > 0)\n", " S = np.array([0.9 + b*np.cos(theta) + 0.35*b*np.sin(2*theta) + 0.2*np.cos(3*theta) for b in bold])\n", " mu, kap = 1.2*np.sin(2*np.pi*t), np.where(bold > 0.5, 5.0, 1.0)\n", " P = np.array([von_mises(theta, mu[k], kap[k]) for k in range(K)])\n", " def draw(n):\n", " Y = np.zeros((n, K))\n", " for k in range(K):\n", " p = P[k]/P[k].sum()\n", " Y[:, k] = S[k][rng.choice(L, size=n, p=p)] + rng.normal(0, sigma, n)\n", " return Y\n", " return dict(Y=draw(N), Y_test=draw(N), S=S, P=P, t=t, bold=bold, sigma=sigma)\n", "\n", "def obs_loglik(Y, S, P, sigma):\n", " tot = 0.0\n", " for k in range(Y.shape[1]):\n", " w = P[k]/P[k].sum()\n", " d = np.exp(-(Y[:, k][:, None] - S[k][None, :])**2/(2*sigma**2))\n", " tot += np.sum(np.log(np.maximum(d @ w, 1e-300)))\n", " return tot\n", "\n", "def kpt(Y, S0, sigma, lam=1e-2, iters=20, damping=1.0, safeguard=False):\n", " K = Y.shape[1]\n", " P = np.ones((K, L)); S = S0.copy()\n", " ll = obs_loglik(Y, S, P, sigma); trace = [ll]\n", " for _ in range(iters):\n", " Pn, Sn = P.copy(), S.copy()\n", " for k in range(K):\n", " lo = -(Y[:, k][:, None] - S[k][None, :])**2/(2*sigma**2)\n", " lo = lo + np.log(np.maximum(P[k]/P[k].sum(), 1e-300))\n", " lo -= lo.max(1, keepdims=True)\n", " W = np.exp(lo); W /= W.sum(1, keepdims=True)\n", " m = (Y[:, k][:, None]*(W @ E_neg)).mean(0)\n", " F, M = to_omega(coeffs(S[k])), to_omega(m)\n", " Pn[k] = simplex_project(synth(from_omega(np.conj(F)*M/(np.abs(F)**2 + lam*Lam + 1e-10))))\n", " Pw = to_omega(coeffs(Pn[k]))\n", " Sn[k] = synth(from_omega(np.conj(Pw)*M/(np.abs(Pw)**2 + lam*Lam + 1e-10)))\n", " step, ok = 1.0, False\n", " for _ in range(6):\n", " w = step*damping\n", " S_try = (1-w)*S + w*Sn\n", " P_try = np.array([simplex_project(r) for r in (1-w)*P + w*Pn])\n", " ll_try = obs_loglik(Y, S_try, P_try, sigma)\n", " if (not safeguard) or ll_try >= ll - 1e-8:\n", " S, P, ll, ok = S_try, P_try, ll_try, True\n", " break\n", " step /= 2\n", " trace.append(ll)\n", " if not ok:\n", " break\n", " return P, S, np.array(trace)\n", "\n", "def anchor_mat(P, S):\n", " d = np.array([np.angle(phase_coeff(r, theta, 1)) for r in P])\n", " return (np.array([shift_periodic(P[k], d[k]) for k in range(len(d))]),\n", " np.array([shift_periodic(S[k], d[k]) for k in range(len(d))]))\n", "\n", "def rel_L2_cone(A, B, t):\n", " e, r = np.mean((A-B)**2, 1), np.mean(A**2, 1)\n", " tz = lambda y: np.trapezoid(t*y, t)\n", " return float(np.sqrt(tz(e)/max(tz(r), 1e-15)))\n", "\n", "sim = simulate()\n", "S_init = np.array([np.real(np.fft.ifft(np.fft.fft(\n", " np.interp(theta, 2*np.pi*np.arange(sim[\"Y\"].shape[0])/sim[\"Y\"].shape[0],\n", " np.sort(sim[\"Y\"][:, k]), period=2*np.pi))*(np.abs(np.fft.fftfreq(L, 1/L)) <= J)))\n", " for k in range(sim[\"Y\"].shape[1])])\n", "Pt_a, St_a = anchor_mat(sim[\"P\"], sim[\"S\"])\n", "\n", "print(f\"{'scheme':24s} {'iters':>5s} {'phase JSD':>10s} {'surface relL2':>14s} {'final loglik':>13s} monotone\")\n", "traces = {}\n", "for name, sg, dp, it in ((\"plain alternating\", False, 1.0, 20), (\"plain alternating\", False, 1.0, 60),\n", " (\"safeguarded + damped\", True, 0.5, 20), (\"safeguarded + damped\", True, 0.5, 60)):\n", " P, S, tr = kpt(sim[\"Y\"], S_init.copy(), sim[\"sigma\"], iters=it, damping=dp, safeguard=sg)\n", " Pa, Sa = anchor_mat(P, S)\n", " mono = bool(np.all(np.diff(tr) >= -1e-6))\n", " traces[f\"{name} ({it})\"] = tr\n", " print(f\"{name:24s} {it:5d} {jsd_rows(Pt_a, Pa).mean():10.4f} {rel_L2_cone(St_a, Sa, sim['t']):14.4f} \"\n", " f\"{tr[-1]:13.1f} {mono}\")\n", "\n", "plt.figure(figsize=(6, 3.2))\n", "for nm, tr in traces.items():\n", " if \"60\" in nm:\n", " plt.plot(tr, label=nm)\n", "plt.xlabel(\"iteration\"); plt.ylabel(\"observed-data log-likelihood\")\n", "plt.legend(fontsize=8); plt.title(\"Unguarded updates can destroy the fit\", fontsize=10)\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Finding 6 \u2014 KPT versus the uniform-phase null, stratified\n", "\n", "Two separate effects make a naive average misleading here.\n", "\n", "**The anchor convention changes the answer.** Measuring the *same fit* with the V6 anchor and with the\n", "corrected anchor gives materially different divergences, because the buggy sign leaves truth and estimate\n", "in different frames. This is shown first below: it accounts for most of the apparent gap between KPT and\n", "the null.\n", "\n", "**The average dilutes a real effect.** Half the time points are rest blocks whose true phase law is nearly\n", "uniform. There the null is nearly correct by construction and no estimator can gain, so including those\n", "times in an average understates the method. Stratifying by the mean resultant length\n", "$R(t)=|\\widehat\\varphi_t(1)|$ \u2014 zero for the uniform law, approaching one as the law concentrates \u2014\n", "separates the regimes. The skill score $1-\\mathrm{JSD}_{\\text{method}}/\\mathrm{JSD}_{\\text{null}}$ is\n", "positive when the method beats the null." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "P_kpt, S_kpt, _ = kpt(sim[\"Y\"], S_init.copy(), sim[\"sigma\"], iters=25, damping=0.5, safeguard=True)\n", "Pk_a, _ = anchor_mat(P_kpt, S_kpt)\n", "K = sim[\"Y\"].shape[1]\n", "\n", "# --- how much of the KPT-vs-null comparison is decided by the anchoring sign alone? ---\n", "def anchor_with(P, buggy):\n", " d = np.array([np.angle(np.mean(r*np.exp(1j*theta))) if buggy\n", " else np.angle(phase_coeff(r, theta, 1)) for r in P])\n", " return np.array([shift_periodic(P[k], d[k]) for k in range(len(d))])\n", "\n", "print(\"Same fit, same data, two anchoring conventions:\")\n", "for nm, bug in ((\"V6 (buggy sign)\", True), (\"corrected sign\", False)):\n", " jk = jsd_rows(anchor_with(sim[\"P\"], bug), anchor_with(P_kpt, bug)).mean()\n", " jn = jsd_rows(anchor_with(sim[\"P\"], bug), np.ones((K, L))).mean()\n", " print(f\" {nm:18s} mean JSD: KPT {jk:.4f} null {jn:.4f} skill {1-jk/jn:+.3f}\")\n", "print(\" The anchor bug alone removes most of the measured advantage.\\n\")\n", "\n", "j_kpt = jsd_rows(Pt_a, Pk_a)\n", "j_null = jsd_rows(Pt_a, np.ones((K, L)))\n", "skill = 1 - j_kpt/np.maximum(j_null, 1e-12) # > 0 means better than the null\n", "R_true = np.abs([phase_coeff(sim[\"P\"][k], theta, 1) for k in range(K)])\n", "on = sim[\"bold\"] > 0.5\n", "\n", "print(f\"overall: mean JSD KPT {j_kpt.mean():.4f} null {j_null.mean():.4f} mean skill {skill.mean():+.3f}\")\n", "print(f\"paired Wilcoxon over all {K} times: p = {stats.wilcoxon(j_null - j_kpt).pvalue:.4f}\\n\")\n", "qs = np.quantile(R_true, [0, 1/3, 2/3, 1.0])\n", "print(f\"{'R(t) tercile':>22s} {'n':>3s} {'JSD KPT':>9s} {'JSD null':>9s} {'skill':>7s} {'win rate':>9s}\")\n", "for i in range(3):\n", " m = (R_true >= qs[i]) & (R_true <= qs[i+1])\n", " print(f\"[{qs[i]:.2f}, {qs[i+1]:.2f}]{'':>10s} {m.sum():3d} {j_kpt[m].mean():9.4f} {j_null[m].mean():9.4f} \"\n", " f\"{skill[m].mean():+7.3f} {np.mean(j_kpt[m] < j_null[m]):9.2f}\")\n", "print()\n", "for nm, m in ((\"ON blocks (kappa=5)\", on), (\"OFF blocks (kappa=1)\", ~on)):\n", " p = stats.wilcoxon(j_null[m] - j_kpt[m]).pvalue\n", " print(f\"{nm}: skill {skill[m].mean():+.3f} win rate {np.mean(j_kpt[m] < j_null[m]):.2f} paired p = {p:.4f}\")\n", "\n", "plt.figure(figsize=(5.4, 3.2))\n", "plt.axhline(0, color=\"grey\", lw=0.8)\n", "plt.scatter(R_true[on], skill[on], s=18, label=\"ON blocks\")\n", "plt.scatter(R_true[~on], skill[~on], s=18, marker=\"s\", label=\"OFF blocks\")\n", "plt.xlabel(r\"$R(t)=|\\hat\\varphi_t(1)|$ (0 = uniform)\"); plt.ylabel(\"skill vs uniform null\")\n", "plt.legend(fontsize=8); plt.title(\"KPT gains only where the truth is concentrated\", fontsize=10)\n", "plt.tight_layout(); plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Finding 7 \u2014 the regularizer is misnamed but harmless here\n", "\n", "V6 calls its penalty Wiener\u2013Sobolev and uses the symbol $\\Lambda(\\omega)=(2\\sin(\\omega/2))^{2p}$. Because\n", "the transformed variable $\\omega$ is dual to the harmonic index $n$, and therefore *is* $\\theta$, this acts\n", "as a pointwise weight $\\int(2\\sin(\\theta/2))^{2p}|\\varphi|^2$ rather than as a Sobolev seminorm\n", "$\\sum_n|n|^{2p}|\\widehat\\varphi(n)|^2$. The weight is not rotation invariant, so it could in principle drag\n", "mass toward $\\theta=0$.\n", "\n", "We tested for that bias and **did not find it** at the manuscript's $\\lambda$: reconstructing a von Mises\n", "law placed at four different rotations gives the correct mode and total-variation errors below $0.01$ in\n", "every case. The naming is imprecise; the numerical behaviour is fine. V8 documents the penalty precisely\n", "and keeps it." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "S_pos = 1.6 + np.cos(theta) + 0.25*np.cos(2*theta) # strictly positive, so (I2) holds\n", "print(f\"{'mu':>6s} {'recovered mode':>15s} {'TV error':>10s}\")\n", "for mu in (0.0, np.pi/2, np.pi, 3*np.pi/2):\n", " phi = von_mises(theta, mu, 4.0)\n", " M = to_omega(coeffs(S_pos*phi)) # noiseless mixed moments: M(theta)=S(theta)phi(theta)\n", " F = to_omega(coeffs(S_pos))\n", " est = synth(from_omega(np.conj(F)*M/(np.abs(F)**2 + 1e-2*Lam + 1e-10)))\n", " est = est/est.mean()\n", " print(f\"{mu:6.2f} {theta[np.argmax(est)]:15.3f} {0.5*np.abs(est-phi).mean():10.4f}\")\n", "print(\"\\nNo rotation-dependent bias at lambda = 1e-2: the mode is recovered exactly in every case.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 8. Auxiliary verifications used by the V8 implementation\n", "\n", "Three identities that the R code relies on." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# (a) exact Euclidean simplex projection\n", "v = rng_global.normal(1, 1, L); p = simplex_project(v)\n", "print(f\"simplex projection: mean = {p.mean():.12f}, min = {p.min():.2e}, \"\n", " f\"idempotent = {np.allclose(p, simplex_project(p))}\")\n", "\n", "# (b) the O(B log B) identity for the CRPS self-term, E|X - X'|\n", "x = rng_global.normal(size=400); B = len(x); xs = np.sort(x)\n", "ident = 2*np.sum((2*np.arange(1, B+1) - B - 1)*xs)/B**2\n", "brute = np.abs(x[:, None] - x[None, :]).mean()\n", "print(f\"E|X-X'| sorted identity = {ident:.10f}, brute force = {brute:.10f}, match = {np.isclose(ident, brute)}\")\n", "\n", "# (c) Valsa-Brancik inverse Laplace transform: the alternating tail needs a genuine Euler\n", "# transformation. Averaging only the last two partial sums leaves an error near 1e-2 on exp(-t);\n", "# binomially weighting the last P partial sums brings it to 1e-9. The R self-test asserts 1e-6,\n", "# so this distinction decides whether the V8 Rmd renders at all.\n", "from scipy.special import comb\n", "\n", "def ilt_naive(Ffun, t, a=8.0, n_terms=200):\n", " n = np.arange(n_terms)\n", " z = (a + 1j*(n + 0.5)*np.pi)/t\n", " part = np.cumsum(((-1.0)**n)*np.imag(Ffun(z)))\n", " return -np.exp(a)/t*0.5*(part[-1] + part[-2])\n", "\n", "def ilt_euler(Ffun, t, a=8.0, n_main=100, n_euler=29):\n", " n = np.arange(n_main + n_euler + 1)\n", " z = (a + 1j*(n + 0.5)*np.pi)/t\n", " part = np.cumsum(((-1.0)**n)*np.imag(Ffun(z)))\n", " w = comb(n_euler, np.arange(n_euler+1))/2.0**n_euler\n", " return -np.exp(a)/t*np.sum(w*part[n_main:])\n", "\n", "cases = ((1.0, lambda z: 1/(z+1), np.exp(-1.0), \"1/(z+1) -> exp(-t)\"),\n", " (2.0, lambda z: 1/z**2, 2.0, \"1/z^2 -> t\"),\n", " (1.5, lambda z: 1/(z**2+1), np.sin(1.5), \"1/(z^2+1) -> sin t\"))\n", "print(f\"{'transform':22s} {'naive err':>12s} {'Euler err':>12s}\")\n", "for t_eval, F, exact, nm in cases:\n", " print(f\"{nm:22s} {abs(ilt_naive(F,t_eval)-exact):12.2e} {abs(ilt_euler(F,t_eval)-exact):12.2e}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 9. Summary\n", "\n", "| Finding | Status | Effect on V6 results |\n", "|---|---|---|\n", "| Anchor sign | Confirmed | Truth and estimate compared in different frames whenever their first moments disagree |\n", "| Reflection degeneracy | Confirmed, exact | $\\varphi$ identifiable only up to reflection; information matrix rank 1 |\n", "| (I2) violated | Confirmed | Deconvolution locally uninformative where the surface crosses zero |\n", "| $\\chi^2_2$ Rayleigh rule | Confirmed invalid on degenerate surfaces | Null rejection rate 1.000; rank-aware version is calibrated |\n", "| Unguarded iteration | Confirmed | Relative surface error grows with iteration count; safeguarding fixes it |\n", "| KPT vs uniform null | Confirmed, nuanced | KPT does win, but the anchor bug hides most of it and the gain is concentrated where the truth is non-uniform |\n", "| Regularizer naming | Confirmed, immaterial | No measurable bias at the manuscript's $\\lambda$ |\n", "| ILT tail summation | Found while testing V8 | Averaging two partial sums gives 1e-2 error on exp(-t); Euler weights give 1e-9 |\n", "\n", "All eight are addressed in `TCIU_SK_Appendix03_KPT_DS_fMRI_V8.Rmd`. Findings 1-5 and 8 change the\n", "reported numbers, finding 6 changes how they should be summarized, and finding 7 changes only the wording.\n", "Findings 1, 2, 5 and 8 are additionally enforced by self-tests there, so a future edit that reintroduces\n", "any of them stops the render instead of silently producing a wrong table." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.12" } }, "nbformat": 4, "nbformat_minor": 5 }