SOCR ≫ TCIU Website ≫ TCIU GitHub ≫

This TCIU Appendix offers a self-contained implementation of the kime phase tompgraphy algorithms, KPT-GEM (expectation maximization) and KPT-FFT (FFT alternating spectral) for estimating the phase distribution and reconstructing the kime-surface representations from observed repeated measurement data from longitudinal processes. The KPT algorihtm is validated using two simulations; double-slit physics experiment and fMRI ON vs. OFF neuroscience experiment. We report quyantitative performance metrics (summary tables) and provide qualitative results as surface/heatmap overlays contrasting for the true vs. estimated phase distributions and kime-surfaces.

The key steps of the KPT algorithm include E-step, Wiener-Sobolev filtering, simplex projection, ans reality constraints, subject to normalization, operator conventions.

V.6 revises and expands the previous V.5 of the TCIU KPT appendix 3. It is aligned with the rewritten manuscript new_SKA_KPT_V6.tex and is organized to mirror the paper’s structure: introduction, mathematical foundations, KPT algorithms, statistical inference, applications, and conclusions. The notebook is fully self-contained, uses the time-specific latent phase model throughout, adopts the normalization convention based on the Haar probability measure \(d\theta/(2\pi)\), and removes the inconsistencies and exploratory sections present in earlier drafts.

The core generative model used everywhere in this protocol is

\[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),\] where each phase law \(\varphi_{t_k}\) is a density on the circle with respect to \(d\theta/(2\pi)\). On a discrete phase grid \(\{\theta_\ell\}_{\ell=0}^{L_\theta-1}\), the corresponding normalization is

\[\frac{1}{L_\theta}\sum_{\ell=0}^{L_\theta-1} \varphi_{t_k}(\theta_\ell) = 1.\] Accordingly, the uniform phase law is represented by a vector of ones.

1 Introduction

This appendix implements the computational procedures described in the revised paper. The emphasis here is practical rather than theorem-proving. In particular, the notebook implements the projected practical algorithms used in the numerical experiments:

  • KPT-GEM (practical phase-update variant): Wiener–Sobolev phase updates with a fixed initialized surface.
  • KPT-FFT (practical alternating variant): projected phase updates alternating with surface updates.

The revised paper distinguishes these computational routines from the exact constrained EM reference algorithm, which is the appropriate object for strict monotonicity claims. That exact reference algorithm is not the default numerical engine in this notebook. This notebook therefore reports empirical convergence diagnostics for the practical projected routines and does not interpret them as globally monotone EM algorithms.

2 Mathematical Foundations

2.1 Model, normalization, and alignment conventions

Throughout the notebook, phase densities are normalized with respect to the Haar probability measure \(d\theta/(2\pi)\). On the discrete grid this means row mean = 1. All simulation, estimation, alignment, prediction, and evaluation steps below use the same convention.

Because the latent phase origin is a gauge degree of freedom, comparisons between true and estimated phase laws are performed after an explicit anchoring step. In the numerical studies below, we use the first trigonometric moment anchor whenever the first harmonic is non-negligible; if the first moment is nearly zero, the row is left unchanged.

2.2 Simulation and fitting parameters copied from the paper

simulation_parameters <- tibble::tribble(
  ~Setting, ~Double_slit, ~fMRI,
  "Observation times",
  "K = 50, equally spaced on [0, 1]",
  "K = 150, sampled every 2 seconds on [0, 298]",
  "Phase grid",
  "L_theta = 256",
  "L_theta = 256",
  "Replicates per time point",
  "N = 100",
  "P x R = 10 x 20 = 200",
  "Noise standard deviation",
  "sigma = 0.15",
  "sigma = 0.20",
  "Latent phase law",
  "Two-component von Mises mixture with kappa_phase = 4",
  "von Mises with kappa_on = 6 and kappa_off = 1",
  "Surface generator",
  "Baseline + Gaussian envelope + first/second harmonics with kappa_wave = 12",
  "Baseline drift + HRF-convolved block design + first/second harmonics",
  "KPT-GEM tuning",
  "J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = FALSE, max_iter = 30, tol = 1e-3",
  "J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = FALSE, max_iter = 25, tol = 8e-4",
  "KPT-FFT tuning",
  "J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = TRUE, max_iter = 30, tol = 1e-3",
  "J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = TRUE, max_iter = 25, tol = 8e-4",
  "Evaluation alignment",
  "First-moment anchor; the same rotation is applied to the surface",
  "First-moment anchor; the same rotation is applied to the surface"
)

