SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

1 Scope and V8 updates

This appendix implements and validates kime-representation that converts repeatedly sampled longitudinal processes into kimesurfaces, i.e. two-dimensional manifolds parameterized by the kime-magnitude \(t\) and the kime-phase \(\theta\). It expands and supersedes the prior version (V6), which implemented a single family of estimators (KPT-GEM and KPT-FFT) and reported reconstruction metrics against a known ground truth. The revised V8 makes three structural changes.

  1. Multiple representation strategies are compared under one protocol. Epoch stitching, rank/monotone transport, the analytic Laplace route, KPT, and a uniform-phase null model are all wrapped in a common interface and scored on the same held-out data.

  2. Identifiability is treated as a measurable property of each data set, not an assumption. Four diagnostics (gauge anchor, rotational aliasing \(d_t\), reflection symmetry, and the lower spectral bound \(\inf|F|\)) are computed and reported before any estimate is interpreted.

  3. Every inferential claim is calibrated by simulation inside the notebook. Type I error, power, and interval coverage are estimated rather than asserted.

Every correction below was found by re-implementing the V6 routines in NumPy and testing them against the manuscript’s own statements. That audit is the notebook KPT_audit_V8.ipynb, whose section numbers are referenced in the last column of the table. The two files are meant to be read together. The notebook identifies some limitations of V6 and implements approporiate revisions.

V6 implemented a single family of estimators (KPT-GEM and KPT-FFT) and reported reconstruction metrics against a known ground truth. V8 keeps that content but makes three structural changes.

  1. Multiple representation strategies are compared under one protocol. Epoch stitching, rank/monotone transport, the analytic Laplace route, KPT, and a uniform-phase null model are all wrapped in a common interface and scored on the same held-out data.
  2. Identifiability is treated as a measurable property of each data set, not an assumption. Four diagnostics (gauge anchor, rotational aliasing \(d_t\), reflection symmetry, and the lower spectral bound \(\inf|F|\)) are computed and reported before any estimate is interpreted.
  3. Every inferential claim is calibrated by simulation inside the notebook. Type I error, power, and interval coverage are estimated rather than asserted.

The table below lists the specific \(V6 \to V8\) modifications, along with the corresponding motivations. The supporting numerical experiments are reproduced in Section 7 and Section 5 so that a reader can confirm them rather than take them on trust.

# Issue in V6 Consequence Resolution here Audit
1 Anchor applied the rotation with the wrong sign: delta = Arg(m) with m the \(e^{+i\theta}\) moment A density with mode \(\mu\) was rotated to mode \(2\mu\) instead of \(0\), so truth and estimate were compared in different frames whenever their first moments disagreed anchor_rows() implements the paper’s rule \(\alpha=\arg\widehat\varphi(n_\star)/n_\star\) with the \(e^{-in\theta}\) convention; verified by unit test §2
2 Both simulators used surfaces even in \(\theta\) (\(S=a\cos\theta+b\cos2\theta+\text{const}\)) \(\varphi\) is identifiable only up to reflection \(\theta\mapsto-\theta\); the information matrix is exactly singular, so any two-dimensional asymptotic inference on the first harmonic is invalid reflection_diagnostic() detects it; default simulators include a genuine sine component; the Rayleigh test is rank-aware §3
3 Identifiability condition (I2), \(\inf_{|z|=1}|F(z,t)|\ge c_F>0\), fails for the V6 surfaces (\(\min_\theta|S|\approx 0\)) The deconvolution is ill-posed exactly where the surface crosses zero, and the regularizer silently determines the answer there spectral_diagnostic() reports \(\min_\theta|S|\) and the condition ratio; simulators expose a well_posed switch; results are reported in both regimes §3
4 The practical alternating scheme ran a fixed number of unsafeguarded updates The iteration can diverge: in our replication the relative surface error grew from 0.73 to 1.86 as iterations went from 20 to 60 Safeguarded damped updates with backtracking on the observed-data log-likelihood; monotone by construction and stable in iteration count §5
5 No null model There was no way to tell whether KPT beat “assume the phase law is uniform” rep_uniform() is scored alongside every other method in every table §6
6 KPT-GEM never updated the surface, yet its surface error was tabulated next to KPT-FFT’s The GEM surface error measured the initializer, not the algorithm Initializer is reported as its own row; GEM is labeled as a fixed-surface phase estimator —
7 Evaluation used only recovery metrics against a truth that is unavailable in practice No way to compare representations on real data Held-out predictive scoring (log-score, CRPS, PIT, interval coverage) is method-agnostic and is the primary comparison —
8 Metrics reported as point values; bootstrap helper defined but never executed No uncertainty on any reported number Bootstrap confidence intervals on all metrics, paired method comparisons, Benjamini–Hochberg control across time points —
9 \(\sigma\) defaulted to 0.5*mad(Y) An arbitrary scale choice propagates into every posterior weight Profile-likelihood estimate with an explicit sensitivity analysis —
10 The regularizer was called Wiener–Sobolev but the symbol \((2\sin(\omega/2))^{2p}\) acts pointwise in \(\theta\), not as a Sobolev seminorm Naming and interpretation mismatch Documented precisely; we verified it is immaterial at the paper’s \(\lambda\) (total-variation error \(\le 0.006\) across rotations) and kept it, with the alternative available §7
11 The numerical inverse Laplace transform averaged only the last two partial sums of the alternating tail Relative error \(\approx 6\times10^{-3}\) on \(e^{-t}\), three orders of magnitude worse than the method allows Genuine Euler transformation: \(M=100\) exact terms plus a binomially weighted average of the next \(29\) partial sums, giving \(6\times10^{-9}\); self-test 9 enforces it §8
12 Averaged divergence used to compare estimators Two effects cancel: the anchor bug inflates the divergence of any good estimate, and rest blocks with a near-uniform truth dilute a real gain Skill scores against an explicit null, stratified by how far the truth is from uniform, with paired tests (Section 9) §6

One result is worth stating before the tables, because the earlier V6 of this appendix was slightly incorrect. Measured with the V6 anchor, KPT appears no better than simply assuming uniform phases. Measured with the corrected anchor, on the same fit and the same data, which is better and in V8, the skill against the null moves from \(-0.10\) to \(+0.20\). The anchoring sign error was important, as it silently destroyed the measured advantage of every phase estimator. Section 9 reports the comparison properly, stratified by how much phase structure there is to find.

1.1 Parallel backend

Bootstrap replicates are independent, so the two resampling steps in this notebook parallelize cleanly. The backend below forks on Linux and macOS and builds a PSOCK cluster on Windows, falling back to serial execution when parallel is unavailable or only one core can be spared.

Reproducibility does not depend on the core count. All resampling indices are drawn up front in the master process and handed to the workers, and the scoring seed is a deterministic function of the replicate number. The same seed therefore gives the same answer on 1 core or 32.

par_backend <- function(n_cores = NULL) {
  if (!requireNamespace("parallel", quietly = TRUE))
    return(list(type = "serial", cores = 1L))
  n <- if (is.null(n_cores)) parallel::detectCores(logical = FALSE) else n_cores
  if (is.na(n) || n < 2L) return(list(type = "serial", cores = 1L))
  n <- max(1L, min(as.integer(n) - 1L, 8L))          # leave a core free; cap at 8
  if (n < 2L) return(list(type = "serial", cores = 1L))
  list(type = if (.Platform$OS.type == "unix") "fork" else "psock", cores = n)
}

## Names of every function in the global environment, for PSOCK export.
par_export_names <- function() {
  g <- ls(envir = globalenv())
  g[vapply(g, function(nm) is.function(get(nm, envir = globalenv())), logical(1))]
}

par_lapply <- function(X, FUN, backend = PAR, export_extra = character()) {
  if (backend$type == "serial") return(lapply(X, FUN))
  if (backend$type == "fork") {
    res <- parallel::mclapply(X, FUN, mc.cores = backend$cores, mc.preschedule = FALSE)
  } else {
    cl <- parallel::makePSOCKcluster(backend$cores)
    on.exit(parallel::stopCluster(cl), add = TRUE)
    parallel::clusterEvalQ(cl, suppressPackageStartupMessages({
      library(dplyr); library(tibble); library(stats)
    }))
    parallel::clusterExport(cl, c(par_export_names(), export_extra), envir = globalenv())
    res <- parallel::parLapply(cl, X, FUN)
  }
  failed <- vapply(res, inherits, logical(1), what = "try-error")
  if (any(failed)) stop("parallel worker error: ", as.character(res[[which(failed)[1]]]))
  res
}

## Set PAR <- par_backend(1) to force serial execution and check the speedup.
PAR <- par_backend()
cat(sprintf("Parallel backend: %s on %d core(s).\n", PAR$type, PAR$cores))
## Parallel backend: psock on 8 core(s).

2 Model, conventions, and the four identifiability diagnostics

2.1 Observation model

At times \(0\le t_1<\cdots<t_K\le T\) the repeated measurements follow the time-specific latent-phase model

\[Y_{j,k}=\mathcal S(t_k,\Theta_j(t_k))+\varepsilon_{j,k},\qquad \Theta_j(t_k)\sim\varphi_{t_k},\qquad \varepsilon_{j,k}\sim\mathcal N(0,\sigma^2),\]

with the phases independent of the noise. All angular integrals use the Haar probability measure \(d\theta/(2\pi)\), so on the discrete grid every phase density satisfies row mean \(=1\) and the uniform law is the vector of ones.

Writing \(\widehat\varphi_t(n)=\int e^{-in\theta}\varphi_t(\theta)\,d\theta/(2\pi)\) and \(f_n(t)\) for the surface coefficients, the mixed moments \(m_n(t)=\mathbb E[Y(t)e^{-in\Theta(t)}]\) satisfy \(m_n=\sum_k f_k\widehat\varphi_t(n-k)\), i.e. \(M(z,t)=F(z,t)\Phi(z,t)\) on \(|z|=1\). Because \(F(e^{i\theta},t)=\mathcal S(t,\theta)\) and \(\Phi(e^{i\theta},t)=\varphi_t(\theta)\), this factorization is pointwise multiplication in \(\theta\): KPT is a deconvolution whose transfer function is the surface itself. That identity is the reason the four diagnostics below are all properties of \(\mathcal S\).

2.2 Identifiability

(D1) Rotational gauge. \((\mathcal S(t,\cdot-\alpha),\varphi_t(\cdot-\alpha))\) reproduces the law of \(Y\) for any \(\alpha(t)\). This is a model symmetry, so a gauge must be fixed before any comparison. We use the paper’s anchor: choose \(\alpha(t)=\arg\widehat\varphi_t(n_\star)/n_\star\) where \(n_\star\) is the smallest harmonic with non-negligible modulus, which makes \(e^{-in_\star\alpha}\widehat\varphi_t(n_\star)\) real and positive. The gauge is per time point, so a kimesurface assembled from independently anchored rows need not be smooth in \(t\); Section 2.5 treats that separately.

(D2) Rotational aliasing. With \(A_t=\{n:f_n(t)\ne0\}\) and \(d_t=\gcd A_t\), if \(d_t\ge2\) the surface is invariant under rotation by \(2\pi/d_t\) and \(\varphi_t\) is identifiable only modulo that finite subgroup.

(D3) Reflection symmetry. If \(\mathcal S(t,\alpha+u)=\mathcal S(t,\alpha-u)\) for all \(u\), equivalently if all \(f_n(t)e^{-in\alpha}\) are real, then \(\varphi_t\) and its reflection about \(\alpha\) induce the same law of \(Y\). This degeneracy is not removed by anchoring, it collapses the first-harmonic information matrix to rank one, and it invalidates \(\chi^2_2\) inference. Both V6 simulators had this property.

(D4) Spectral lower bound. Identifiability of the deconvolution requires \(\inf_\theta|\mathcal S(t,\theta)|\ge c_F>0\). Where the surface crosses zero the inverse problem is locally uninformative and the regularizer, not the data, determines \(\widehat\varphi\).

## ---- grids, quadrature, numerical guards ------------------------------------
theta_grid <- function(L) seq(0, 2*pi, length.out = L + 1L)[1L:L]
make_nvec  <- function(J) seq.int(-J, J, by = 1L)

freq_index <- function(L) {
  if (L %% 2L == 0L) c(0:(L/2L), -(L/2L - 1L):-1L)
  else               c(0:((L - 1L)/2L), -((L - 1L)/2L):-1L)
}

trapz_vec <- function(x, y) sum(0.5 * (y[-1L] + y[-length(y)]) * diff(x))

softmax_rows <- function(log_mat) {          # row-wise, overflow safe
  m <- apply(log_mat, 1L, max)
  w <- exp(log_mat - m)
  w / pmax(rowSums(w), .Machine$double.eps)
}

normalize_rows_prob <- function(P) {         # rows sum to 1 (probability weights)
  P <- pmax(P, 1e-15)
  sweep(P, 1L, rowSums(P), "/")
}

## ---- exact Euclidean projection onto {x >= 0, mean(x) = 1} -------------------
euclidean_simplex_project <- function(v, z = length(v)) {
  n <- length(v)
  u <- sort(v, decreasing = TRUE)
  cssv <- cumsum(u) - z
  rho <- max(which(u - cssv / seq_len(n) > 0))
  pmax(v - cssv[rho] / rho, 0)
}
row_project_simplex <- function(P) {
  P <- as.matrix(P)
  out <- t(apply(P, 1L, euclidean_simplex_project, z = ncol(P)))
  dimnames(out) <- dimnames(P)
  out
}

