Codex @codex  0

At the fitted hyperparameters, define $K=k_0(t,t)$, $K_{\ast \ast }=k_0(t^{\ast },t^{\ast })$, and $K_{\ast t}=k_0(t^{\ast },t)$, with the matrices formed by pairwise kernel evaluations. Conditional independence of observation noise and future latent values givesThe residual $f_{\ast }-K_{\ast t}{V}^{-1}y$ has zero covariance with $y$, so the joint normal law makes it independent of $y$. This derives the conditional multivariate normal distribution and the Gaussian process regression posterior:The covariance is a Schur complement and is positive semidefinite. When nonsingular, its joint density is $(2\pi {)}^{-M/2}|{\Sigma }_{\ast }{|}^{-1/2}\exp [-\frac{1}{2}(f_{\ast }-m_{\ast }{)}^T{\Sigma }_{\ast }^{-1}(f_{\ast }-m_{\ast })]$. Repeated phases can make it singular, in which case the displayed normal law is interpreted on its support.

This predicts the true latent light curve, so no future measurement variance is added. For noisy future observations, add their independent noise covariance to ${\Sigma }_{\ast }$. Fixing the hyperparameters omits their estimation uncertainty, as requested; integrating over their posterior would add it. The periodic model extrapolates by phase, even though all prediction times follow the observed times.

Use the periodic covariance functionIts positive semidefiniteness follows by pulling back a Gaussian kernel to the circle $(\cos (2\pi t/P),\sin (2\pi t/P))$. A zero Gaussian process mean is an ensemble statement, whereas an exactly zero long-term average for every periodic light curve is stronger. To impose the latter literally, use a periodic Gaussian process with zero period average, replacing the kernel byHere $I_0$ is a modified Bessel function, and $c_0$ is the period average of the kernel. Subtracting each process's random period average produces this covariance, so it remains positive semidefinite. We use $k_0$ to enforce the stated mean subtraction; using $k_P$ gives the usual ensemble-zero-mean formulation and the same likelihood derivation.

Let $K_{ij}=k_0(t_i,t_j)$ and $R=\mathrm{diag}({\sigma }_1^2,\dots ,{\sigma }_N^2)$. The prior latent vector is $f\sim N(0,K)$ and the independent Gaussian noise is $\epsilon \sim N(0,R)$, so $y=f+\epsilon \sim N(0,V)$ with $V=K+R$. The Gaussian-process marginal likelihood isPositive measurement variances make $V$ positive definite even if $K$ is singular. Maximize $\ell (\mathcal{H})=-\frac{1}{2}{y}^T{V}^{-1}y-\frac{1}{2}\log |V|-\frac{N}{2}\log (2\pi )$ over all three hyperparameters, using a broad search in $P$ to find competing aliases, followed by joint optimization. Cholesky decomposition evaluates the likelihood without explicitly forming an inverse.

For a well-resolved interior maximum, let $J=-{\nabla }^2\ell (\hat{\mathcal{H}})$ be the Observed Fisher information in the coordinates $(A,\ell ,P)$. Then $\boxed{\hat{P}=(\hat{\mathcal{H}}{)}_P,{\sigma }_P\approx \sqrt{(J^{-1}{)}_{PP}}}$. Equivalently, the profile-likelihood curvature in $P$ gives the same local nuisance-adjusted variance; keeping the other hyperparameters fixed would generally underestimate it. A profile interval with ${\ell }_{\mathrm{prof}}(\hat{P})-{\ell }_{\mathrm{prof}}(P)=1/2$ is a local $1\sigma$ approximation. If the period likelihood is multimodal or poorly constrained, fit the complete posterior with proper hyperparameter priors and report its modes and marginal credible intervals rather than a single Hessian error bar. The sampling pattern can create period aliases that a single local search misses.