knitr::kable(
  simulation_parameters,
  caption = "Simulation and fitting parameters used throughout this notebook. These match the revised paper."
)
Simulation and fitting parameters used throughout this notebook. These match the revised paper.
Setting Double_slit fMRI
Observation times K = 50, equally spaced on [0, 1] K = 150, sampled every 2 seconds on [0, 298]
Phase grid L_theta = 256 L_theta = 256
Replicates per time point N = 100 P x R = 10 x 20 = 200
Noise standard deviation sigma = 0.15 sigma = 0.20
Latent phase law Two-component von Mises mixture with kappa_phase = 4 von Mises with kappa_on = 6 and kappa_off = 1
Surface generator Baseline + Gaussian envelope + first/second harmonics with kappa_wave = 12 Baseline drift + HRF-convolved block design + first/second harmonics
KPT-GEM tuning J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = FALSE, max_iter = 30, tol = 1e-3 J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = FALSE, max_iter = 25, tol = 8e-4
KPT-FFT tuning J = 16, lambda_phi = lambda_F = 5e-3, p_phi = p_F = 1, update_surface = TRUE, max_iter = 30, tol = 1e-3 J = 12, lambda_phi = lambda_F = 1e-2, p_phi = p_F = 1, update_surface = TRUE, max_iter = 25, tol = 8e-4
Evaluation alignment First-moment anchor; the same rotation is applied to the surface First-moment anchor; the same rotation is applied to the surface

3 Kime-Phase Tomography (KPT)

3.1 Core utilities

# ---- Grid helpers -------------------------------------------------------------
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))
}

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

# ---- 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))
  theta <- cssv[rho] / rho
  pmax(v - theta, 0)
}

row_project_simplex <- function(P) {
  P <- as.matrix(P)
  out <- t(apply(P, 1L, euclidean_simplex_project, z = ncol(P)))
  rownames(out) <- rownames(P)
  colnames(out) <- colnames(P)
  out
}

# ---- Fourier utilities --------------------------------------------------------
discrete_fourier_coeffs <- function(x_theta, theta, nvec) {
  Eneg <- exp(-1i * outer(theta, nvec))
  as.vector(drop(t(x_theta) %*% Eneg) / length(theta))
}

synthesize_from_coeffs <- function(coeffs_n, theta, nvec) {
  Epos <- exp(1i * outer(nvec, theta))
  as.vector(Re(drop(coeffs_n %*% Epos)))
}

enforce_reality_coeffs <- function(coeffs_n, nvec) {
  coeffs_n <- as.complex(coeffs_n)
  coeffs_n[nvec == 0] <- Re(coeffs_n[nvec == 0])
  J <- max(abs(nvec))
  for (j in seq_len(J)) {
    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)
  for (i in seq_along(nvec)) {
    m <- (nvec[i] + N) %% N
    out[m + 1L] <- x_n[i]
  }
  out
}

dft_order_to_coeffs <- function(x_m, nvec) {
  N <- length(nvec)
  out <- complex(length = N)
  for (i in seq_along(nvec)) {
    m <- (nvec[i] + N) %% N
    out[i] <- x_m[m + 1L]
  }
  out
}

coeffs_to_Xomega_fft <- function(x_n, nvec) {
  fft(coeffs_to_dft_order(x_n, nvec), inverse = FALSE)
}

Xomega_to_coeffs_fft <- function(X_w, nvec) {
  N <- length(nvec)
  dft_order_to_coeffs(fft(X_w, inverse = TRUE) / N, nvec)
}

sobolev_weights <- function(nvec, p = 1L) {
  abs(nvec)^(2 * p)
}

# ---- Periodic interpolation and circular shifts -------------------------------
periodic_interp <- function(theta, y, theta_new) {
  stopifnot(length(theta) == length(y))
  theta_ext <- c(theta, 2*pi)
  y_ext <- c(y, y[1L])
  xout <- theta_new %% (2*pi)
  approx(theta_ext, y_ext, xout = xout, rule = 2L, ties = mean)$y
}

shift_periodic_row <- function(x, delta) {
  L <- length(x)
  n <- freq_index(L)
  X <- fft(x)
  Re(fft(X * exp(-1i * n * delta), inverse = TRUE) / L)
}

shift_periodic_matrix <- function(X, delta_vec) {
  stopifnot(nrow(X) == length(delta_vec))
  out <- matrix(0, nrow = nrow(X), ncol = ncol(X))
  for (k in seq_len(nrow(X))) out[k, ] <- shift_periodic_row(X[k, ], delta_vec[k])
  out
}

# ---- Alignment helpers --------------------------------------------------------
first_trig_moment <- function(phi_mat, theta, harmonic = 1L) {
  apply(phi_mat, 1L, function(row) mean(row * exp(1i * harmonic * theta)))
}

anchor_by_first_moment <- function(phi_mat, theta, harmonic = 1L, tol = 1e-8) {
  m <- first_trig_moment(phi_mat, theta, harmonic = harmonic)
  delta <- ifelse(Mod(m) > tol, Arg(m) / harmonic, 0)
  phi_rot <- shift_periodic_matrix(phi_mat, delta)
  phi_rot <- row_project_simplex(phi_rot)
  list(phi = phi_rot, delta = delta, moment = m)
}

apply_anchor_to_surface <- function(S_mat, delta_vec) {
  shift_periodic_matrix(S_mat, delta_vec)
}