## ---- Fourier utilities on the phase circle ----------------------------------
discrete_fourier_coeffs <- function(x_theta, theta, nvec) {
  as.vector(drop(t(x_theta) %*% exp(-1i * outer(theta, nvec))) / length(theta))
}
synthesize_from_coeffs <- function(coeffs_n, theta, nvec) {
  as.vector(Re(drop(coeffs_n %*% exp(1i * outer(nvec, theta)))))
}
enforce_reality_coeffs <- function(coeffs_n, nvec) {
  coeffs_n <- as.complex(coeffs_n)
  coeffs_n[nvec == 0] <- Re(coeffs_n[nvec == 0])
  for (j in seq_len(max(abs(nvec)))) {
    pos <- which(nvec == j); neg <- which(nvec == -j)
    if (length(pos) == 1L && length(neg) == 1L) {
      avg <- 0.5 * (coeffs_n[pos] + Conj(coeffs_n[neg]))
      coeffs_n[pos] <- avg; coeffs_n[neg] <- Conj(avg)
    }
  }
  coeffs_n
}
coeffs_to_dft_order <- function(x_n, nvec) {
  N <- length(nvec); out <- complex(length = N)
  out[((nvec + N) %% N) + 1L] <- x_n
  out
}
dft_order_to_coeffs <- function(x_m, nvec) {
  N <- length(nvec)
  x_m[((nvec + N) %% N) + 1L]
}
coeffs_to_Xomega_fft <- function(x_n, nvec) fft(coeffs_to_dft_order(x_n, nvec))
Xomega_to_coeffs_fft <- function(X_w, nvec) {
  N <- length(nvec)
  dft_order_to_coeffs(fft(X_w, inverse = TRUE) / N, nvec)
}

## ---- periodic interpolation and exact spectral circular shift ---------------
periodic_interp <- function(theta, y, theta_new) {
  stats::approx(c(theta, 2*pi), c(y, y[1L]), xout = theta_new %% (2*pi),
                rule = 2L, ties = mean)$y
}
shift_periodic_row <- function(x, delta) {   # returns x(theta - delta)
  L <- length(x); n <- freq_index(L)
  Re(fft(fft(x) * exp(-1i * n * delta), inverse = TRUE) / L)
}
shift_periodic_matrix <- function(X, delta_vec) {
  stopifnot(nrow(X) == length(delta_vec))
  out <- X
  for (k in seq_len(nrow(X))) out[k, ] <- shift_periodic_row(X[k, ], delta_vec[k])
  out
}

2.3 Gauge fixing (D1), corrected

The V6 routine rotated by \(+\arg\big(\tfrac1L\sum\varphi(\theta)e^{+i\theta}\big)\). Since \(\tfrac1L\sum\varphi e^{+i\theta}=\widehat\varphi(-1)=\overline{\widehat\varphi(1)}\), that is a rotation by \(-\arg\widehat\varphi(1)\), which moves a density with mode \(\mu\) to mode \(2\mu\) rather than to \(0\). The corrected rule below is verified by unit test in Section 12.

## phihat(n) = mean(phi * exp(-i n theta)) under the Haar convention
phase_coeff <- function(row, theta, n) mean(row * exp(-1i * n * theta))

## Paper's rule: alpha = Arg(phihat(n_star)) / n_star, applied as phi(theta - alpha).
anchor_rows <- function(phi_mat, theta, harmonic = NULL, n_max = 8L, tol = 1e-3) {
  K <- nrow(phi_mat)
  delta <- numeric(K); n_used <- integer(K)
  for (k in seq_len(K)) {
    if (!is.null(harmonic)) {
      n_star <- as.integer(harmonic)
    } else {
      mods <- vapply(seq_len(n_max), function(n) Mod(phase_coeff(phi_mat[k, ], theta, n)), numeric(1))
      ## n_star = smallest harmonic whose modulus is above the noise floor; if the law is
      ## numerically uniform no preferred origin exists and we leave the row unrotated.
      n_star <- if (any(mods > tol)) which(mods > tol)[1L] else 0L
    }
    n_used[k] <- n_star
    delta[k] <- if (n_star >= 1L) Arg(phase_coeff(phi_mat[k, ], theta, n_star)) / n_star else 0
  }
  phi_rot <- row_project_simplex(shift_periodic_matrix(phi_mat, delta))
  list(phi = phi_rot, delta = delta, harmonic = n_used)
}
apply_anchor_to_surface <- function(S_mat, delta_vec) shift_periodic_matrix(S_mat, delta_vec)

## Convenience: anchor a (phi, S) pair with a single call.
anchor_pair <- function(phi_mat, S_mat, theta, harmonic = NULL) {
  a <- anchor_rows(phi_mat, theta, harmonic = harmonic)
  list(phi = a$phi, S = apply_anchor_to_surface(S_mat, a$delta), delta = a$delta)
}

2.4 The remaining three diagnostics (D2–D4)

  • (D2) rotational aliasing: d_t = gcd of the active harmonic indices.

  • (D3) reflection symmetry: smallest normalized asymmetry over all reflection axes.

  • (D4) spectral lower bound: (I2) holds only if \(\min|S|\) is bounded away from \(0.\)

gcd2 <- function(a, b) { while (b) { t <- b; b <- a %% b; a <- t }; abs(a) }

## (D2) rotational aliasing: d_t = gcd of the active harmonic indices
aliasing_diagnostic <- function(S_mat, theta, nvec, rel_tol = 0.02) {
  K <- nrow(S_mat); d <- integer(K)
  for (k in seq_len(K)) {
    f <- discrete_fourier_coeffs(S_mat[k, ], theta, nvec)
    act <- abs(nvec[Mod(f) > rel_tol * max(Mod(f))])
    act <- act[act > 0]
    d[k] <- if (!length(act)) 0L else Reduce(gcd2, act)
  }
  tibble(k = seq_len(K), d_t = d,
         identifiable_mod = ifelse(d > 1, sprintf("2*pi/%d", d), "full circle"))
}

## (D3) reflection symmetry: smallest normalized asymmetry over all reflection axes
reflection_diagnostic <- function(S_mat, theta, n_axis = 180L) {
  L <- length(theta); K <- nrow(S_mat)
  axes <- seq(0, pi, length.out = n_axis)
  u <- theta
  out <- numeric(K); best <- numeric(K)
  for (k in seq_len(K)) {
    s <- S_mat[k, ]; v <- stats::var(s)
    r <- vapply(axes, function(a) {
      lhs <- periodic_interp(theta, s, a + u)
      rhs <- periodic_interp(theta, s, a - u)
      mean((lhs - rhs)^2) / max(4 * v, 1e-15)
    }, numeric(1))
    out[k] <- min(r); best[k] <- axes[which.min(r)]
  }
  tibble(k = seq_len(K), asymmetry = out, axis = best,
         reflection_degenerate = out < 1e-6)
}

## (D4) spectral lower bound: (I2) holds only if min|S| is bounded away from 0
spectral_diagnostic <- function(S_mat) {
  tibble(k = seq_len(nrow(S_mat)),
         min_absS = apply(abs(S_mat), 1L, min),
         max_absS = apply(abs(S_mat), 1L, max)) |>
    mutate(cond_ratio = max_absS / pmax(min_absS, 1e-12),
           I2_ok = min_absS > 0.05 * max_absS)
}

## One call that returns a compact identifiability report for a surface.
identifiability_report <- function(S_mat, theta, nvec, label = "") {
  al <- aliasing_diagnostic(S_mat, theta, nvec)
  rf <- reflection_diagnostic(S_mat, theta)
  sp <- spectral_diagnostic(S_mat)
  tibble(
    Surface = label,
    `median d_t (D2)` = stats::median(al$d_t),
    `rows with d_t>1` = sum(al$d_t > 1),
    `median asymmetry (D3)` = signif(stats::median(rf$asymmetry), 3),
    `reflection-degenerate rows` = sum(rf$reflection_degenerate),
    `min over t of min|S| (D4)` = signif(min(sp$min_absS), 3),
    `rows failing (I2)` = sum(!sp$I2_ok)
  )
}

2.5 Cross-time registration of the surface

Because the gauge (D1) is fixed independently at each \(t_k\), the anchored rows of a kimesurface carry no common phase origin and the assembled two-dimensional object can be discontinuous in \(t\) even when the underlying surface is smooth. Anchoring is the right convention for comparing phase laws. It is the wrong convention for displaying or integrating a surface. In V8, we separate the two conventions. For surface work we register rows sequentially by maximizing the circular cross-correlation with the previous row, which selects the rotation that makes the surface as smooth in \(t\) as the data allow.

## Circular cross-correlation shift that best aligns 'row' to 'ref' (spectral, O(L log L)).
best_circular_shift <- function(ref, row) {
  L <- length(ref); n <- freq_index(L)
  cc <- Re(fft(Conj(fft(ref)) * fft(row), inverse = TRUE) / L)
  lag <- which.max(cc) - 1L
  (2*pi * lag) / L
}
register_surface_rows <- function(S_mat, phi_mat = NULL) {
  K <- nrow(S_mat); delta <- numeric(K)
  S_out <- S_mat
  for (k in seq_len(K)[-1L]) {
    d <- best_circular_shift(S_out[k - 1L, ], S_mat[k, ])
    delta[k] <- d
    S_out[k, ] <- shift_periodic_row(S_mat[k, ], -d)
  }
  phi_out <- if (is.null(phi_mat)) NULL else shift_periodic_matrix(phi_mat, -delta)
  list(S = S_out, phi = phi_out, delta = delta)
}

3 Kime-representation strategies

Each strategy below maps replicated longitudinal data \(Y\in\mathbb R^{N\times K}\) to a common object. A surface \(\widehat{\mathcal S}\) on the \((t_k,\theta_\ell)\) grid, a phase law \(\widehat\varphi\) on the same grid (row mean \(1\)), a noise scale \(\widehat\sigma\), and enough information to form a predictive distribution for a new replicate. This common interface makes the comparison in Section 8 meaningful. The strategies disagree about what the phase means, but they all imply a predictive law that can be scored on held-out replicates.

Strategy What it assumes What it estimates Cost
rep_uniform phases uniform; surface irrelevant beyond the induced marginal marginal of \(Y\) at each \(t\) \(O(NK)\)
rep_stitch replicate order carries the phase; phases equally spaced surface by epoch stitching \(O(NK\log N)\)
rep_transport \(\mathcal S(t,\cdot)\) monotone on a half-circle surface by rank/quantile transport \(O(NK\log N)\)
rep_laplace replicate structure is nuisance; the ensemble mean is the signal analytic kimesurface \(F=\mathcal Lf\), invertible \(O(K^2+K\,n_{\rm ILT})\)
rep_kpt the generative model itself \(\varphi_t\) and \(\mathcal S\) jointly, by deconvolution \(O(I_{\rm it}K(NL+J\log J))\)

3.1 Noise scale

V6 used 0.5 * mad(Y), an arbitrary scale that propagates into every posterior weight. V8 profiles the observed-data log-likelihood over \(\sigma\) given a working representation, and reports a sensitivity analysis in Section 8.

observed_loglik_gaussian <- function(Y, Sgrid, phi_grid, sigma) {
  K <- ncol(Y); ll <- 0
  for (k in seq_len(K)) {
    w <- phi_grid[k, ]; w <- pmax(w, 0); w <- w / sum(w)
    dens <- outer(Y[, k], Sgrid[k, ], function(y, mu) stats::dnorm(y, mu, sigma))
    ll <- ll + sum(log(pmax(drop(dens %*% w), 1e-300)))
  }
  ll
}

estimate_sigma_profile <- function(Y, Sgrid, phi_grid,
                                   grid = NULL, n_grid = 24L) {
  if (is.null(grid)) {
    s0 <- stats::sd(as.vector(Y))
    grid <- exp(seq(log(0.02 * s0), log(1.5 * s0), length.out = n_grid))
  }
  ll <- vapply(grid, function(s) observed_loglik_gaussian(Y, Sgrid, phi_grid, s), numeric(1))
  list(sigma = grid[which.max(ll)], grid = grid, loglik = ll)
}

3.2 A common constructor and predictive interface

Next, we define a generic kime-representation function, new_representation(), along with a predictive density of a new replicate, predictive_density(), a predictive CDF needed for PIT, predictive_cdf(), and a mechanism to draw from the predictive density at time index \(k\).

new_representation <- function(method, Sgrid, phi, sigma, theta, time, runtime = NA_real_, extra = list()) {
  stopifnot(nrow(Sgrid) == length(time), ncol(Sgrid) == length(theta))
  stopifnot(all(dim(phi) == dim(Sgrid)))
  structure(list(method = method, Sgrid = Sgrid, phi = phi, sigma = sigma,
                 theta = theta, time = time, runtime = runtime, extra = extra),
            class = "kime_representation")
}

## Predictive density of a new replicate at time index k, evaluated at y.
predictive_density <- function(rep, k, y) {
  w <- rep$phi[k, ]; w <- pmax(w, 0); w <- w / sum(w)
  dens <- outer(y, rep$Sgrid[k, ], function(yy, mu) stats::dnorm(yy, mu, rep$sigma))
  drop(dens %*% w)
}
## Predictive CDF at y (needed for PIT).
predictive_cdf <- function(rep, k, y) {
  w <- rep$phi[k, ]; w <- pmax(w, 0); w <- w / sum(w)
  cdf <- outer(y, rep$Sgrid[k, ], function(yy, mu) stats::pnorm(yy, mu, rep$sigma))
  drop(cdf %*% w)
}
## Draws from the predictive at time index k.
predictive_sample <- function(rep, k, B = 1000L) {
  w <- rep$phi[k, ]; w <- pmax(w, 0); w <- w / sum(w)
  idx <- sample.int(length(rep$theta), B, replace = TRUE, prob = w)
  rep$Sgrid[k, idx] + stats::rnorm(B, 0, rep$sigma)
}

3.3 Strategy 1: uniform-phase null model

The baseline reference other strategies are expected to outperform. The phase law is uniform by fiat. The kiemsurface (surface) is whatever reproduces the observed marginal. We use a band-limited quantile surface, so that pushing the uniform law through \(\widehat{\mathcal S}(t_k,\cdot)\) reproduces the empirical distribution of \(Y_{\cdot,k}\).