# ---- Surface initialization ---------------------------------------------------
initialize_surface <- function(Y, time, theta, method = c("empirical", "constant")) {
  method <- match.arg(method)
  K <- ncol(Y)
  L <- length(theta)
  if (method == "constant") return(matrix(mean(Y), nrow = K, ncol = L))

  mu_k <- colMeans(Y)
  sd_k <- apply(Y, 2L, stats::sd)

  mu_s <- if (K >= 5L) stats::smooth.spline(time, mu_k, spar = 0.6)$y else mu_k
  sd_s <- if (K >= 5L) stats::smooth.spline(time, sd_k, spar = 0.6)$y else sd_k

  S <- matrix(0, nrow = K, ncol = L)
  for (k in seq_len(K)) {
    S[k, ] <- mu_s[k] + sd_s[k] * cos(theta) + 0.25 * sd_s[k] * cos(2 * theta)
  }
  S
}

# ---- Criterion tracking (observed-data log-likelihood) ------------------------
observed_loglik_gaussian <- function(Y, Sgrid, phi_grid, sigma) {
  N <- nrow(Y)
  K <- ncol(Y)
  ll <- 0
  for (k in seq_len(K)) {
    weights <- pmax(phi_grid[k, ], 0)
    weights <- weights / sum(weights)
    dens_mat <- outer(Y[, k], Sgrid[k, ], function(y, mu) stats::dnorm(y, mean = mu, sd = sigma))
    mix_dens <- drop(dens_mat %*% weights)
    ll <- ll + sum(log(pmax(mix_dens, 1e-300)))
  }
  ll
}

# ---- Metrics -----------------------------------------------------------------
normalize_rows_prob <- function(P) {
  P <- pmax(P, 1e-15)
  sweep(P, 1L, rowSums(P), "/")
}

phase_metrics_df <- function(phi_true, phi_est, time, base2 = TRUE) {
  P <- normalize_rows_prob(phi_true)
  Q <- normalize_rows_prob(phi_est)
  log_fun <- if (base2) function(x) log2(pmax(x, 1e-15)) else function(x) log(pmax(x, 1e-15))
  K <- nrow(P)
  out <- tibble(
    time = time,
    KL_true_est = numeric(K),
    KL_est_true = numeric(K),
    JSD_bits = numeric(K),
    Hellinger = numeric(K),
    TV = numeric(K)
  )
  for (k in seq_len(K)) {
    M <- 0.5 * (P[k, ] + Q[k, ])
    out$KL_true_est[k] <- sum(P[k, ] * log_fun(P[k, ] / Q[k, ]))
    out$KL_est_true[k] <- sum(Q[k, ] * log_fun(Q[k, ] / P[k, ]))
    out$JSD_bits[k] <- 0.5 * sum(P[k, ] * log_fun(P[k, ] / M)) +
      0.5 * sum(Q[k, ] * log_fun(Q[k, ] / M))
    out$Hellinger[k] <- sqrt(0.5 * sum((sqrt(P[k, ]) - sqrt(Q[k, ]))^2))
    out$TV[k] <- 0.5 * sum(abs(P[k, ] - Q[k, ]))
  }
  out
}

surface_L2_cone <- function(S_true, S_est, time) {
  stopifnot(all(dim(S_true) == dim(S_est)))
  err_theta <- rowMeans((S_true - S_est)^2)
  true_theta <- rowMeans(S_true^2)
  integrand_err <- time * err_theta
  integrand_true <- time * true_theta
  abs_sq <- trapz_vec(time, integrand_err)
  true_sq <- trapz_vec(time, integrand_true)
  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_theta)
  )
}

# ---- Predictive summaries -----------------------------------------------------
kpt_predictive_bands <- function(fit, B = 4000L,
                                 probs = c(0.05, 0.25, 0.5, 0.75, 0.95),
                                 seed = 123) {
  set.seed(seed)
  theta <- fit$theta
  Sgrid <- fit$anchored$Sgrid
  phi_grid <- fit$anchored$varphi_theta
  sigma <- fit$parameters$sigma
  time <- fit$time

  out <- vector("list", length(time))
  for (k in seq_along(time)) {
    w <- normalize_rows_prob(matrix(phi_grid[k, ], nrow = 1L))[1L, ]
    idx <- sample.int(length(theta), size = B, replace = TRUE, prob = w)
    mu <- Sgrid[k, idx]
    ysim <- mu + rnorm(B, mean = 0, sd = sigma)
    qs <- stats::quantile(ysim, probs = probs, names = FALSE)
    out[[k]] <- c(time = time[k], mean = mean(ysim), setNames(qs, paste0("q", probs)))
  }
  as.data.frame(do.call(rbind, out))
}

3.2 Practical KPT algorithms

# The practical estimator implemented here mirrors the revised paper:
# - phase densities live on the discrete simplex with row mean = 1
# - KPT-GEM uses projected phase updates with a fixed initialized surface
# - KPT-FFT alternates projected phase updates and surface updates
# - convergence is tracked using real differences, not overwritten objects
# - any optional log-likelihood trace is explicitly labeled as observed log-likelihood

kpt_fit <- function(Y, time = NULL,
                    J = 12, L_theta = 256,
                    lambda_phi = 1e-2, p_phi = 1,
                    lambda_F = 1e-2, p_F = 1,
                    algorithm = c("GEM", "FFT"),
                    update_surface = NULL,
                    max_iter = 30, tol = 1e-3,
                    sigma = NULL,
                    init_surface = c("empirical", "constant"),
                    track_loglik = FALSE,
                    verbose = TRUE) {
  algorithm <- match.arg(algorithm)
  init_surface <- match.arg(init_surface)
  if (is.null(update_surface)) update_surface <- identical(algorithm, "FFT")

  I <- nrow(Y)
  K <- ncol(Y)
  if (is.null(time)) time <- seq_len(K)

  theta <- theta_grid(L_theta)
  nvec <- make_nvec(J)
  Nn <- length(nvec)
  omega <- 2*pi * (0:(Nn - 1L)) / Nn
  Lambda_phi <- (2 * sin(omega / 2))^(2 * p_phi)
  Lambda_F <- (2 * sin(omega / 2))^(2 * p_F)

  if (is.null(sigma)) {
    sigma <- 0.5 * stats::mad(as.vector(Y), center = 0) + 1e-6
  }

  phi_grid <- matrix(1, nrow = K, ncol = L_theta)
  Sgrid <- initialize_surface(Y, time, theta, method = init_surface)

  f_n <- matrix(0 + 0i, nrow = K, ncol = Nn)
  for (k in seq_len(K)) {
    f_n[k, ] <- enforce_reality_coeffs(discrete_fourier_coeffs(Sgrid[k, ], theta, nvec), nvec)
  }

  delta_phi <- delta_f <- numeric(max_iter)
  obs_loglik <- rep(NA_real_, max_iter)

  for (it in seq_len(max_iter)) {
    phi_old <- phi_grid
    f_old <- f_n

    # E-step: posterior weights and mixed moments
    m_hat <- matrix(0 + 0i, nrow = K, ncol = Nn)
    for (k in seq_len(K)) {
      prior_k <- normalize_rows_prob(matrix(phi_grid[k, ], nrow = 1L))[1L, ]
      loglik_mat <- outer(Y[, k], Sgrid[k, ], function(y, mu) -(y - mu)^2 / (2 * sigma^2))
      logw <- sweep(loglik_mat, 2L, log(pmax(prior_k, 1e-15)), "+")
      W <- softmax_rows(logw)
      Eneg <- exp(-1i * outer(theta, nvec))
      xi <- W %*% Eneg
      m_hat[k, ] <- colMeans(Y[, k] * xi)
    }

    # Phase update: regularized deconvolution, then exact Euclidean projection
    phi_hat_n <- matrix(0 + 0i, nrow = K, ncol = Nn)
    phi_grid_new <- matrix(0, nrow = K, ncol = L_theta)
    for (k in seq_len(K)) {
      F_w <- coeffs_to_Xomega_fft(f_n[k, ], nvec)
      M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
      denom <- Mod(F_w)^2 + lambda_phi * Lambda_phi + 1e-10
      Phi_w <- Conj(F_w) * M_w / denom
      phi_hat_n[k, ] <- Xomega_to_coeffs_fft(Phi_w, nvec)
      phi_grid_new[k, ] <- synthesize_from_coeffs(phi_hat_n[k, ], theta, nvec)
    }
    phi_grid_new <- row_project_simplex(phi_grid_new)

    # Recompute phase coefficients from the projected densities
    phi_n_projected <- matrix(0 + 0i, nrow = K, ncol = Nn)
    for (k in seq_len(K)) {
      phi_n_projected[k, ] <- discrete_fourier_coeffs(phi_grid_new[k, ], theta, nvec)
    }

    # Optional surface update
    if (isTRUE(update_surface)) {
      f_n_new <- matrix(0 + 0i, nrow = K, ncol = Nn)
      Sgrid_new <- matrix(0, nrow = K, ncol = L_theta)
      for (k in seq_len(K)) {
        Phi_w <- coeffs_to_Xomega_fft(phi_n_projected[k, ], nvec)
        M_w <- coeffs_to_Xomega_fft(m_hat[k, ], nvec)
        denom <- Mod(Phi_w)^2 + lambda_F * Lambda_F + 1e-10
        F_w <- Conj(Phi_w) * M_w / denom
        f_n_new[k, ] <- enforce_reality_coeffs(Xomega_to_coeffs_fft(F_w, nvec), nvec)
        Sgrid_new[k, ] <- synthesize_from_coeffs(f_n_new[k, ], theta, nvec)
      }
      f_n <- f_n_new
      Sgrid <- Sgrid_new
    }

    phi_grid <- phi_grid_new

    delta_phi[it] <- sqrt(mean((phi_grid - phi_old)^2))
    delta_f[it] <- sqrt(mean(Mod(f_n - f_old)^2))

    if (isTRUE(track_loglik)) {
      obs_loglik[it] <- observed_loglik_gaussian(Y, Sgrid, phi_grid, sigma)
    }

    if (isTRUE(verbose)) {
      if (isTRUE(update_surface)) {
        message(sprintf("[%s] iter %02d: dphi = %.3e, df = %.3e", algorithm, it, delta_phi[it], delta_f[it]))
      } else {
        message(sprintf("[%s] iter %02d: dphi = %.3e", algorithm, it, delta_phi[it]))
      }
    }

    if (max(delta_phi[it], delta_f[it]) < tol) break
  }

  fit <- list(
    varphi_theta = phi_grid,
    Sgrid = Sgrid,
    phi_n = phi_n_projected,
    f_n = f_n,
    theta = theta,
    time = time,
    convergence = list(
      iterations = it,
      delta_phi = delta_phi[seq_len(it)],
      delta_f = delta_f[seq_len(it)],
      observed_loglik = obs_loglik[seq_len(it)]
    ),
    parameters = list(
      algorithm = algorithm,
      J = J,
      L_theta = L_theta,
      lambda_phi = lambda_phi,
      p_phi = p_phi,
      lambda_F = lambda_F,
      p_F = p_F,
      sigma = sigma,
      update_surface = update_surface,
      max_iter = max_iter,
      tol = tol,
      normalization = "row mean = 1"
    )
  )

  anchored <- anchor_by_first_moment(fit$varphi_theta, fit$theta)
  fit$anchored <- list(
    varphi_theta = anchored$phi,
    Sgrid = apply_anchor_to_surface(fit$Sgrid, anchored$delta),
    delta = anchored$delta
  )
  fit
}