band_limit_row <- function(g, J) {
  L <- length(g); n <- freq_index(L)
  Re(fft(fft(g) * (abs(n) <= J), inverse = TRUE) / L)
}

rep_uniform <- function(Y, time, theta, J = 10L, sigma = NULL) {
  t0 <- proc.time()[["elapsed"]]
  K <- ncol(Y); N <- nrow(Y); L <- length(theta)
  S <- matrix(0, K, L)
  u <- (seq_len(N) - 0.5) / N                         # plotting positions
  for (k in seq_len(K)) {
    q <- stats::quantile(Y[, k], probs = u, names = FALSE, type = 7)
    ## place the sorted values around the circle, then band-limit
    S[k, ] <- band_limit_row(stats::approx(2*pi * u, q, xout = theta, rule = 2L)$y, J)
  }
  phi <- matrix(1, K, L)
  if (is.null(sigma)) sigma <- estimate_sigma_profile(Y, S, phi)$sigma
  new_representation("Uniform-null", S, phi, sigma, theta, time,
                     proc.time()[["elapsed"]] - t0)
}

3.4 Strategy 2: epoch stitching

The classical TCIU construction: the repetition index is the phase index. Replicate \(j\) at time \(t_k\) is assigned phase \(2\pi(j-1)/N\) in acquisition order, and the surface is the band-limited interpolant of the resulting scatter. This is faithful to how block-design data are actually collected, and it assumes that acquisition order is informative about phase, an assumption worth stating, because if the order is arbitrary the resulting surface is noise.

rep_stitch <- function(Y, time, theta, J = 10L, sigma = NULL, order_index = NULL) {
  t0 <- proc.time()[["elapsed"]]
  K <- ncol(Y); N <- nrow(Y); L <- length(theta)
  ph <- 2*pi * (seq_len(N) - 1) / N
  S <- matrix(0, K, L)
  for (k in seq_len(K)) {
    ord <- if (is.null(order_index)) seq_len(N) else order_index
    S[k, ] <- band_limit_row(periodic_interp(ph, Y[ord, k], theta), J)
  }
  phi <- matrix(1, K, L)
  if (is.null(sigma)) sigma <- estimate_sigma_profile(Y, S, phi)$sigma
  new_representation("Stitch", S, phi, sigma, theta, time,
                     proc.time()[["elapsed"]] - t0)
}

3.5 Strategy 3: rank / monotone transport

Sort the replicates and transport them onto the circle monotonically. If \(\mathcal S(t,\cdot)\) is unimodal, this recovers the surface up to reflection and rotation, which is exactly the ambiguity described by (D1) and (D3). It is the cheapest estimator that uses the values rather than the acquisition order.

rep_transport <- function(Y, time, theta, J = 10L, sigma = NULL) {
  t0 <- proc.time()[["elapsed"]]
  K <- ncol(Y); N <- nrow(Y); L <- length(theta)
  ph <- 2*pi * (seq_len(N) - 1) / N
  S <- matrix(0, K, L)
  for (k in seq_len(K)) {
    S[k, ] <- band_limit_row(periodic_interp(ph, sort(Y[, k]), theta), J)
  }
  phi <- matrix(1, K, L)
  if (is.null(sigma)) sigma <- estimate_sigma_profile(Y, S, phi)$sigma
  new_representation("Rank-transport", S, phi, sigma, theta, time,
                     proc.time()[["elapsed"]] - t0)
}

3.6 Strategy 4: analytic Laplace kimesurface with numerical inversion

The analytic route encodes the ensemble mean \(\bar y(t)\) as a holomorphic function on the kime plane, \(F(\kappa)=\int_0^\infty f(\tau)e^{-\kappa\tau}d\tau\), evaluated on the polar grid \(\kappa=te^{i\theta}\). Lerch uniqueness and Bromwich inversion make the map a bijection, so nothing is lost about the mean, but by construction the replicate-to-replicate phase structure is not represented. We therefore score this route on two things: the fidelity of its own round trip (which is what it guarantees) and, for comparability, the predictive law it implies, namely mean signal plus total replicate spread. This article (DOI: 10.1093/imamat/hxae033) provides background on Laplace transform for kimesurface reconstuction.

## Discrete Laplace transform of a sampled signal on an arbitrary complex grid.
discrete_laplace <- function(f_tau, tau, kappa) {
  dtau <- c(diff(tau), tail(diff(tau), 1L))
  vapply(kappa, function(z) sum(f_tau * exp(-z * tau) * dtau), complex(1))
}

## Valsa-Brancik numerical inverse Laplace transform.
## The alternating tail is summed by a genuine Euler transformation: M exact terms followed by a
## binomially weighted average of the next P partial sums. Averaging only the last two partial
## sums, as a naive implementation does, leaves an error of order 1e-2 on exp(-t); the Euler
## weights bring it to 1e-9 (self-test 9 checks this).
ilt_valsa <- function(Ffun, t, a = 8, n_main = 100L, n_euler = 29L) {
  w <- choose(n_euler, 0:n_euler) / 2^n_euler
  vapply(t, function(tt) {
    if (tt <= 0) return(NA_real_)
    n <- 0:(n_main + n_euler)
    z <- (a + 1i * (n + 0.5) * pi) / tt
    partial <- cumsum(((-1)^n) * Im(Ffun(z)))
    -exp(a) / tt * sum(w * partial[(n_main + 1L):(n_main + n_euler + 1L)])
  }, numeric(1))
}

rep_laplace <- function(Y, time, theta, J = 10L, a_ilt = 8, n_main = 100L, n_euler = 29L) {
  t0 <- proc.time()[["elapsed"]]
  K <- ncol(Y); L <- length(theta)
  mu_t <- colMeans(Y)
  tau <- time - min(time)
  ## Laplace kimesurface on the polar grid (stored for display and round-trip checks)
  t_pos <- pmax(tau, .Machine$double.eps)
  Fsurf <- matrix(0 + 0i, K, L)
  for (k in seq_len(K)) Fsurf[k, ] <- discrete_laplace(mu_t, tau, t_pos[k] * exp(1i * theta))
  ## round-trip: invert the transform of the mean signal back to the time domain
  Ffun <- function(z) discrete_laplace(mu_t, tau, z)
  mu_rt <- ilt_valsa(Ffun, pmax(tau, 1e-6), a = a_ilt, n_main = n_main, n_euler = n_euler)
  rt_err <- sqrt(mean((mu_rt - mu_t)^2, na.rm = TRUE)) / max(stats::sd(mu_t), 1e-12)
  ## predictive object: the mean signal carries no theta dependence
  S <- matrix(mu_t, nrow = K, ncol = L)   # row k is constant in theta
  phi <- matrix(1, K, L)
  sigma <- sqrt(mean(apply(Y, 2L, stats::var)))       # total replicate spread
  new_representation("Laplace-ILT", S, phi, sigma, theta, time,
                     proc.time()[["elapsed"]] - t0,
                     extra = list(F_surface = Fsurf, roundtrip_rel_rmse = rt_err,
                                  mu_roundtrip = mu_rt, mu_hat = mu_t))
}

3.7 Strategy 5: kime-phase tomography, safeguarded

The KPT update alternates an E-step (posterior phase weights), a regularized deconvolution of the phase law, exact projection onto the simplex, and optionally a surface update. V6 ran these updates unguarded. Replicating that scheme we observed the relative surface error grow from \(0.73\) at 20 iterations to \(1.86\) at 60, the iteration diverges, which is consistent with the paper’s own warning that the projected routines are not monotone EM.

V8 keeps the same updates but wraps them in a safeguard: the proposed iterate is damped toward the current one and accepted only if the observed-data log-likelihood does not decrease, with step halving otherwise. The resulting sequence is monotone by construction and stable in the iteration count. In the same configuration it reached relative surface error \(0.53\) and stayed there from 20 to 60 iterations.

initialize_surface <- function(Y, time, theta, J = 10L,
                               method = c("transport", "moment", "constant")) {
  method <- match.arg(method)
  K <- ncol(Y); L <- length(theta)
  if (method == "constant") return(matrix(mean(Y), K, L))
  if (method == "transport") {
    N <- nrow(Y); ph <- 2*pi * (seq_len(N) - 1) / N
    S <- matrix(0, K, L)
    for (k in seq_len(K)) S[k, ] <- band_limit_row(periodic_interp(ph, sort(Y[, k]), theta), J)
    return(S)
  }
  mu_k <- colMeans(Y); sd_k <- apply(Y, 2L, stats::sd)
  if (K >= 5L) {
    mu_k <- stats::smooth.spline(time, mu_k, spar = 0.6)$y
    sd_k <- stats::smooth.spline(time, sd_k, spar = 0.6)$y
  }
  outer(mu_k, rep(1, L)) + outer(sd_k, cos(theta)) + 0.25 * outer(sd_k, cos(2 * theta))
}

kpt_fit <- function(Y, time, theta,
                    J = 10L, lambda_phi = 1e-2, p_phi = 1,
                    lambda_F = 1e-2, p_F = 1,
                    update_surface = TRUE,
                    max_iter = 25L, tol = 1e-4,
                    damping = 0.5, safeguard = TRUE, max_halving = 6L,
                    sigma = NULL, init_surface = "transport",
                    verbose = FALSE) {
  t0 <- proc.time()[["elapsed"]]
  N <- nrow(Y); K <- ncol(Y); L <- length(theta)
  nvec <- make_nvec(J); Nn <- length(nvec)
  omega <- 2*pi * (0:(Nn - 1L)) / Nn
  Lam_phi <- (2 * sin(omega / 2))^(2 * p_phi)
  Lam_F   <- (2 * sin(omega / 2))^(2 * p_F)
  E_neg <- exp(-1i * outer(theta, nvec))

  S <- initialize_surface(Y, time, theta, J = J, method = init_surface)
  phi <- matrix(1, K, L)
  if (is.null(sigma)) sigma <- estimate_sigma_profile(Y, S, phi)$sigma

  ll_cur <- observed_loglik_gaussian(Y, S, phi, sigma)
  trace <- numeric(0); accepted <- logical(0); dphi_tr <- numeric(0)

  for (it in seq_len(max_iter)) {
    phi_new <- phi; S_new <- S
    for (k in seq_len(K)) {
      ## E-step: posterior phase weights for each replicate
      prior <- pmax(phi[k, ], 0); prior <- prior / sum(prior)
      logw <- outer(Y[, k], S[k, ], function(y, mu) -(y - mu)^2 / (2 * sigma^2))
      logw <- sweep(logw, 2L, log(pmax(prior, 1e-300)), "+")
      W <- softmax_rows(logw)
      m_hat <- colMeans(Y[, k] * (W %*% E_neg))            # mixed moments
      f_n <- enforce_reality_coeffs(discrete_fourier_coeffs(S[k, ], theta, nvec), nvec)

      ## phase update: regularized deconvolution then exact projection
      F_w <- coeffs_to_Xomega_fft(f_n, nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat, nvec)
      Phi_w <- Conj(F_w) * M_w / (Mod(F_w)^2 + lambda_phi * Lam_phi + 1e-10)
      phi_prop <- synthesize_from_coeffs(Xomega_to_coeffs_fft(Phi_w, nvec), theta, nvec)
      phi_new[k, ] <- euclidean_simplex_project(phi_prop, z = L)

      ## optional surface update, using the projected phase law
      if (isTRUE(update_surface)) {
        phi_n <- discrete_fourier_coeffs(phi_new[k, ], theta, nvec)
        P_w <- coeffs_to_Xomega_fft(phi_n, nvec)
        Fw2 <- Conj(P_w) * M_w / (Mod(P_w)^2 + lambda_F * Lam_F + 1e-10)
        f_upd <- enforce_reality_coeffs(Xomega_to_coeffs_fft(Fw2, nvec), nvec)
        S_new[k, ] <- synthesize_from_coeffs(f_upd, theta, nvec)
      }
    }

    ## safeguard: damped step, accepted only if the observed log-likelihood does not fall
    step <- 1; ok <- FALSE
    for (h in seq_len(max_halving)) {
      w <- step * damping
      S_try <- (1 - w) * S + w * S_new
      phi_try <- row_project_simplex((1 - w) * phi + w * phi_new)
      ll_try <- observed_loglik_gaussian(Y, S_try, phi_try, sigma)
      if (!safeguard || ll_try >= ll_cur - 1e-8) {
        dphi <- sqrt(mean((phi_try - phi)^2))
        S <- S_try; phi <- phi_try; ll_cur <- ll_try; ok <- TRUE
        break
      }
      step <- step / 2
    }
    trace <- c(trace, ll_cur); accepted <- c(accepted, ok)
    dphi_tr <- c(dphi_tr, if (ok) dphi else 0)
    if (verbose) message(sprintf("[KPT] iter %02d  loglik = %.2f  accepted = %s", it, ll_cur, ok))
    if (!ok) break                                   # no admissible step: stationary
    if (dphi < tol) break
  }

  list(S = S, phi = phi, sigma = sigma, theta = theta, time = time,
       loglik_trace = trace, accepted = accepted, delta_phi = dphi_tr,
       iterations = length(trace), runtime = proc.time()[["elapsed"]] - t0,
       monotone = all(diff(trace) >= -1e-6),
       parameters = list(J = J, lambda_phi = lambda_phi, lambda_F = lambda_F,
                         damping = damping, safeguard = safeguard,
                         update_surface = update_surface, init_surface = init_surface))
}

## Wrappers matching the paper's two named variants.
rep_kpt <- function(Y, time, theta, variant = c("FFT", "GEM"), ...) {
  variant <- match.arg(variant)
  fit <- kpt_fit(Y, time, theta, update_surface = identical(variant, "FFT"), ...)
  new_representation(paste0("KPT-", variant), fit$S, fit$phi, fit$sigma,
                     theta, time, fit$runtime, extra = fit)
}