fit_kpt_gem <- function(Y, time, ...) {
  kpt_fit(
    Y = Y,
    time = time,
    algorithm = "GEM",
    update_surface = FALSE,
    ...
  )
}

fit_kpt_fft <- function(Y, time, ...) {
  kpt_fit(
    Y = Y,
    time = time,
    algorithm = "FFT",
    update_surface = TRUE,
    ...
  )
}

3.3 Plotting helper functions

plot_phase_slice <- function(theta, truth, est_list, title = "", subtitle = "") {
  df <- tibble(theta = theta, Truth = truth)
  for (nm in names(est_list)) df[[nm]] <- est_list[[nm]]
  df_long <- tidyr::pivot_longer(df, -theta, names_to = "Series", values_to = "Density")
  ggplot(df_long, aes(theta, Density, color = Series)) +
    geom_line(linewidth = 1) +
    labs(title = title, subtitle = subtitle, x = expression(theta), y = expression(varphi(theta))) +
    theme_minimal(base_size = 12)
}

plot_surface_heatmap <- function(Sgrid, time, theta, title = "") {
  df <- expand.grid(time = time, theta = theta)
  df$value <- as.vector(Sgrid)
  ggplot(df, aes(theta, time, fill = value)) +
    geom_raster(interpolate = FALSE) +
    scale_y_continuous(expand = c(0, 0)) +
    labs(title = title, x = expression(theta), y = "time") +
    theme_minimal(base_size = 12)
}

plot_phase_metric_time <- function(metric_df, title = "") {
  df <- metric_df |>
    tidyr::pivot_longer(-time, names_to = "Metric", values_to = "Value")
  ggplot(df, aes(time, Value, color = Metric)) +
    geom_line(linewidth = 0.9) +
    labs(title = title, x = "time", y = "metric") +
    theme_minimal(base_size = 12)
}

plot_predictive_overlay <- function(Y, time, pred_list, sample_frac = 0.2) {
  n <- nrow(Y)
  K <- ncol(Y)
  obs_df <- expand.grid(rep = seq_len(n), time = time)
  obs_df$y <- as.vector(Y)
  if (sample_frac < 1) {
    keep <- sample.int(nrow(obs_df), size = max(1L, floor(sample_frac * nrow(obs_df))))
    obs_df <- obs_df[keep, ]
  }

  pred_df <- purrr::imap_dfr(pred_list, function(df, nm) {
    tibble(
      time = df$time,
      mean = df$mean,
      q05 = df$q0.05,
      q25 = df$q0.25,
      q75 = df$q0.75,
      q95 = df$q0.95,
      Method = nm
    )
  })

  ggplot() +
    geom_point(data = obs_df, aes(time, y), alpha = 0.18, shape = 1) +
    geom_ribbon(data = pred_df, aes(time, ymin = q05, ymax = q95, fill = Method), alpha = 0.15) +
    geom_ribbon(data = pred_df, aes(time, ymin = q25, ymax = q75, fill = Method), alpha = 0.25) +
    geom_line(data = pred_df, aes(time, mean, color = Method), linewidth = 1) +
    labs(title = "Posterior predictive overlays", x = "time", y = "response") +
    theme_minimal(base_size = 12)
}

4 Statistical Inference

The revised paper recommends bootstrap-based uncertainty quantification whenever the phase projection or weak spectral identifiability may materially affect finite-sample inference. The following helper implements a row bootstrap: entire replicate trajectories are resampled with replacement, the estimator is re-fitted, and the anchored phase laws are collected.

This is computationally expensive and is therefore provided as an optional protocol component. It is not run by default in the full experiments below.

bootstrap_phase_bands <- function(Y, time, fit_fun, B = 25, seed = 123, probs = c(0.025, 0.5, 0.975), ...) {
  set.seed(seed)
  N <- nrow(Y)
  fits <- vector("list", B)
  for (b in seq_len(B)) {
    idx <- sample.int(N, size = N, replace = TRUE)
    fits[[b]] <- fit_fun(Y[idx, , drop = FALSE], time = time, ...)
  }

  phi_array <- simplify2array(lapply(fits, function(f) f$anchored$varphi_theta))
  # [K, L, B]
  bands <- apply(phi_array, c(1, 2), quantile, probs = probs, na.rm = TRUE)
  list(fits = fits, bands = bands, probs = probs)
}

Example usage for local uncertainty studies:

# Example (not run by default because it can be slow):
# # # Generate data
# # ds <- simulate_double_slit_aligned(N = 100, K = 50)
# # ds <- simulate_double_slit()
# 
# ds_boot <- bootstrap_phase_bands(
#   Y = ds$Y,
#   time = ds$time,
#   fit_fun = fit_kpt_gem,
#   B = 50,
#   J = 16,
#   L_theta = 256,
#   lambda_phi = 5e-3,
#   p_phi = 1,
#   lambda_F = 5e-3,
#   p_F = 1,
#   sigma = 0.15,
#   max_iter = 30,
#   tol = 1e-3,
#   verbose = FALSE
# )

5 Applications

5.1 Simulation study I: phase-modulated double-slit experiment

5.1.1 Double-slit simulator

dvonmises <- function(theta, mu, kappa) {
  exp(kappa * cos(theta - mu)) / (2*pi * besselI(kappa, 0))
}

simulate_double_slit <- function(N = 100, K = 50, L_theta = 256,
                                 kappa_phase = 4, noise_sd = 0.15,
                                 kappa_wave = 12, seed = 2025) {
  set.seed(seed)
  theta <- theta_grid(L_theta)
  time <- seq(0, 1, length.out = K)

  S_true <- matrix(0, nrow = K, ncol = L_theta)
  phi_true <- matrix(0, nrow = K, ncol = L_theta)

  for (k in seq_len(K)) {
    t <- time[k]
    envelope <- 0.8 * exp(-20 * (t - 0.5)^2)
    baseline <- 0.2 * sin(2*pi * t)
    S_true[k, ] <- baseline +
      envelope * cos(kappa_wave * t + theta) +
      0.24 * envelope * cos(2 * (kappa_wave * t + theta) + pi / 6)

    mu1 <- (2*pi * (0.2 + 0.6 * t)) %% (2*pi)
    mu2 <- (2*pi * (0.9 - 0.6 * t)) %% (2*pi)
    w <- 0.5 + 0.4 * sin(2*pi * t)
    p1 <- dvonmises(theta, mu1, kappa_phase)
    p2 <- dvonmises(theta, mu2, kappa_phase)
    dens <- w * p1 + (1 - w) * p2
    phi_true[k, ] <- dens / mean(dens)
  }

  Y <- matrix(0, nrow = N, ncol = K)
  for (k in seq_len(K)) {
    probs <- normalize_rows_prob(matrix(phi_true[k, ], nrow = 1L))[1L, ]
    idx <- sample.int(L_theta, size = N, replace = TRUE, prob = probs)
    Y[, k] <- S_true[k, idx] + rnorm(N, 0, noise_sd)
  }

  truth_anchor <- anchor_by_first_moment(phi_true, theta)
  list(
    Y = Y,
    time = time,
    theta = theta,
    phi_true = phi_true,
    S_true = S_true,
    phi_true_anchored = truth_anchor$phi,
    S_true_anchored = apply_anchor_to_surface(S_true, truth_anchor$delta),
    params = list(
      N = N, K = K, L_theta = L_theta,
      kappa_phase = kappa_phase, noise_sd = noise_sd,
      kappa_wave = kappa_wave, seed = seed
    )
  )
}

ds <- simulate_double_slit()

5.1.2 Fit KPT-GEM and KPT-FFT to the double-slit data

fit_ds_gem <- fit_kpt_gem(
  Y = ds$Y,
  time = ds$time,
  J = 16,
  L_theta = length(ds$theta),
  lambda_phi = 5e-3,
  p_phi = 1,
  lambda_F = 5e-3,
  p_F = 1,
  sigma = 0.15,
  max_iter = 30,
  tol = 1e-3,
  track_loglik = FALSE,
  verbose = TRUE
)

fit_ds_fft <- fit_kpt_fft(
  Y = ds$Y,
  time = ds$time,
  J = 16,
  L_theta = length(ds$theta),
  lambda_phi = 5e-3,
  p_phi = 1,
  lambda_F = 5e-3,
  p_F = 1,
  sigma = 0.15,
  max_iter = 30,
  tol = 1e-3,
  track_loglik = FALSE,
  verbose = TRUE
)