The GEM variant holds the surface fixed at its initializer, so its surface error is a property of initialize_surface() rather than of the phase updates. V8 therefore also scores the initializer on its own (Initializer row in the tables) so that the two can be told apart, in V6 they could not.

4 Evaluation framework

Two families of criteria are used in the assessment that answer different questions.

  • Recovery metrics compare an estimate to a known truth and are available only in simulation.
  • Predictivemetrics score a representation on replicates it has never seen and are available on real data. V8 treats the predictive scores as primary, because they are the only criteria that transfer to an experiment where \(\varphi_t\) is unknown, which is every real experiment.
phase_metrics_df <- function(phi_true, phi_est, time, base2 = TRUE) {
  P <- normalize_rows_prob(phi_true); Q <- normalize_rows_prob(phi_est)
  lg <- if (base2) function(x) log2(pmax(x, 1e-15)) else function(x) log(pmax(x, 1e-15))
  K <- nrow(P)
  M <- 0.5 * (P + Q)
  tibble(
    time        = time,
    KL_true_est = rowSums(P * lg(P / Q)),
    KL_est_true = rowSums(Q * lg(Q / P)),
    JSD_bits    = 0.5 * rowSums(P * lg(P / M)) + 0.5 * rowSums(Q * lg(Q / M)),
    Hellinger   = sqrt(0.5 * rowSums((sqrt(P) - sqrt(Q))^2)),
    TV          = 0.5 * rowSums(abs(P - Q))
  )
}

## L2 error on the kime cone: the natural measure is t dt dtheta/(2 pi)
surface_L2_cone <- function(S_true, S_est, time) {
  stopifnot(all(dim(S_true) == dim(S_est)))
  err <- rowMeans((S_true - S_est)^2)
  ref <- rowMeans(S_true^2)
  abs_sq  <- trapz_vec(time, time * err)
  true_sq <- trapz_vec(time, time * ref)
  list(abs_L2 = sqrt(max(abs_sq, 0)),
       rel_L2 = sqrt(max(abs_sq, 0) / max(true_sq, 1e-15)),
       integrand = tibble(time = time, mse_theta = err))
}
## Sample-based CRPS, fair estimator: E|X - y| - 0.5 E|X - X'|.
## The second term is common to all observations at a time point, so it is computed once;
## it is evaluated in O(B log B) by the sorted-sample identity rather than by an outer product.
crps_self_term <- function(draws) {
  B <- length(draws); x <- sort(draws)
  ## E|X - X'| = (2 / B^2) * sum_i (2i - B - 1) x_(i)
  e <- 2 * sum((2 * seq_len(B) - B - 1) * x) / (B^2)
  0.5 * e * B / (B - 1)
}
crps_from_sample <- function(draws, y, self_term = NULL) {
  if (is.null(self_term)) self_term <- crps_self_term(draws)
  mean(abs(draws - y)) - self_term
}

## Held-out predictive scoring, shared by every representation.
score_representation <- function(rep, Y_test, B_draw = 500L, level = 0.90, seed = 1L) {
  set.seed(seed)
  K <- ncol(Y_test); N <- nrow(Y_test)
  lo_p <- (1 - level) / 2; hi_p <- 1 - lo_p
  ls <- crps <- cover <- pit <- vector("list", K)
  for (k in seq_len(K)) {
    y <- Y_test[, k]
    ls[[k]]  <- log(pmax(predictive_density(rep, k, y), 1e-300))
    pit[[k]] <- predictive_cdf(rep, k, y)
    dr <- predictive_sample(rep, k, B_draw)
    qs <- stats::quantile(dr, c(lo_p, hi_p), names = FALSE)
    cover[[k]] <- (y >= qs[1]) & (y <= qs[2])
    st <- crps_self_term(dr)
    crps[[k]] <- vapply(y, function(yy) crps_from_sample(dr, yy, self_term = st), numeric(1))
  }
  pit_v <- unlist(pit)
  tibble(
    method      = rep$method,
    logscore    = mean(unlist(ls)),
    CRPS        = mean(unlist(crps)),
    coverage    = mean(unlist(cover)),
    nominal     = level,
    PIT_KS      = suppressWarnings(as.numeric(stats::ks.test(pmin(pmax(pit_v, 1e-9), 1 - 1e-9), "punif")$statistic)),
    sigma_hat   = rep$sigma,
    runtime_sec = rep$runtime
  )
}

## Per-time-point log-score, used for paired comparisons between methods.
logscore_by_time <- function(rep, Y_test) {
  vapply(seq_len(ncol(Y_test)),
         function(k) mean(log(pmax(predictive_density(rep, k, Y_test[, k]), 1e-300))),
         numeric(1))
}

4.1 Uncertainty on the metrics

V6 reported point values, while V8 attaches a replicate bootstrap interval to every metric, and compares methods pairwise on the same time points, which removes the between-time variability that otherwise swamps the comparison.

## Replicate bootstrap for any representation. Indices are drawn in the master process so that the
## result is independent of the number of workers; `builder` is called as builder(Y, time, theta, ...).
## Pass the builder itself rather than a wrapper closure: a closure defined at top level does not carry
## the contents of the global environment to PSOCK workers.
# bootstrap_metric_ci <- function(Y_train, Y_test, time, theta, builder,
#                                 B = 60L, probs = c(0.025, 0.975), seed = 11L,
#                                 backend = PAR, B_draw = 200L, ...) {
#   N <- nrow(Y_train)
#   set.seed(seed)
#   idx_list <- lapply(seq_len(B), function(b) sample.int(N, N, replace = TRUE))
#   dots <- list(...)
#   one <- function(b) {
#     rp <- do.call(builder, c(list(Y_train[idx_list[[b]], , drop = FALSE], time, theta), dots))
#     score_representation(rp, Y_test, B_draw = B_draw, seed = b)
#   }
#   res <- bind_rows(par_lapply(seq_len(B), one, backend = backend))
bootstrap_metric_ci <- function(Y_train, Y_test, time, theta, builder,
                                B = 60L, probs = c(0.025, 0.975), seed = 11L,
                                backend = PAR, B_draw = 200L, ...) {
  ## Force all args so serialization ships values, not caller-scope promises.
  force(Y_train); force(Y_test); force(time); force(theta); force(builder)
  force(B_draw); force(probs); force(seed); force(backend)

  N <- nrow(Y_train)
  set.seed(seed)
  idx_list <- lapply(seq_len(B), function(b) sample.int(N, N, replace = TRUE))
  dots <- list(...)              # list(...) already forces ... elements
  # one <- function(b) {
  #   rp <- do.call(builder, c(list(Y_train[idx_list[[b]], , drop = FALSE], time, theta), dots))
  #   score_representation(rp, Y_test, B_draw = B_draw, seed = b)
  # }
  one <- local({
    ## Everything referenced by `one` becomes a plain value in this env.
    builder  <- builder
    idx_list <- idx_list
    Y_train  <- Y_train
    time     <- time
    theta    <- theta
    dots     <- dots
    Y_test   <- Y_test
    B_draw   <- B_draw
    function(b) {
      rp <- do.call(builder,
                    c(list(Y_train[idx_list[[b]], , drop = FALSE], time, theta), dots))
      score_representation(rp, Y_test, B_draw = B_draw, seed = b)
    }
  })
  res <- bind_rows(par_lapply(seq_len(B), one, backend = backend))
 
  tibble(method = res$method[1],
         logscore_lo = stats::quantile(res$logscore, probs[1], names = FALSE),
         logscore_hi = stats::quantile(res$logscore, probs[2], names = FALSE),
         CRPS_lo = stats::quantile(res$CRPS, probs[1], names = FALSE),
         CRPS_hi = stats::quantile(res$CRPS, probs[2], names = FALSE),
         B = B, cores = backend$cores)
}

## Paired comparison of two representations across time points, with BH control.
paired_method_test <- function(rep_a, rep_b, Y_test, alpha = 0.05) {
  la <- logscore_by_time(rep_a, Y_test)
  lb <- logscore_by_time(rep_b, Y_test)
  d <- la - lb
  tt <- stats::t.test(d)
  wl <- suppressWarnings(stats::wilcox.test(d))
  tibble(
    comparison = paste(rep_a$method, "vs", rep_b$method),
    mean_diff  = mean(d),
    ci_lo      = tt$conf.int[1], ci_hi = tt$conf.int[2],
    p_t        = tt$p.value,
    p_wilcoxon = wl$p.value
  )
}

bh_adjust <- function(tbl, p_col = "p_t") {
  tbl[[paste0(p_col, "_BH")]] <- stats::p.adjust(tbl[[p_col]], method = "BH")
  tbl
}

5 Statistical inference, and its calibration

5.1 Testing phase uniformity: a rank-aware Rayleigh statistic

At a fixed time the null is \(H_0:\varphi_t\equiv1\). The paper’s statistic de-biases the regularized first harmonic and forms \(R_N=N\,\tilde u^\top\widehat\Sigma^{-1}\tilde u\) with \(\tilde u=(\Re\tilde\varphi(1), \Im\tilde\varphi(1))\), claiming a \(\chi^2_2\) limit. That limit requires \(\widehat\Sigma\) to be nonsingular, and for a reflection-symmetric surface it is exactly singular. Every per-replicate contribution is real, so the imaginary part carries no information. Applying the \(\chi^2_2\) rule anyway rejects essentially always. Replicating that situation we obtained empirical size \(1.000\) at every nominal level.

V8 computes the effective rank of \(\widehat\Sigma\) by eigenvalue thresholding and refers the statistic to \(\chi^2_r\). Calibration is then correct in both regimes, which the next chunk verifies by simulation rather than assertion.

rayleigh_uniformity <- function(y, S_row, sigma, theta, nvec,
                                harmonic = 1L, tau = 1e-4, rank_tol = 1e-8) {
  N <- length(y)
  E_neg <- exp(-1i * outer(theta, nvec))
  ## posterior weights under the null (uniform prior)
  W <- softmax_rows(outer(y, S_row, function(yy, mu) -(yy - mu)^2 / (2 * sigma^2)))
  xi <- y * (W %*% E_neg)                             # N x Nn replicate contributions
  f_n <- discrete_fourier_coeffs(S_row, theta, nvec)
  F_w <- coeffs_to_Xomega_fft(f_n, nvec)
  Hinv <- Conj(F_w) / (Mod(F_w)^2 + tau)              # de-biased inverse filter
  Nn <- length(nvec)
  ord <- ((nvec + Nn) %% Nn) + 1L
  Xi <- matrix(0 + 0i, Nn, N)
  Xi[ord, ] <- t(xi)                                  # replicates in DFT order, one per column
  Cc <- mvfft(mvfft(Xi) * Hinv, inverse = TRUE) / Nn  # column-wise, vectorised over replicates
  c_j <- Cc[ord[which(nvec == harmonic)], ]
  u <- cbind(Re(c_j), Im(c_j))
  ubar <- colMeans(u)
  Sig  <- stats::cov(u)
  eg <- eigen(Sig, symmetric = TRUE)
  keep <- eg$values > rank_tol * max(max(eg$values), 1e-300)
  if (!any(keep)) return(list(stat = 0, df = 0L, p_value = 1, rank_ratio = 0))
  z <- drop(crossprod(eg$vectors[, keep, drop = FALSE], ubar))
  stat <- N * sum(z^2 / eg$values[keep])
  list(stat = stat, df = sum(keep),
       p_value = stats::pchisq(stat, df = sum(keep), lower.tail = FALSE),
       rank_ratio = min(eg$values) / max(eg$values),
       phi1 = mean(c_j))
}

## Parametric bootstrap calibration under H0, as recommended in the paper.
rayleigh_bootstrap_p <- function(y, S_row, sigma, theta, nvec, B = 199L, seed = 5L, ...) {
  set.seed(seed)
  obs <- rayleigh_uniformity(y, S_row, sigma, theta, nvec, ...)$stat
  N <- length(y); L <- length(theta)
  null_stat <- vapply(seq_len(B), function(b) {
    idx <- sample.int(L, N, replace = TRUE)           # uniform phases
    yb <- S_row[idx] + stats::rnorm(N, 0, sigma)
    rayleigh_uniformity(yb, S_row, sigma, theta, nvec, ...)$stat
  }, numeric(1))
  (1 + sum(null_stat >= obs)) / (B + 1)
}
## Monte Carlo validation: size under H0 and power under von Mises alternatives.
dvonmises <- function(theta, mu, kappa) exp(kappa * cos(theta - mu)) / (2*pi * besselI(kappa, 0))
vm_density_row <- function(theta, mu, kappa) { d <- dvonmises(theta, mu, kappa); d / mean(d) }

validate_rayleigh <- function(S_row, theta, nvec, sigma = 0.2, N = 150L,
                              B = CFG$B_calib, kappas = if (FAST) c(0, 0.5, 1, 2) else c(0, 0.25, 0.5, 1, 2),
                              levels = c(0.10, 0.05, 0.01), seed = 99L) {
  set.seed(seed)
  L <- length(theta)
  run <- function(kappa) {
    p <- vapply(seq_len(B), function(b) {
      idx <- if (kappa == 0) sample.int(L, N, replace = TRUE)
             else sample.int(L, N, replace = TRUE,
                             prob = vm_density_row(theta, 1.0, kappa) / L)
      y <- S_row[idx] + stats::rnorm(N, 0, sigma)
      rayleigh_uniformity(y, S_row, sigma, theta, nvec)$p_value
    }, numeric(1))
    p
  }
  res <- lapply(kappas, run); names(res) <- paste0("kappa=", kappas)
  size_tbl <- bind_rows(lapply(names(res), function(nm) {
    p <- res[[nm]]
    tibble(scenario = nm,
           level = levels,
           rejection_rate = vapply(levels, function(a) mean(p <= a), numeric(1)),
           mc_se = sqrt(levels * (1 - levels) / length(p)))
  }))
  list(p_values = res, table = size_tbl)
}