5.1.3 Double-slit evaluation after alignment

phase_ds_gem <- phase_metrics_df(ds$phi_true_anchored, fit_ds_gem$anchored$varphi_theta, ds$time)
phase_ds_fft <- phase_metrics_df(ds$phi_true_anchored, fit_ds_fft$anchored$varphi_theta, ds$time)

surf_ds_gem <- surface_L2_cone(ds$S_true_anchored, fit_ds_gem$anchored$Sgrid, ds$time)
surf_ds_fft <- surface_L2_cone(ds$S_true_anchored, fit_ds_fft$anchored$Sgrid, ds$time)

summary_ds <- tibble(
  Algorithm = c("KPT-GEM", "KPT-FFT"),
  Mean_JSD_bits = c(mean(phase_ds_gem$JSD_bits), mean(phase_ds_fft$JSD_bits)),
  Mean_Hellinger = c(mean(phase_ds_gem$Hellinger), mean(phase_ds_fft$Hellinger)),
  Mean_TV = c(mean(phase_ds_gem$TV), mean(phase_ds_fft$TV)),
  Rel_L2_surface = c(surf_ds_gem$rel_L2, surf_ds_fft$rel_L2)
)

knitr::kable(summary_ds, digits = 4, caption = "Double-slit experiment: aligned reconstruction metrics.")
Double-slit experiment: aligned reconstruction metrics.
Algorithm Mean_JSD_bits Mean_Hellinger Mean_TV Rel_L2_surface
KPT-GEM 0.2629 0.4474 0.4861 1.0896
KPT-FFT 0.5213 0.6490 0.6846 1.9756

5.1.4 Double-slit representative phase laws and surfaces

rep_idx_ds <- c(1, 13, 25, 37, 50)
for (k in rep_idx_ds) {
  print(
    plot_phase_slice(
      theta = ds$theta,
      truth = ds$phi_true_anchored[k, ],
      est_list = list(
        `KPT-GEM` = fit_ds_gem$anchored$varphi_theta[k, ],
        `KPT-FFT` = fit_ds_fft$anchored$varphi_theta[k, ]
      ),
      title = "Double-slit phase recovery",
      subtitle = sprintf("time = %.3f", ds$time[k])
    )
  )
}

plot_surface_heatmap(ds$S_true_anchored, ds$time, ds$theta, title = "Double-slit truth: anchored kime surface")

plot_surface_heatmap(fit_ds_gem$anchored$Sgrid, ds$time, ds$theta, title = "Double-slit KPT-GEM: anchored estimated surface")

plot_surface_heatmap(fit_ds_fft$anchored$Sgrid, ds$time, ds$theta, title = "Double-slit KPT-FFT: anchored estimated surface")

5.1.5 Double-slit phase metrics over time

plot_phase_metric_time(phase_ds_gem, title = "Double-slit KPT-GEM phase metrics versus time")

plot_phase_metric_time(phase_ds_fft, title = "Double-slit KPT-FFT phase metrics versus time")

5.1.6 Posterior predictive time overlays

pred_ds_gem <- kpt_predictive_bands(fit_ds_gem, B = 3000, seed = 11)
pred_ds_fft <- kpt_predictive_bands(fit_ds_fft, B = 3000, seed = 12)
plot_predictive_overlay(ds$Y, ds$time, list(`KPT-GEM` = pred_ds_gem, `KPT-FFT` = pred_ds_fft), sample_frac = 0.25)

5.1.7 Predictive mapping from the kime representation to screen position

The revised paper keeps the double-slit screen-position mapping and now treats it as a proper normalized likelihood. The implementation below follows the paper directly.

ds_physics <- function(lambda = 532e-9,
                       d = 20e-6,
                       L = 0.5,
                       x_window = c(-0.02, 0.02),
                       nx = 800,
                       envelope_sigma = 0.005) {
  x <- seq(x_window[1], x_window[2], length.out = nx)
  envelope <- exp(-0.5 * (x / envelope_sigma)^2)
  alpha <- 2*pi * d / (lambda * L)
  list(
    lambda = lambda,
    d = d,
    L = L,
    x = x,
    envelope = envelope,
    alpha = alpha
  )
}

rectify_surface <- function(u) {
  pmax(u - min(u), 0)
}

screen_density_from_kpt <- function(Sbar, phibar, theta, physics, mu = 0, rectify = rectify_surface) {
  x <- physics$x
  alpha <- physics$alpha
  env <- physics$envelope
  phibar <- phibar / mean(phibar)

  intensity <- numeric(length(x))
  for (i in seq_along(x)) {
    theta_query <- (alpha * x[i] + mu - theta) %% (2*pi)
    sval <- periodic_interp(theta, Sbar, theta_query)
    intensity[i] <- env[i] * mean(rectify(sval) * phibar)
  }
  intensity <- pmax(intensity, 0)
  density <- intensity / max(trapz_vec(x, intensity), 1e-15)
  tibble(x = x, intensity = intensity, density = density)
}