5.2 Comparing phase laws across times or conditions

Rotation-invariant Sobolev discrepancies compare two phase laws without re-fixing the gauge \[D_w(\varphi,\psi)=\sum_{n\ne0}w_n|\widehat\varphi(n)-\widehat\psi(n)|^2\] with decaying weights, \(w_n\). Because the gauge is per time point, we also provide the rotation-invariant variant built from \(|\widehat\varphi(n)|\), which is what should be used when the two laws were anchored independently.

sobolev_discrepancy <- function(phi_a, phi_b, theta, n_max = 8L, s = 1,
                                rotation_invariant = TRUE) {
  w <- (seq_len(n_max))^(-2 * s)
  ca <- vapply(seq_len(n_max), function(n) phase_coeff(phi_a, theta, n), complex(1))
  cb <- vapply(seq_len(n_max), function(n) phase_coeff(phi_b, theta, n), complex(1))
  if (rotation_invariant) sum(w * (Mod(ca) - Mod(cb))^2)
  else                    sum(w * Mod(ca - cb)^2)
}

## Permutation test for a two-group difference in phase law (e.g. ON vs OFF blocks).
permutation_phase_test <- function(phi_mat, group, theta, n_perm = 499L,
                                   n_max = 8L, s = 1, seed = 7L) {
  set.seed(seed)
  g <- as.factor(group); stopifnot(nlevels(g) == 2L)
  agg <- function(idx) colMeans(phi_mat[idx, , drop = FALSE])
  obs <- sobolev_discrepancy(agg(which(g == levels(g)[1])),
                             agg(which(g == levels(g)[2])), theta, n_max, s)
  null <- vapply(seq_len(n_perm), function(b) {
    gp <- sample(g)
    sobolev_discrepancy(agg(which(gp == levels(g)[1])),
                        agg(which(gp == levels(g)[2])), theta, n_max, s)
  }, numeric(1))
  list(statistic = obs, p_value = (1 + sum(null >= obs)) / (n_perm + 1), null = null)
}

5.3 Simultaneous bootstrap bands and their empirical coverage

The paper states a bootstrap band result on the discrete phase grid. V8 implements the band and then checks its coverage by simulation, which is the part that was missing.

# bootstrap_phase_bands <- function(Y, time, theta, builder, B = 60L,
#                                   level = 0.90, seed = 3L, backend = PAR, ...) {
#   N <- nrow(Y)
#   set.seed(seed)
#   idx_list <- lapply(seq_len(B), function(b) sample.int(N, N, replace = TRUE))
#   dots <- list(...)
#   one <- function(b) {
#     rp <- do.call(builder, c(list(Y[idx_list[[b]], , drop = FALSE], time, theta), dots))
#     anchor_rows(rp$phi, theta)$phi
#   }
#   # one <- local({
#   #   ## Everything referenced by `one` becomes a plain value in this env.
#   #   builder  <- builder
#   #   idx_list <- idx_list
#   #   Y_train  <- Y_train
#   #   time     <- time
#   #   theta    <- theta
#   #   dots     <- dots
#   #   Y_test   <- Y_test
#   #   B_draw   <- B_draw
#   #   function(b) {
#   #     rp <- do.call(builder,
#   #                   c(list(Y_train[idx_list[[b]], , drop = FALSE], time, theta), dots))
#   #     score_representation(rp, Y_test, B_draw = B_draw, seed = b)
#   #   }
#   # })
#   
#   arr <- simplify2array(par_lapply(seq_len(B), one, backend = backend))   # K x L x B
#   centre <- apply(arr, c(1, 2), stats::median)
#   sdev   <- apply(arr, c(1, 2), stats::sd)
#   ## simultaneous width: sup-t calibration over the grid
#   zb <- vapply(seq_len(B), function(b) max(abs(arr[, , b] - centre) / pmax(sdev, 1e-12)), numeric(1))
#   qz <- stats::quantile(zb, level, names = FALSE)
#   list(centre = centre, sd = sdev, crit = qz,
#        lower = centre - qz * sdev, upper = centre + qz * sdev,
#        level = level, B = B, cores = backend$cores)
# }
bootstrap_phase_bands <- function(Y, time, theta, builder, B = 60L,
                                  level = 0.90, seed = 3L, backend = PAR, ...) {
  N <- nrow(Y)
  set.seed(seed)
  idx_list <- lapply(seq_len(B), function(b) sample.int(N, N, replace = TRUE))
  dots <- list(...)

  one <- local({
    Y        <- Y
    idx_list <- idx_list
    time     <- time
    theta    <- theta
    builder  <- builder
    dots     <- dots
    function(b) {
      rp <- do.call(builder,
                    c(list(Y[idx_list[[b]], , drop = FALSE], time, theta), dots))
      anchor_rows(rp$phi, theta)$phi
    }
  })

  arr <- simplify2array(par_lapply(seq_len(B), one, backend = backend))
  centre <- apply(arr, c(1, 2), stats::median)
  sdev   <- apply(arr, c(1, 2), stats::sd)
  zb <- vapply(seq_len(B),
               function(b) max(abs(arr[, , b] - centre) / pmax(sdev, 1e-12)),
               numeric(1))
  qz <- stats::quantile(zb, level, names = FALSE)
  list(centre = centre, sd = sdev, crit = qz,
       lower = centre - qz * sdev, upper = centre + qz * sdev,
       level = level, B = B, cores = backend$cores)
}

6 Simulators

Both simulators expose two switches that turn the V6 pathologies on and off, so that every claim in Section 1 can be reproduced rather than taken on faith.

  • well_posed = TRUE adds a DC offset large enough that \(\min_\theta|\mathcal S(t,\theta)|>0\), so condition (I2) holds. With FALSE the surface crosses zero, as in V6.
  • break_reflection = TRUE adds a genuine \(\sin\) harmonic, so the surface is not even and the phase law is identifiable beyond reflection. With FALSE the surface is even, as in V6.

Setting both to FALSE reproduces the V6 generative model exactly.

canonical_hrf <- function(t) {
  a1 <- 6; a2 <- 16; b1 <- 1; b2 <- 1; cc <- 1/6
  h <- (t^(a1 - 1) * b1^a1 * exp(-b1 * t)) / gamma(a1) -
       cc * (t^(a2 - 1) * b2^a2 * exp(-b2 * t)) / gamma(a2)
  h[t < 0] <- 0
  h
}

simulate_fmri <- function(N = CFG$N_rep, K = CFG$K_fmri, TR = 2,
                          L_theta = CFG$L_theta, noise_sd = 0.2,
                          kappa_on = 6, kappa_off = 1, mu_amp = 1.2,
                          well_posed = TRUE, break_reflection = TRUE,
                          n_split = 2L, seed = 2026) {
  set.seed(seed)
  theta <- theta_grid(L_theta)
  time <- seq(0, by = TR, length.out = K)
  block_len <- max(1L, round(30 / TR))
  design <- rep_len(rep(c(0, 1), each = block_len), K)

  # hrf <- canonical_hrf(seq(0, 30, by = TR))
  # bold <- as.numeric(stats::filter(design, hrf, method = "convolution", sides = 1))
  hrf <- canonical_hrf(seq(0, 30, by = TR))
  hrf <- hrf[seq_len(min(length(hrf), length(design)))]
  bold <- as.numeric(stats::filter(design, hrf, 
                                   method = "convolution", sides = 1))
  
  bold[is.na(bold)] <- 0
  bold <- if (diff(range(bold)) > 0) 0.2 + 0.8 * (bold - min(bold)) / diff(range(bold)) else rep(0.6, K)

  baseline <- 0.05 * sin(2*pi * time / max(time))
  dc <- if (well_posed) 1.6 else 0
  c_sin <- if (break_reflection) 0.35 else 0
  mu_theta <- mu_amp * sin(2*pi * time / max(time))

  S_true <- matrix(0, K, L_theta); phi_true <- matrix(0, K, L_theta)
  for (k in seq_len(K)) {
    S_true[k, ] <- dc + baseline[k] + bold[k] * cos(theta) +
                   0.25 * bold[k] * cos(2 * theta) + c_sin * bold[k] * sin(2 * theta)
    phi_true[k, ] <- vm_density_row(theta, mu_theta[k],
                                    if (design[k] == 1) kappa_on else kappa_off)
  }

  draw <- function(n) {
    Y <- matrix(0, n, K)
    for (k in seq_len(K)) {
      idx <- sample.int(L_theta, n, replace = TRUE, prob = phi_true[k, ] / L_theta)
      Y[, k] <- S_true[k, idx] + stats::rnorm(n, 0, noise_sd)
    }
    Y
  }
  splits <- lapply(seq_len(n_split), function(i) draw(N))

  list(Y = splits[[1]], Y_test = if (n_split > 1) splits[[2]] else NULL,
       time = time, theta = theta, design = design, bold = bold,
       S_true = S_true, phi_true = phi_true, sigma = noise_sd,
       settings = list(well_posed = well_posed, break_reflection = break_reflection,
                       kappa_on = kappa_on, kappa_off = kappa_off, N = N, K = K))
}

simulate_double_slit <- function(N = CFG$N_rep, K = 50L, L_theta = CFG$L_theta,
                                 kappa_phase = 4, noise_sd = 0.15, kappa_wave = 12,
                                 well_posed = TRUE, break_reflection = TRUE,
                                 n_split = 2L, seed = 2025) {
  set.seed(seed)
  theta <- theta_grid(L_theta); time <- seq(0, 1, length.out = K)
  dc <- if (well_posed) 1.4 else 0
  c_sin <- if (break_reflection) 0.3 else 0
  S_true <- matrix(0, K, L_theta); phi_true <- matrix(0, K, L_theta)
  for (k in seq_len(K)) {
    tt <- time[k]
    env <- 0.8 * exp(-20 * (tt - 0.5)^2)
    S_true[k, ] <- dc + env * cos(theta) + 0.3 * env * cos(2 * theta) +
                   c_sin * env * sin(2 * theta)
    d <- 0.5 * dvonmises(theta, -0.6 + 2*pi*0.1*tt, kappa_phase) +
         0.5 * dvonmises(theta,  1.1 + 2*pi*0.1*tt, kappa_phase)
    phi_true[k, ] <- d / mean(d)
  }
  draw <- function(n) {
    Y <- matrix(0, n, K)
    for (k in seq_len(K)) {
      idx <- sample.int(L_theta, n, replace = TRUE, prob = phi_true[k, ] / L_theta)
      Y[, k] <- S_true[k, idx] + stats::rnorm(n, 0, noise_sd)
    }
    Y
  }
  splits <- lapply(seq_len(n_split), function(i) draw(N))
  list(Y = splits[[1]], Y_test = if (n_split > 1) splits[[2]] else NULL,
       time = time, theta = theta, S_true = S_true, phi_true = phi_true,
       sigma = noise_sd,
       settings = list(well_posed = well_posed, break_reflection = break_reflection))
}

7 The identifiability diagnostics on V6 vs. V8 generative settings

nvec_diag <- make_nvec(CFG$J)
sim_v6 <- simulate_fmri(well_posed = FALSE, break_reflection = FALSE, seed = 4001)
sim_v8 <- simulate_fmri(well_posed = TRUE,  break_reflection = TRUE,  seed = 4001)

diag_tbl <- bind_rows(
  identifiability_report(sim_v6$S_true, sim_v6$theta, nvec_diag, "V6 setting (even, zero-crossing)"),
  identifiability_report(sim_v8$S_true, sim_v8$theta, nvec_diag, "V8 setting (generic, DC-offset)")
)
knitr::kable(diag_tbl, caption = "Identifiability diagnostics (D2-D4) for the two generative settings.")
Table 7.1: Identifiability diagnostics (D2-D4) for the two generative settings.
Surface median d_t (D2) rows with d_t>1 median asymmetry (D3) reflection-degenerate rows min over t of min|S| (D4) rows failing (I2)
V6 setting (even, zero-crossing) 1 0 0.0000 150 1.14e-05 150
V8 setting (generic, DC-offset) 1 0 0.0751 0 4.83e-01 0

The V6 row shows the two structural problems together: every surface row is reflection-degenerate, so the phase law is identifiable only up to \(\theta\mapsto-\theta\), and \(\min_\theta|\mathcal S|\) is numerically zero, so condition (I2) fails and the deconvolution is locally uninformative. The V8 row is clean on both counts. Neither condition is visible in any reconstruction metric, which is why they must be checked directly.

k_show <- round(nrow(sim_v6$S_true) * 0.35)
bind_rows(
  tibble(theta = sim_v6$theta, S = sim_v6$S_true[k_show, ], Setting = "V6 (even, crosses zero)"),
  tibble(theta = sim_v8$theta, S = sim_v8$S_true[k_show, ], Setting = "V8 (generic, bounded away from zero)")
) |>
  ggplot(aes(theta, S, colour = Setting)) +
  geom_hline(yintercept = 0, linewidth = 0.3, colour = "grey50") +
  geom_line(linewidth = 0.8) +
  scale_x_continuous(breaks = c(0, pi, 2*pi), labels = c("0", expression(pi), expression(2*pi))) +
  labs(title = "Surface slice at a fixed time",
       subtitle = "The V6 slice is symmetric about theta = 0 and touches zero; both break identifiability",
       x = expression(theta), y = expression(S(t[k], theta))) +
  theme(legend.position = "bottom")

8 Head-to-head comparison of kime-representation strategies

sim <- sim_v8
th <- sim$theta; tm <- sim$time; Ytr <- sim$Y; Yte <- sim$Y_test
Jc <- CFG$J

reps <- list()
reps[["Uniform-null"]]   <- rep_uniform(Ytr, tm, th, J = Jc)
reps[["Stitch"]]         <- rep_stitch(Ytr, tm, th, J = Jc)
reps[["Rank-transport"]] <- rep_transport(Ytr, tm, th, J = Jc)
reps[["Laplace-ILT"]]    <- rep_laplace(Ytr, tm, th, J = Jc)
reps[["KPT-GEM"]]        <- rep_kpt(Ytr, tm, th, variant = "GEM", J = Jc,
                                    max_iter = CFG$max_iter, sigma = sim$sigma)
reps[["KPT-FFT"]]        <- rep_kpt(Ytr, tm, th, variant = "FFT", J = Jc,
                                    max_iter = CFG$max_iter, sigma = sim$sigma)
## the initializer, scored on its own so that KPT-GEM can be interpreted
S_init <- initialize_surface(Ytr, tm, th, J = Jc, method = "transport")
reps[["Initializer"]]    <- new_representation("Initializer", S_init,
                                               matrix(1, length(tm), length(th)),
                                               sim$sigma, th, tm, 0)
pred_tbl <- bind_rows(lapply(reps, score_representation, Y_test = Yte)) |>
  arrange(desc(logscore))
# knitr::kable(pred_tbl, digits = 4,
#              caption = "Held-out predictive scores (higher log-score is better; lower CRPS is better; coverage should match the nominal 0.90; PIT_KS near 0 indicates calibrated predictive distributions).")
# knitr::kable(
#   recovery, 
#   digits = 4,
#   col.names = c(
#     "method", 
#     "mean_JSD_bits $\\downarrow$", 
#     "mean_Hellinger $\\downarrow$", 
#     "mean_TV $\\downarrow$", 
#     "surface_rel_L2 $\\downarrow$"
#   ),
#   caption = "Recovery against the known truth after per-time gauge anchoring. Phase metrics for the uniform-phase strategies measure the divergence between the truth and a flat law, and are included as the reference every phase estimator must beat."
# )

# <!-- **Format Adaptations** -->
# 
# <!-- * **PDF / LaTeX Output:** The code above uses LaTeX math syntax `$\\downarrow$` since all these metrics represent errors, distances, or divergences where lower values indicate closer recovery to the truth. -->
# <!-- * **HTML / Word Output:** For HTML or Word documents, swap the LaTeX string for Unicode arrows: -->
# col.names = c("method", "mean_JSD_bits ↓", "mean_Hellinger ↓", "mean_TV ↓", "surface_rel_L2 ↓")

knitr::kable(
  pred_tbl, 
  digits = 4,
  col.names = c(
    "method", 
    "logscore $\\uparrow$", 
    "CRPS $\\downarrow$", 
    "coverage", 
    "nominal", 
    "PIT_KS $\\approx 0$", 
    "sigma_hat", 
    "runtime_sec"
  ),
  caption = "Held-out predictive scores (higher log-score is better; lower CRPS is better; coverage should match the nominal 0.90; PIT_KS near 0 indicates calibrated predictive distributions)."
)
Table 8.1: Held-out predictive scores (higher log-score is better; lower CRPS is better; coverage should match the nominal 0.90; PIT_KS near 0 indicates calibrated predictive distributions).
method logscore \(\uparrow\) CRPS \(\downarrow\) coverage nominal PIT_KS \(\approx 0\) sigma_hat runtime_sec
KPT-GEM -0.3836 0.2166 0.9016 0.9 0.0075 0.2000 89.25
Rank-transport -0.3878 0.2166 0.8952 0.9 0.0045 0.1167 65.11
Uniform-null -0.3883 0.2166 0.8923 0.9 0.0061 0.1167 14.04
KPT-FFT -0.3925 0.2182 0.9082 0.9 0.0517 0.2000 60.92
Initializer -0.4012 0.2173 0.9337 0.9 0.0296 0.2000 0.00
Stitch -0.4787 0.2223 0.9021 0.9 0.0178 0.3599 13.44
Laplace-ILT -0.5143 0.2244 0.8995 0.9 0.0260 0.4029 3.46 
truth_anchor <- anchor_pair(sim$phi_true, sim$S_true, th)
recovery <- bind_rows(lapply(names(reps), function(nm) {
  r <- reps[[nm]]
  a <- anchor_pair(r$phi, r$Sgrid, th)
  pm <- phase_metrics_df(truth_anchor$phi, a$phi, tm)
  sl <- surface_L2_cone(truth_anchor$S, a$S, tm)
  tibble(method = nm, mean_JSD_bits = mean(pm$JSD_bits),
         mean_Hellinger = mean(pm$Hellinger), mean_TV = mean(pm$TV),
         surface_rel_L2 = sl$rel_L2)
})) |> arrange(surface_rel_L2)

# knitr::kable(recovery, digits = 4,
#              caption = "Recovery against the known truth after per-time gauge anchoring. Phase metrics for the uniform-phase strategies measure the divergence between the truth and a flat law, and are included as the reference every phase estimator must beat.")

knitr::kable(recovery, digits = 4, col.names = c(
    "method", 
    "mean_JSD_bits $\\downarrow$", 
    "mean_Hellinger $\\downarrow$", 
    "mean_TV $\\downarrow$", 
    "surface_rel_L2 $\\downarrow$"
  ),
  caption = "Recovery against the known truth after per-time gauge anchoring. Phase metrics for the uniform-phase strategies measure the divergence between the truth and a flat law, and are included as the reference every phase estimator must beat."
)
Table 8.2: Recovery against the known truth after per-time gauge anchoring. Phase metrics for the uniform-phase strategies measure the divergence between the truth and a flat law, and are included as the reference every phase estimator must beat.
method mean_JSD_bits \(\downarrow\) mean_Hellinger \(\downarrow\) mean_TV \(\downarrow\) surface_rel_L2 \(\downarrow\)
KPT-FFT 0.2039 0.3630 0.3930 0.3280
Laplace-ILT 0.2759 0.4352 0.4856 0.4174
Stitch 0.2759 0.4352 0.4856 0.4271
KPT-GEM 0.2380 0.3979 0.4283 0.4357
Uniform-null 0.2759 0.4352 0.4856 0.5042
Rank-transport 0.2759 0.4352 0.4856 0.5060
Initializer 0.2759 0.4352 0.4856 0.5060 
t_boot <- system.time({
  boot_tbl <- bind_rows(
    bootstrap_metric_ci(Ytr, Yte, tm, th, rep_uniform,   B = CFG$B_boot, J = Jc),
    bootstrap_metric_ci(Ytr, Yte, tm, th, rep_transport, B = CFG$B_boot, J = Jc),
    bootstrap_metric_ci(Ytr, Yte, tm, th, rep_kpt, B = max(8L, CFG$B_boot %/% 6L),
                        variant = "FFT", J = Jc, max_iter = 10L, sigma = sim$sigma)
  )
})
knitr::kable(boot_tbl, digits = 4,
             caption = "Replicate-bootstrap 95% intervals for the predictive scores.")
Table 8.3: Replicate-bootstrap 95% intervals for the predictive scores.
method logscore_lo logscore_hi CRPS_lo CRPS_hi B cores
Uniform-null -0.4282 -0.4029 0.2170 0.2182 200 8
Rank-transport -0.4282 -0.4027 0.2170 0.2182 200 8
KPT-FFT -0.4028 -0.3970 0.2182 0.2194 33 8 
cat(sprintf("elapsed %.1f s on %d core(s)\n", t_boot[["elapsed"]], PAR$cores))
## elapsed 5439.8 s on 8 core(s)
paired <- bind_rows(
  paired_method_test(reps[["KPT-FFT"]],        reps[["Uniform-null"]], Yte),
  paired_method_test(reps[["KPT-FFT"]],        reps[["Rank-transport"]], Yte),
  paired_method_test(reps[["Rank-transport"]], reps[["Uniform-null"]], Yte),
  paired_method_test(reps[["Stitch"]],         reps[["Uniform-null"]], Yte),
  paired_method_test(reps[["Laplace-ILT"]],    reps[["Uniform-null"]], Yte)
) |> bh_adjust("p_t")
knitr::kable(paired, digits = 4,
             caption = "Paired comparisons of per-time-point held-out log-score, with Benjamini-Hochberg adjustment across the five comparisons.")
Table 8.4: Paired comparisons of per-time-point held-out log-score, with Benjamini-Hochberg adjustment across the five comparisons.
comparison mean_diff ci_lo ci_hi p_t p_wilcoxon p_t_BH
KPT-FFT vs Uniform-null -0.0042 -0.0073 -0.0011 0.0083 0.0001 0.0103
KPT-FFT vs Rank-transport -0.0047 -0.0075 -0.0019 0.0011 0.0000 0.0019
Rank-transport vs Uniform-null 0.0005 -0.0002 0.0013 0.1838 0.4852 0.1838
Stitch vs Uniform-null -0.0905 -0.1036 -0.0774 0.0000 0.0000 0.0000
Laplace-ILT vs Uniform-null -0.1261 -0.1468 -0.1053 0.0000 0.0000 0.0000 
on_idx <- which(sim$design == 1); off_idx <- which(sim$design == 0)
cond_tbl <- bind_rows(lapply(names(reps), function(nm) {
  a <- anchor_pair(reps[[nm]]$phi, reps[[nm]]$Sgrid, th)
  pm <- phase_metrics_df(truth_anchor$phi, a$phi, tm)
  tibble(method = nm,
         JSD_ON  = mean(pm$JSD_bits[on_idx]),
         JSD_OFF = mean(pm$JSD_bits[off_idx]))
}))
knitr::kable(cond_tbl, digits = 4,
             caption = "Phase divergence split by block. The OFF blocks have kappa = 1, a nearly flat phase law, so the uniform null is hard to beat there; the ON blocks (kappa = 6) are where a phase estimator can earn its keep.")
Table 8.5: Phase divergence split by block. The OFF blocks have kappa = 1, a nearly flat phase law, so the uniform null is hard to beat there; the ON blocks (kappa = 6) are where a phase estimator can earn its keep.
method JSD_ON JSD_OFF
Uniform-null 0.4739 0.0778
Stitch 0.4739 0.0778
Rank-transport 0.4739 0.0778
Laplace-ILT 0.4739 0.0778
KPT-GEM 0.4122 0.0638
KPT-FFT 0.3648 0.0431
Initializer 0.4739 0.0778

9 KPT benefits

The Jensen-Shannon Divergence (\(\mathrm{JSD}\)) is a symmetric, smoothed measure of the similarity between two probability distributions. Because it is bounded between \(0\) and \(1\) (when using base-\(2\) logarithms), it provides a stable distance metric.

Simply averaging a divergence over all time points answers the wrong question. In fMRI, half the time points are rest blocks whose true phase law is nearly uniform and there the null model is nearly correct by construction and no estimator can gain. Thus, including those times in an average dilutes a real effect with a set of times where no effect is possible.

In V8, 2 revisions address that problem.

  • Report a skill score relative to the null, \[skill = 1-\mathrm{JSD}_{\text{method}}/\mathrm{JSD}_{\text{null}},\] which is positive exactly when the method beats the null.
  • Stratifying by how far the truth actually is from uniform, as measured by the mean resultant length \[R(t)=|\widehat\varphi_t(1)|,\] which is \(0\) for the uniform law and approaches \(1\) as the law concentrates.

KPT inference is paired across time points, so the between-time variability that dominates, the raw averages cancels.

The skill measures the relative performance of a predictive method compared to a baseline (null) model using the Jensen-Shannon Divergence (\(\mathrm{JSD}\)), which measures distributional error, or distance with lower JSD is better. The skill formulation scales the error reduction relative to the baseline.

  • Perfect skill, \(skill = 1\), indicates that the method has a \(\mathrm{JSD}_{\text{method}}=0\), meaning its predicted distribution perfectly matches the reference distribution.
  • Positive skill, \(0 < skill < 1\), suggests that the method outperformes the null. The closer the skill is to \(1\), the greater the reduction in divergence. For example, a skill score of \(0.6\) means the method reduced the baseline’s divergence error by \(60\%\).
  • Trivial skill, \(skill = 0\), when the method performs identically to the null model (\(\mathrm{JSD}_{\text{method}} = \mathrm{JSD}_{\text{null}}\)). There is no accuracy benefit to using the new method over the baseline.
  • Negative skill, \(skill < 0\), the method has a higher divergence than the null model (\(\mathrm{JSD}_{\text{method}} > \mathrm{JSD}_{\text{null}}\)). This indicates that the method is actively performing worse than a naive baseline (e.g., random guessing, historical mean, or persistence).
K_all <- length(tm)
phi_null <- matrix(1, K_all, length(th))
jsd_null <- phase_metrics_df(truth_anchor$phi, phi_null, tm)$JSD_bits

jsd_by_method <- lapply(reps, function(r) {
  a <- anchor_pair(r$phi, r$Sgrid, th)
  phase_metrics_df(truth_anchor$phi, a$phi, tm)$JSD_bits
})

## how far from uniform is the truth at each time point?
R_true <- vapply(seq_len(K_all),
                 function(k) Mod(phase_coeff(sim$phi_true[k, ], th, 1L)), numeric(1))
## NOTE: the resultant length of a von Mises law depends only on kappa, so with the two-level
## kappa_on / kappa_off design R_true takes just two distinct values and quantile breaks would not be
## unique. ntile() splits on ranks and is tie-safe. See the note after the figure.
tercile <- factor(dplyr::ntile(R_true, 3), levels = 1:3, labels = c("low", "medium", "high"))

## irreducible floor: the best any J-harmonic-truncated estimator could do
nvec_J <- make_nvec(Jc)
phi_trunc <- t(vapply(seq_len(K_all), function(k) {
  euclidean_simplex_project(
    synthesize_from_coeffs(discrete_fourier_coeffs(truth_anchor$phi[k, ], th, nvec_J), th, nvec_J),
    z = length(th))
}, numeric(length(th))))
jsd_floor <- phase_metrics_df(truth_anchor$phi, phi_trunc, tm)$JSD_bits