simulate_screen_hits <- function(Sbar, phibar, theta, physics, n_hits = 1000, mu = 0, seed = 1) {
  set.seed(seed)
  dens_df <- screen_density_from_kpt(Sbar, phibar, theta, physics, mu = mu)
  probs <- dens_df$density / sum(dens_df$density)
  sample(dens_df$x, size = n_hits, replace = TRUE, prob = probs)
}

posterior_mu_discrete <- function(x_hits, Sbar, phibar, theta, physics, mu_grid = theta) {
  prior <- pmax(phibar, 1e-15)
  prior <- prior / sum(prior)

  logpost <- numeric(length(mu_grid))
  for (m in seq_along(mu_grid)) {
    dens_df <- screen_density_from_kpt(Sbar, phibar, theta, physics, mu = mu_grid[m])
    px <- approx(dens_df$x, dens_df$density, xout = x_hits, rule = 2L)$y
    logpost[m] <- log(prior[m]) + sum(log(pmax(px, 1e-300)))
  }
  logpost <- logpost - max(logpost)
  post <- exp(logpost)
  post <- post / sum(post)
  tibble(mu = mu_grid, posterior = post)
}

plot_screen_density <- function(hit_df, density_list) {
  dens_df <- purrr::imap_dfr(density_list, function(df, nm) dplyr::mutate(df, Method = nm))
  ggplot() +
    geom_histogram(data = hit_df, aes(x, y = after_stat(density)), bins = 60,
                   fill = "grey80", color = "grey50", alpha = 0.5) +
    geom_line(data = dens_df, aes(x, density, color = Method), linewidth = 1) +
    labs(title = "Double-slit screen-position density", x = "screen position x", y = "density") +
    theme_minimal(base_size = 12)
}

Below we run the double-slit screen analysis.

physics_ds <- ds_physics()

Sbar_true <- colMeans(ds$S_true_anchored)
phibar_true <- colMeans(ds$phi_true_anchored)
Sbar_gem  <- colMeans(fit_ds_gem$anchored$Sgrid)
phibar_gem <- colMeans(fit_ds_gem$anchored$varphi_theta)
Sbar_fft  <- colMeans(fit_ds_fft$anchored$Sgrid)
phibar_fft <- colMeans(fit_ds_fft$anchored$varphi_theta)

screen_true <- screen_density_from_kpt(Sbar_true, phibar_true, ds$theta, physics_ds, mu = 0)
screen_gem  <- screen_density_from_kpt(Sbar_gem, phibar_gem, ds$theta, physics_ds, mu = 0)
screen_fft  <- screen_density_from_kpt(Sbar_fft, phibar_fft, ds$theta, physics_ds, mu = 0)

x_hits <- simulate_screen_hits(Sbar_true, phibar_true, ds$theta, physics_ds, n_hits = 1500, mu = 0, seed = 99)
hit_df <- tibble(x = x_hits)

plot_screen_density(hit_df, list(Truth = screen_true, `KPT-GEM` = screen_gem, `KPT-FFT` = screen_fft))

Next, we estimate and plot the posterior distribution for the source offset \(\mu\).

mu_post_gem <- posterior_mu_discrete(x_hits, Sbar_gem, phibar_gem, ds$theta, physics_ds)
mu_post_fft <- posterior_mu_discrete(x_hits, Sbar_fft, phibar_fft, ds$theta, physics_ds)

mu_df <- bind_rows(
  mutate(mu_post_gem, Method = "KPT-GEM"),
  mutate(mu_post_fft, Method = "KPT-FFT")
)

ggplot(mu_df, aes(mu, posterior, color = Method)) +
  geom_line(linewidth = 1) +
  labs(title = "Posterior distribution for the source offset mu", x = expression(mu), y = "posterior probability") +
  theme_minimal(base_size = 12)

6 Conclusions

This revised protocol resolves the main issues identified in the previous notebook version.

  1. It uses the time-specific latent phase model everywhere.
  2. It enforces the correct normalization convention based on \(d\theta/(2\pi)\) and the discrete row-mean-equals-one rule.
  3. It distinguishes the exact EM reference algorithm from the practical projected routines used numerically.
  4. It fixes the earlier convergence-diagnostic bug by computing differences before overwriting state variables.
  5. It replaces the earlier truncation-and-renormalization rule with the exact Euclidean projection onto the discrete simplex.
  6. It stores and uses sigma consistently through the predictive functions.
  7. It implements the phase-dependent fMRI simulation used in the revised paper.
  8. It uses one synchronized simulation-parameter table matching the manuscript.
  9. It keeps the double-slit screen-position mapping and turns it into a proper normalized likelihood, including a posterior for the source offset \(\mu\).
  10. It removes unrelated experimental debris and external dependencies so the protocol is self-contained and locally executable.

V.6 of hte KPT notebook provides a clean computational companion to the revised paper. For more realistic evaluation, we can add larger-scale bootstrap studies, repeated-seed Monte Carlo summaries, or sensitivity analyses for the regularization parameters.

SOCR Resource Visitor number Web Analytics SOCR Email