skill_long <- bind_rows(lapply(names(jsd_by_method), function(nm) {
  tibble(method = nm, time_index = seq_len(K_all), R_true = R_true, tercile = tercile,
         block = ifelse(sim$design == 1, "ON", "OFF"),
         jsd = jsd_by_method[[nm]], jsd_null = jsd_null,
         skill = 1 - jsd_by_method[[nm]] / pmax(jsd_null, 1e-12))
}))

skill_tbl <- skill_long |>
  group_by(method) |>
  summarise(mean_skill = mean(skill),
            skill_low_R = mean(skill[tercile == "low"]),
            skill_high_R = mean(skill[tercile == "high"]),
            win_rate = mean(jsd < jsd_null),
            win_rate_high_R = mean(jsd[tercile == "high"] < jsd_null[tercile == "high"]),
            .groups = "drop") |>
  arrange(desc(skill_high_R))
knitr::kable(skill_tbl, digits = 3,
             caption = "Skill against the uniform-phase null, overall and split by how far the true phase law is from uniform. Skill is 1 - JSD/JSD_null: 0 means no better than assuming uniform phases, 1 means perfect recovery. The low-R tercile is where the null is nearly correct by construction.")
Table 9.1: Skill against the uniform-phase null, overall and split by how far the true phase law is from uniform. Skill is 1 - JSD/JSD_null: 0 means no better than assuming uniform phases, 1 means perfect recovery. The low-R tercile is where the null is nearly correct by construction.
method mean_skill skill_low_R skill_high_R win_rate win_rate_high_R
KPT-FFT 0.338 0.466 0.219 0.980 1.00
KPT-GEM 0.155 0.143 0.132 0.913 0.98
Initializer 0.000 0.000 0.000 0.000 0.00
Laplace-ILT 0.000 0.000 0.000 0.000 0.00
Rank-transport 0.000 0.000 0.000 0.000 0.00
Stitch 0.000 0.000 0.000 0.000 0.00
Uniform-null 0.000 0.000 0.000 0.000 0.00 
cat(sprintf("Irreducible J = %d truncation floor: mean JSD = %.4f (the ceiling on achievable skill).\n",
            Jc, mean(jsd_floor)))
## Irreducible J = 16 truncation floor: mean JSD = 0.0000 (the ceiling on achievable skill).
## For each method and stratum, test whether the per-time divergence is smaller than the null's.
strata <- list(all = rep(TRUE, K_all),
               ON = sim$design == 1, OFF = sim$design == 0,
               high_R = tercile == "high", low_R = tercile == "low")

paired_vs_null <- bind_rows(lapply(names(jsd_by_method), function(nm) {
  bind_rows(lapply(names(strata), function(sn) {
    m <- strata[[sn]]
    d <- jsd_null[m] - jsd_by_method[[nm]][m]          # positive = method beats the null
    if (all(abs(d) < 1e-12)) {
      tibble(method = nm, stratum = sn, n_times = sum(m), mean_gain = 0, p_value = 1)
    } else {
      tibble(method = nm, stratum = sn, n_times = sum(m), mean_gain = mean(d),
             p_value = suppressWarnings(stats::wilcox.test(d, exact = FALSE))$p.value)
    }
  }))
})) |>
  filter(!method %in% c("Uniform-null", "Laplace-ILT", "Stitch", "Initializer")) |>
  mutate(p_BH = stats::p.adjust(p_value, method = "BH")) |>
  arrange(stratum, p_BH)

knitr::kable(paired_vs_null, digits = 4,
             caption = "Paired Wilcoxon tests of per-time-point JSD against the uniform-phase null, within strata, with Benjamini-Hochberg adjustment across all method-stratum combinations. Positive mean_gain favors the method.")
Table 9.2: Paired Wilcoxon tests of per-time-point JSD against the uniform-phase null, within strata, with Benjamini-Hochberg adjustment across all method-stratum combinations. Positive mean_gain favors the method.
method stratum n_times mean_gain p_value p_BH
KPT-FFT OFF 75 0.0347 0e+00 0.000
KPT-GEM OFF 75 0.0140 0e+00 0.000
Rank-transport OFF 75 0.0000 1e+00 1.000
KPT-GEM ON 75 0.0617 0e+00 0.000
KPT-FFT ON 75 0.1092 0e+00 0.000
Rank-transport ON 75 0.0000 1e+00 1.000
KPT-FFT all 150 0.0719 0e+00 0.000
KPT-GEM all 150 0.0379 0e+00 0.000
Rank-transport all 150 0.0000 1e+00 1.000
KPT-FFT high_R 50 0.1038 0e+00 0.000
KPT-GEM high_R 50 0.0626 0e+00 0.000
Rank-transport high_R 50 0.0000 1e+00 1.000
KPT-FFT low_R 50 0.0362 0e+00 0.000
KPT-GEM low_R 50 0.0111 7e-04 0.001
Rank-transport low_R 50 0.0000 1e+00 1.000 
skill_long |>
  filter(method %in% c("KPT-FFT", "KPT-GEM", "Rank-transport")) |>
  ggplot(aes(tercile, skill)) +
  geom_hline(yintercept = 0, linewidth = 0.4, colour = "grey40") +
  geom_boxplot(outlier.size = 0.8, fill = "grey92", linewidth = 0.4) +
  geom_jitter(aes(colour = block), width = 0.12, height = 0, size = 1.3, alpha = 0.8) +
  facet_wrap(~ method) +
  labs(title = "Skill against the uniform-phase null, by how concentrated the truth is",
       subtitle = "terciles of R(t) = |phi-hat_t(1)|; R = 0 is the uniform law, where no estimator can gain",
       x = "R(t) tercile", y = "skill vs null") +
  theme(legend.position = "bottom", legend.title = element_blank())

A caveat about the stratification. The mean resultant length of a von Mises law depends only on \(\kappa\), so with the two-level kappa_on / kappa_off design R_true takes exactly two values and the terciles nearly coincide with the ON/OFF blocks; the stratification then adds little beyond the block split. For a genuine gradient, give the simulator a continuous concentration by replacing the if (design[k] == 1) kappa_on else kappa_off argument in simulate_fmri() with kappa_off + (kappa_on - kappa_off) * bold[k], so that concentration tracks the haemodynamic response rather than the block indicator.

Reading the output. Report skill_high_R and the paired test within the high-\(R\) stratum as the headline, with the overall average alongside, and state the truncation floor so the reader knows the ceiling. A phase estimator can only be rewarded where there is phase structure to find.

fit_fft <- reps[["KPT-FFT"]]$extra
tibble(iteration = seq_along(fit_fft$loglik_trace), loglik = fit_fft$loglik_trace) |>
  ggplot(aes(iteration, loglik)) +
  geom_line(linewidth = 0.8) + geom_point(size = 1.2) +
  labs(title = sprintf("Safeguarded KPT: observed-data log-likelihood (monotone = %s)",
                       fit_fft$monotone),
       subtitle = sprintf("%d iterations, damping = %.2f; every accepted step is non-decreasing by construction",
                          fit_fft$iterations, fit_fft$parameters$damping),
       x = "iteration", y = "observed-data log-likelihood")

lp <- reps[["Laplace-ILT"]]$extra
tibble(time = tm, Original = lp$mu_hat, Recovered = lp$mu_roundtrip) |>
  pivot_longer(-time, names_to = "series", values_to = "value") |>
  ggplot(aes(time, value, colour = series, linetype = series)) +
  geom_line(linewidth = 0.7) +
  labs(title = "Laplace route: transform and numerical inverse of the ensemble mean",
       subtitle = sprintf("relative round-trip RMSE = %.3g (this is what the analytic route guarantees; it carries no replicate-phase information)",
                          lp$roundtrip_rel_rmse),
       x = "time", y = "signal") +
  theme(legend.position = "bottom", legend.title = element_blank())

k_slices <- round(stats::quantile(on_idx, c(0.25, 0.75), names = FALSE))
plot_df <- bind_rows(lapply(k_slices, function(k) {
  bind_rows(
    tibble(theta = th, value = truth_anchor$phi[k, ], series = "Truth", k = k),
    tibble(theta = th, value = anchor_pair(reps[["KPT-FFT"]]$phi, reps[["KPT-FFT"]]$Sgrid, th)$phi[k, ],
           series = "KPT-FFT", k = k),
    tibble(theta = th, value = rep(1, length(th)), series = "Uniform null", k = k)
  )
}))
ggplot(plot_df, aes(theta, value, colour = series, linetype = series)) +
  geom_line(linewidth = 0.8) + facet_wrap(~ paste("time index", k)) +
  scale_x_continuous(breaks = c(0, pi, 2*pi), labels = c("0", expression(pi), expression(2*pi))) +
  labs(title = "Anchored phase laws at two activation time points",
       x = expression(theta), y = expression(varphi[t](theta))) +
  theme(legend.position = "bottom", legend.title = element_blank())

10 Inference: calibration, power, and application

10.1 Does the uniformity test hold its nominal level?

k_ref <- on_idx[round(length(on_idx) / 2)]
val_v8 <- validate_rayleigh(sim_v8$S_true[k_ref, ], th, nvec_diag,
                            sigma = sim$sigma, N = CFG$N_rep, B = CFG$B_calib)
val_v6 <- validate_rayleigh(sim_v6$S_true[k_ref, ], th, nvec_diag,
                            sigma = sim$sigma, N = CFG$N_rep, B = CFG$B_calib)

size_tbl <- bind_rows(
  val_v8$table |> filter(scenario == "kappa=0") |> mutate(surface = "V8 (generic)"),
  val_v6$table |> filter(scenario == "kappa=0") |> mutate(surface = "V6 (reflection-degenerate)")
) |> select(surface, level, rejection_rate, mc_se)
knitr::kable(size_tbl, digits = 4,
             caption = "Empirical type I error of the rank-aware uniformity test under H0. Rejection rates should sit within about two Monte Carlo standard errors of the nominal level.")
Table 10.1: Empirical type I error of the rank-aware uniformity test under H0. Rejection rates should sit within about two Monte Carlo standard errors of the nominal level.
surface level rejection_rate mc_se
V8 (generic) 0.10 0.1040 0.0067
V8 (generic) 0.05 0.0455 0.0049
V8 (generic) 0.01 0.0095 0.0022
V6 (reflection-degenerate) 0.10 0.1170 0.0067
V6 (reflection-degenerate) 0.05 0.0680 0.0049
V6 (reflection-degenerate) 0.01 0.0150 0.0022 
rk <- rayleigh_uniformity(sim$Y[, k_ref], sim_v8$S_true[k_ref, ], sim$sigma, th, nvec_diag)
rk6 <- rayleigh_uniformity(sim_v6$Y[, k_ref], sim_v6$S_true[k_ref, ], sim$sigma, th, nvec_diag)
tibble(surface = c("V8 (generic)", "V6 (reflection-degenerate)"),
       effective_df = c(rk$df, rk6$df),
       eigenvalue_ratio = signif(c(rk$rank_ratio, rk6$rank_ratio), 3)) |>
  knitr::kable(caption = "Effective rank of the first-harmonic information matrix. The degenerate surface yields rank 1: referring the statistic to chi-square with 2 degrees of freedom, as an unguarded implementation would, rejects almost always under the null.")
Table 10.1: Effective rank of the first-harmonic information matrix. The degenerate surface yields rank 1: referring the statistic to chi-square with 2 degrees of freedom, as an unguarded implementation would, rejects almost always under the null.
surface effective_df eigenvalue_ratio
V8 (generic) 2 0.0565
V6 (reflection-degenerate) 1 0.0000

10.2 Power

val_v8$table |>
  filter(level == 0.05) |>
  mutate(kappa = as.numeric(sub("kappa=", "", scenario))) |>
  ggplot(aes(kappa, rejection_rate)) +
  geom_hline(yintercept = 0.05, linetype = 2, colour = "grey40") +
  geom_line(linewidth = 0.8) + geom_point(size = 1.6) +
  scale_y_continuous(limits = c(0, 1)) +
  labs(title = "Power of the uniformity test at the 5% level",
       subtitle = sprintf("N = %d replicates, sigma = %.2f; kappa = 0 is the null", CFG$N_rep, sim$sigma),
       x = expression(kappa~"(phase concentration)"), y = "rejection rate")

10.3 Asymptotic vs. bootstrap calibration, and testing across the design

p_asym <- vapply(seq_along(tm), function(k)
  rayleigh_uniformity(Ytr[, k], reps[["KPT-FFT"]]$Sgrid[k, ], reps[["KPT-FFT"]]$sigma, th, nvec_diag)$p_value,
  numeric(1))
p_bh <- stats::p.adjust(p_asym, method = "BH")

k_boot <- c(on_idx[1], off_idx[1])
boot_cmp <- tibble(
  time_index = k_boot,
  block = ifelse(sim$design[k_boot] == 1, "ON", "OFF"),
  p_asymptotic = signif(p_asym[k_boot], 3),
  p_bootstrap = signif(vapply(k_boot, function(k)
    rayleigh_bootstrap_p(Ytr[, k], reps[["KPT-FFT"]]$Sgrid[k, ], reps[["KPT-FFT"]]$sigma,
                         th, nvec_diag, B = 99L), numeric(1)), 3)
)
knitr::kable(boot_cmp, caption = "Asymptotic versus parametric-bootstrap p-values at one ON and one OFF time point.")
Table 10.2: Asymptotic versus parametric-bootstrap p-values at one ON and one OFF time point.
time_index block p_asymptotic p_bootstrap
16 ON 0e+00 0.01
1 OFF 2e-07 0.01 
tibble(block = c("ON", "OFF"),
       n_times = c(length(on_idx), length(off_idx)),
       prop_rejected_BH_05 = c(mean(p_bh[on_idx] < 0.05), mean(p_bh[off_idx] < 0.05))) |>
  knitr::kable(digits = 3,
               caption = "Proportion of time points declared non-uniform after Benjamini-Hochberg control across all K times. The design predicts many rejections during ON blocks (kappa = 6) and few during OFF blocks (kappa = 1).")
Table 10.2: Proportion of time points declared non-uniform after Benjamini-Hochberg control across all K times. The design predicts many rejections during ON blocks (kappa = 6) and few during OFF blocks (kappa = 1).
block n_times prop_rejected_BH_05
ON 75 0.880
OFF 75 0.947 
perm <- permutation_phase_test(anchor_rows(reps[["KPT-FFT"]]$phi, th)$phi,
                               ifelse(sim$design == 1, "ON", "OFF"), th)
cat(sprintf("Rotation-invariant Sobolev discrepancy between ON and OFF phase laws: D = %.4g, permutation p = %.4f\n",
            perm$statistic, perm$p_value))
## Rotation-invariant Sobolev discrepancy between ON and OFF phase laws: D = 0.001866, permutation p = 0.0020

10.4 Simultaneous bootstrap bands

t_band <- system.time({
  bands <- bootstrap_phase_bands(Ytr, tm, th, builder = rep_transport,
                                 B = CFG$B_boot, level = 0.90, J = Jc)
})
kb <- on_idx[round(length(on_idx) / 2)]
tibble(theta = th, centre = bands$centre[kb, ],
       lower = pmax(bands$lower[kb, ], 0), upper = bands$upper[kb, ],
       truth = truth_anchor$phi[kb, ]) |>
  ggplot(aes(theta)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80", alpha = 0.7) +
  geom_line(aes(y = centre), linewidth = 0.8) +
  geom_line(aes(y = truth), linetype = 2, linewidth = 0.7, colour = "firebrick") +
  scale_x_continuous(breaks = c(0, pi, 2*pi), labels = c("0", expression(pi), expression(2*pi))) +
  labs(title = "Simultaneous 90% bootstrap band for the anchored phase law",
       subtitle = sprintf("solid: bootstrap median; dashed: truth; sup-t calibration, B = %d on %d core(s)",
                          bands$B, bands$cores),
       x = expression(theta), y = expression(varphi[t](theta)))

cat(sprintf("elapsed %.1f s on %d core(s)\n", t_band[["elapsed"]], PAR$cores))
## elapsed 2390.2 s on 8 core(s)

11 Double-slit application

ds <- simulate_double_slit()
ds_reps <- list(
  `Uniform-null`   = rep_uniform(ds$Y, ds$time, ds$theta, J = Jc),
  `Rank-transport` = rep_transport(ds$Y, ds$time, ds$theta, J = Jc),
  `KPT-FFT`        = rep_kpt(ds$Y, ds$time, ds$theta, variant = "FFT", J = Jc,
                             max_iter = CFG$max_iter, sigma = ds$sigma)
)
ds_truth <- anchor_pair(ds$phi_true, ds$S_true, ds$theta)
bind_rows(lapply(names(ds_reps), function(nm) {
  a <- anchor_pair(ds_reps[[nm]]$phi, ds_reps[[nm]]$Sgrid, ds$theta)
  pm <- phase_metrics_df(ds_truth$phi, a$phi, ds$time)
  sc <- score_representation(ds_reps[[nm]], ds$Y_test)
  tibble(method = nm, mean_JSD_bits = mean(pm$JSD_bits),
         surface_rel_L2 = surface_L2_cone(ds_truth$S, a$S, ds$time)$rel_L2,
         logscore = sc$logscore, CRPS = sc$CRPS, coverage = sc$coverage)
})) |>
  knitr::kable(digits = 4,
               caption = "Double-slit simulation: the phase law is a two-component von Mises mixture, so it is genuinely far from uniform and a phase estimator has room to improve on the null.")
Table 11.1: Double-slit simulation: the phase law is a two-component von Mises mixture, so it is genuinely far from uniform and a phase estimator has room to improve on the null.
method mean_JSD_bits surface_rel_L2 logscore CRPS coverage
Uniform-null 0.2058 0.2971 -0.0626 0.1715 0.8791
Rank-transport 0.2058 0.2991 -0.0620 0.1715 0.8827
KPT-FFT 0.1322 0.1886 -0.0677 0.1730 0.9166

12 Self-tests

These assertions run at render time and stop the build on failure. They encode the invariants that the V6 code silently violated.

tests_run <- 0L
check <- function(cond, what) {
  tests_run <<- tests_run + 1L
  if (!isTRUE(cond)) stop("SELF-TEST FAILED: ", what, call. = FALSE)
  invisible(TRUE)
}
tt_theta <- theta_grid(64L)

## 1. simplex projection: non-negative, mean one, idempotent
v <- stats::rnorm(64, 1, 1); pv <- euclidean_simplex_project(v, z = 64)
check(all(pv >= 0), "simplex projection non-negativity")
check(abs(mean(pv) - 1) < 1e-10, "simplex projection mean equals one")
check(max(abs(pv - euclidean_simplex_project(pv, z = 64))) < 1e-10, "simplex projection idempotent")

## 2. circular shift returns x(theta - delta)
g <- vm_density_row(tt_theta, 0.7, 4)
check(max(abs(shift_periodic_row(g, 0.5) - vm_density_row(tt_theta, 1.2, 4))) < 1e-8,
      "shift_periodic_row implements x(theta - delta)")

## 3. anchoring sends the mode to zero (the V6 sign error sent it to 2*mu)
for (mu in c(0.4, 2.0, 4.5)) {
  a <- anchor_rows(matrix(vm_density_row(tt_theta, mu, 4), nrow = 1L), tt_theta)
  mode_at <- tt_theta[which.max(a$phi[1, ])]
  check(min(abs(mode_at - c(0, 2*pi))) < 2 * (2*pi / 64),
        sprintf("anchor maps mode %.2f to 0", mu))
}
## and anchoring is idempotent
a1 <- anchor_rows(matrix(vm_density_row(tt_theta, 2.0, 4), nrow = 1L), tt_theta)
a2 <- anchor_rows(a1$phi, tt_theta)
check(max(abs(a1$phi - a2$phi)) < 1e-6, "anchoring is idempotent")

## 4. Fourier round trip on the phase grid
nv <- make_nvec(8L)
x <- 1 + 0.4 * cos(tt_theta) + 0.2 * sin(3 * tt_theta)
check(max(abs(synthesize_from_coeffs(discrete_fourier_coeffs(x, tt_theta, nv), tt_theta, nv) - x)) < 1e-8,
      "Fourier analysis/synthesis round trip")

## 5. reflection diagnostic flags an even surface and clears a generic one
S_even <- matrix(1.6 + cos(tt_theta) + 0.25 * cos(2 * tt_theta), nrow = 1L)
S_gen  <- matrix(1.6 + cos(tt_theta) + 0.35 * sin(2 * tt_theta), nrow = 1L)
check(reflection_diagnostic(S_even, tt_theta)$reflection_degenerate[1], "reflection diagnostic detects even surface")
check(!reflection_diagnostic(S_gen, tt_theta)$reflection_degenerate[1], "reflection diagnostic clears generic surface")

## 6. aliasing diagnostic recovers d_t = 2 for a surface built from even harmonics only
S_alias <- matrix(2 + cos(2 * tt_theta) + 0.3 * cos(4 * tt_theta), nrow = 1L)
check(aliasing_diagnostic(S_alias, tt_theta, nv)$d_t[1] == 2L, "aliasing diagnostic recovers d_t = 2")

## 7. safeguarded KPT is monotone in the observed-data log-likelihood
sim_small <- simulate_fmri(N = 60L, K = 12L, L_theta = 64L, n_split = 1L, seed = 5)
fit_small <- kpt_fit(sim_small$Y, sim_small$time, sim_small$theta, J = 6L,
                     max_iter = 8L, sigma = sim_small$sigma)
check(fit_small$monotone, "safeguarded KPT log-likelihood is non-decreasing")

## 8. predictive density integrates to one and the CDF is consistent with it
rp <- rep_uniform(sim_small$Y, sim_small$time, sim_small$theta, J = 6L)
yg <- seq(min(sim_small$Y) - 4 * rp$sigma, max(sim_small$Y) + 4 * rp$sigma, length.out = 2001)
check(abs(trapz_vec(yg, predictive_density(rp, 1L, yg)) - 1) < 1e-3, "predictive density integrates to one")
check(max(abs(predictive_cdf(rp, 1L, c(-1e6, 1e6)) - c(0, 1))) < 1e-8, "predictive CDF limits")

## 9. ILT inverts transforms whose answers are known in closed form
check(abs(ilt_valsa(function(z) 1 / (z + 1), 1.0) - exp(-1)) < 1e-6, "ILT recovers exp(-t)")
check(abs(ilt_valsa(function(z) 1 / z^2, 2.0) - 2) < 1e-5, "ILT recovers t")

## 10. the uniformity test is exactly degenerate on an even surface
r_even <- rayleigh_uniformity(as.vector(S_even)[sample.int(64L, 80L, TRUE)] + stats::rnorm(80, 0, 0.2),
                              as.vector(S_even), 0.2, tt_theta, nv)
check(r_even$df == 1L, "uniformity test detects rank-one information on an even surface")

## 11. the parallel helper agrees with plain lapply and reports a usable backend
check(identical(par_lapply(1:4, function(i) i^2, backend = list(type = "serial", cores = 1L)),
                lapply(1:4, function(i) i^2)), "par_lapply serial path matches lapply")
check(PAR$cores >= 1L && PAR$type %in% c("serial", "fork", "psock"), "parallel backend is well formed")

cat(sprintf("All %d self-tests passed.\n", tests_run))
## All 20 self-tests passed.

13 Reproducibility

tibble(
  Setting = c("Phase grid L_theta", "Harmonic truncation J", "Time points K",
              "Replicates per split N", "Bootstrap B", "Calibration replications",
              "KPT max iterations", "KPT damping", "Safeguard", "Parallel backend",
              "Master seed", "FAST mode"),
  Value = c(CFG$L_theta, CFG$J, CFG$K_fmri, CFG$N_rep, CFG$B_boot, CFG$B_calib,
            CFG$max_iter, 0.5, "log-likelihood backtracking",
            sprintf("%s x %d", PAR$type, PAR$cores), 20260912, FAST)
) |> knitr::kable(caption = "Configuration actually used in this render. Set FAST <- FALSE in the setup chunk for paper-scale runs.")
Table 13.1: Configuration actually used in this render. Set FAST <- FALSE in the setup chunk for paper-scale runs.
Setting Value
Phase grid L_theta 256
Harmonic truncation J 16
Time points K 150
Replicates per split N 200
Bootstrap B 200
Calibration replications 2000
KPT max iterations 40
KPT damping 0.5
Safeguard log-likelihood backtracking
Parallel backend psock x 8
Master seed 20260912
FAST mode FALSE 
sessionInfo()
## R version 4.3.3 (2024-02-29 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] tibble_3.2.1  knitr_1.51    purrr_1.0.2   tidyr_1.3.1   dplyr_1.1.4  
## [6] ggplot2_4.0.1
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5        cli_3.6.3          rlang_1.1.5        xfun_0.52         
##  [5] otel_0.2.0         generics_0.1.3     S7_0.2.1           jsonlite_1.8.9    
##  [9] labeling_0.4.3     glue_1.8.0         htmltools_0.5.8.1  sass_0.4.9        
## [13] scales_1.4.0       rmarkdown_2.31     grid_4.3.3         evaluate_1.0.3    
## [17] jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.10        lifecycle_1.0.5   
## [21] bookdown_0.40      compiler_4.3.3     RColorBrewer_1.1-3 pkgconfig_2.0.3   
## [25] rstudioapi_0.18.0  farver_2.1.2       digest_0.6.37      R6_2.6.1          
## [29] tidyselect_1.2.1   parallel_4.3.3     pillar_1.10.1      magrittr_2.0.3    
## [33] bslib_0.9.0        withr_3.0.2        tools_4.3.3        gtable_0.3.6      
## [37] cachem_1.1.0

14 Interpretation and limitations

What the comparison does and does not establish. The predictive tables rank representations on held-out replicates from the same generative model. That is the right criterion for choosing a representation for a given data set, and it is available on real data, but it does not certify that any of these strategies recovers a physically meaningful phase. Recovery tables do speak to that, and they are only available in simulation.

KPT vs the null. Safeguarded KPT improves both surface recovery and phase recovery over the uniform-phase null, but the margin is strongly heterogeneous across time and is easy to lose. Two effects matter. First, the V6 anchoring sign error by itself destroys most of the apparent advantage: on one fixed fit, measuring with the V6 anchor gives a skill of \(+0.02\) against the null, while measuring the same fit with the corrected anchor gives \(+0.16\). Second, half the time points are rest blocks with \(\kappa=1\), where the truth is nearly uniform and the null is nearly correct by construction, so averaging over all times dilutes a real effect. Section 9 reports the stratified comparison that keeps both effects visible.

The diagnostics are prerequisites, not decorations. With a reflection-degenerate surface, no estimator can identify \(\varphi\) beyond its symmetric part, and first-harmonic inference has one degree of freedom rather than two. With \(\min_\theta|\mathcal S|\approx0\), the deconvolution is locally uninformative no matter how many replicates are collected. Both conditions are properties of the surface and should be reported alongside any reconstruction.

Where the evidence for each correction lives. The claims in the changelog are not self-certifying: each was established by a separate experiment in KPT_audit_V8.ipynb, and the self-tests in Section 12 guard the fixed behaviour at render time. If a future edit breaks one of them, the render stops rather than silently reporting a wrong number.

Remaining gaps. The noise model is homoscedastic Gaussian and the phases are independent across time; neither is realistic for fMRI, where noise is autocorrelated and phases plausibly evolve smoothly. The registration routine of Section 2.5 chooses rotations greedily and could be replaced by a global dynamic-programming alignment. Finally, the asymptotic uniformity test is calibrated here only for the specific surfaces studied; the bootstrap route is the safer default whenever \(\min_\theta|\mathcal S|\) is small.