diff --git a/ptgp/conditionals.py b/ptgp/conditionals.py
index b3fb84c..aeba4ca 100644
--- a/ptgp/conditionals.py
+++ b/ptgp/conditionals.py
@@ -55,7 +55,7 @@ def conditional_unwhitened(A, Kmn, Knn, f, q_sqrt=None, full_cov=False):
 def base_conditional(Kmn, Kmm, Knn, f, q_sqrt=None, white=False, full_cov=False):
     """Back-compat wrapper. Materialises Kmm; use the helpers above directly
     when you have a structured Kuu_solve / Kuu_sqrt_solve."""
-    # Add jitter to keep Kmm PSD under float noise — matches GPflow / PyMC default.
+    # Add jitter to keep Kmm PSD under float noise; matches GPflow / PyMC default.
     # Re-annotate after the addition: PyTensor canonicalizes ``Kmm + c·I`` into a
     # ``set_subtensor`` on the diagonal, which our PSD-inference rules don't see
     # through. The mathematical identity (PSD + c·I PSD ⇒ PSD) is sound.
diff --git a/ptgp/gp/svgp.py b/ptgp/gp/svgp.py
index cdfa52b..84c9ed3 100644
--- a/ptgp/gp/svgp.py
+++ b/ptgp/gp/svgp.py
@@ -77,7 +77,7 @@ def init_variational_params(M, q_mu_init=None, q_sqrt_init=None):
 
     ``q_sqrt`` is stored as a flat vector of length ``M·(M+1)/2`` and
     materialised by filling a lower-triangular matrix and applying
-    ``pt.softplus`` to the diagonal — guarantees a true Cholesky factor
+    ``pt.softplus`` to the diagonal, which guarantees a true Cholesky factor
     (lower-triangular with strictly positive diagonal) at every optimizer step.
     ``q_mu`` is stored as a length-``M`` flat vector.
 
@@ -203,7 +203,7 @@ def extra_init(self):
     def predict_marginal(self, X, incl_lik=False):
         """Posterior marginal mean and variance at each point in X.
 
-        Returns the per-point posterior — correlations between test
+        Returns the per-point posterior. Correlations between test
         points are discarded. Use ``predict_joint`` for the full (N, N)
         covariance or ``predict_f_samples`` to draw smooth function samples.
 
@@ -265,7 +265,7 @@ def predict_f_samples(self, X, epsilon, jitter=1e-6):
         """Draw samples of the latent f at X from the joint posterior.
 
         Samples are produced via a Cholesky transform of caller-supplied
-        iid-standard-normal noise. The caller owns the RNG — the function
+        iid-standard-normal noise. The caller owns the RNG; the function
         itself is deterministic.
 
         Parameters
diff --git a/ptgp/gp/vfe.py b/ptgp/gp/vfe.py
index 3c2550a..1d6f155 100644
--- a/ptgp/gp/vfe.py
+++ b/ptgp/gp/vfe.py
@@ -8,7 +8,7 @@
 class VFE:
     """Variational Free Energy (SGPR) sparse Gaussian Process.
 
-    Uses Titsias' collapsed bound — inducing variables are analytically
+    Uses Titsias' collapsed bound; inducing variables are analytically
     integrated out. The observation model is Gaussian; parameterize the noise
     via ``sigma``.
 
diff --git a/ptgp/idata.py b/ptgp/idata.py
index c1954d5..d0248fb 100644
--- a/ptgp/idata.py
+++ b/ptgp/idata.py
@@ -23,18 +23,18 @@ def to_idata(
 
     Groups produced (empty groups are omitted):
 
-    - ``point_estimate`` — every optimized value in constrained space, on its
+    - ``point_estimate``: every optimized value in constrained space, on its
       own dims with no ``chain``/``draw`` axes (nothing was sampled). Includes
       PyMC RVs (``sigma``, ``eta``, …) *and* non-PyMC extras (``Z``, ``q_mu``,
       ``q_sqrt_flat``, …).
-    - ``unconstrained_point_estimate`` — the same values in value-var
+    - ``unconstrained_point_estimate``: the same values in value-var
       (transformed) space.
-    - ``optimizer_result`` — scalar fields from the scipy ``OptimizeResult``
+    - ``optimizer_result``: scalar fields from the scipy ``OptimizeResult``
       (``fun``, ``success``, ``status``, ``message``, ``nit``, ``nfev``,
       ``njev``) plus per-iteration trajectory DataArrays from ``history``
       (``elbo``, ``trace_penalty``, …) on an ``iteration`` dimension. There is
       no ``sample_stats``: nothing was sampled.
-    - ``observed_data`` / ``constant_data`` — pulled from the PyMC model.
+    - ``observed_data`` / ``constant_data``: pulled from the PyMC model.
 
     Parameters
     ----------
@@ -103,7 +103,7 @@ def _collect_optimized(model, shared_params, shared_extras):
 
     PyMC RVs contribute to both dicts: the constrained dict uses ``rv.name`` and
     applies the backward transform; the unconstrained dict uses ``vv.name`` and
-    keeps the value-var-space value. Non-PyMC extras have no transform — they
+    keeps the value-var-space value. Non-PyMC extras have no transform, so they
     appear in both dicts under ``sv.name`` so callers can find them in either
     group.
     """
@@ -154,7 +154,7 @@ def _build_optimizer_result(result, history, phase_labels):
     populates ``elbo``, ``trace_penalty``, …).
 
     If a history field name collides with a result-scalar name, the trajectory
-    wins — it's strictly more informative than the terminal value.
+    wins, it's strictly more informative than the terminal value.
     """
     data_vars = {}
     coords = {}
diff --git a/ptgp/inducing.py b/ptgp/inducing.py
index 4f9192f..53101b8 100644
--- a/ptgp/inducing.py
+++ b/ptgp/inducing.py
@@ -111,9 +111,9 @@ class KernelHealthDiagnostics:
         Per-data-point residual conditional variance ``Kff_diag - Q_diag``.
         Large values identify points poorly covered by ``Z``.
     total_variance : float
-        ``tr(Kff)`` — the kernel diagonal sum on ``X`` before conditioning.
+        ``tr(Kff)``, the kernel diagonal sum on ``X`` before conditioning.
     nystrom_residual : float
-        ``sum(d_final) = tr(Kff - Q)`` — total unexplained variance.
+        ``sum(d_final) = tr(Kff - Q)``, total unexplained variance.
     kuu_min_eigenvalue : float
         Smallest eigenvalue of ``K(Z, Z) + jitter * I``.
     kuu_max_eigenvalue : float
@@ -159,7 +159,7 @@ class RandomSubsampleDiagnostics:
         The ``M`` argument.
     M_returned : int
         Rows in the returned ``Z``. Always equals ``M_requested`` for this
-        routine — kept for API symmetry with the other init routines.
+        routine, kept for API symmetry with the other init routines.
     N_candidates : int
         Rows in the input ``X``.
     n_unique : int
@@ -223,7 +223,7 @@ class KMeansDiagnostics:
     dedup_tol : float
         The Euclidean-distance threshold used for deduplication.
     inertia : float
-        ``sum_i ||X[i] - centroid[label[i]]||^2`` — k-means' standard
+        ``sum_i ||X[i] - centroid[label[i]]||^2``, k-means' standard
         within-cluster sum of squares. Smaller is tighter.
     pairwise_min_distance : float
         Minimum Euclidean distance between returned centroids. ``nan`` if
@@ -529,7 +529,7 @@ def greedy_variance_init(
     Implements the "ConditionalVariance" initialization of Burt et al. (2020),
     *Convergence of Sparse Variational Inference in GP Regression*. At each
     step, the next inducing point is the row of ``X`` with largest remaining
-    conditional variance given the already-selected points — equivalent to
+    conditional variance given the already-selected points, equivalent to
     running a partial pivoted Cholesky decomposition of ``K(X, X)`` with the
     standard max-diagonal pivot rule. Selected points are a **subset of X**;
     this is discrete subset selection, not continuous optimization.
@@ -539,7 +539,7 @@ def greedy_variance_init(
     Recommended workflow
     --------------------
     Burt et al. show that with a good greedy initialization, ``Z`` typically
-    does **not** need to be gradient-optimized during training — for most
+    does **not** need to be gradient-optimized during training. For most
     problems the frozen subset is within noise of jointly-optimized ``Z`` at a
     tiny fraction of the compute. The standard recipe:
 
@@ -682,7 +682,7 @@ def compute_inducing_diagnostics(
     Computes the same kernel-derived metrics that
     :func:`greedy_variance_init` reports (Kuu eigenvalues, per-point
     residual conditional variance, total variance), but for a ``Z`` that
-    you already have — from k-means, manual placement, post-training, or
+    you already have, from k-means, manual placement, post-training, or
     elsewhere. No greedy selection is performed; no ``trace_curve``
     history is available.
 
diff --git a/ptgp/inducing_fourier.py b/ptgp/inducing_fourier.py
index bac59d9..755d0f1 100644
--- a/ptgp/inducing_fourier.py
+++ b/ptgp/inducing_fourier.py
@@ -360,7 +360,7 @@ def Kuu_logdet(self, kernel):
     def Kuu_sqrt_solve(self, kernel, rhs):
         """Apply ``R^{-1}`` where ``R @ R.T = Kuu = diag(d) + U @ U.T``. ``rhs`` is ``(M, K)``.
 
-        Uses ``delta = -1 / (sqrt(1 + lam) * (1 + sqrt(1 + lam)))`` — algebraically
+        Uses ``delta = -1 / (sqrt(1 + lam) * (1 + sqrt(1 + lam)))``, algebraically
         equivalent to ``1/sqrt(1+lam) - 1`` divided by ``lam`` but with no division by
         ``lam`` or ``sqrt(lam)``, so it stays finite as ``lam -> 0``.
         """
diff --git a/ptgp/kernels/base.py b/ptgp/kernels/base.py
index e51ae43..760236b 100644
--- a/ptgp/kernels/base.py
+++ b/ptgp/kernels/base.py
@@ -38,7 +38,7 @@ def __call__(self, X, Y=None):
         return self._eval(X, Y)
 
     def _eval(self, X, Y):
-        """Raw kernel matrix — subclasses implement."""
+        """Raw kernel matrix; subclasses implement."""
         raise NotImplementedError
 
     def diag(self, X):
diff --git a/ptgp/kernels/categorical.py b/ptgp/kernels/categorical.py
index ca2c848..e271703 100644
--- a/ptgp/kernels/categorical.py
+++ b/ptgp/kernels/categorical.py
@@ -10,7 +10,7 @@ class Overlap(Kernel):
 
     where D is the number of active categorical columns. The kernel returns
     the fraction of active columns whose levels match. With one active column,
-    this is simply ``1[x == y]``.
+    this is ``1[x == y]``.
 
     Values are expected to be non-negative integer level codes stored in a
     float matrix (the library-wide convention). Entries are cast to ``int64``
@@ -50,7 +50,7 @@ class LowRankCategorical(Kernel):
     kernel. Operates on a single categorical column. Level codes are expected
     to be non-negative integers stored in a float matrix (the library-wide
     convention). Codes are cast to ``int64`` inside ``_eval`` and ``diag``;
-    they are not range-checked — passing a code >= ``num_levels`` is user
+    they are not range-checked; passing a code >= ``num_levels`` is user
     error.
 
     Parameters
diff --git a/ptgp/kernels/combination.py b/ptgp/kernels/combination.py
index 62a6356..f9819fa 100644
--- a/ptgp/kernels/combination.py
+++ b/ptgp/kernels/combination.py
@@ -27,7 +27,7 @@ def _eval(self, X, Y):
         return self.k1(X, Y) + self.k2(X, Y)
 
     def diag(self, X):
-        """Return the per-point variance — sum of the two component diagonals."""
+        """Per-point variance: sum of the two component diagonals."""
         return self.k1.diag(X) + self.k2.diag(X)
 
 
@@ -56,7 +56,7 @@ def _eval(self, X, Y):
         return self.k1 * self.k2(X, Y)
 
     def diag(self, X):
-        """Return the per-point variance — product of component diagonals, with scalar passthrough."""
+        """Per-point variance: product of component diagonals, with scalar passthrough."""
         k1_is_kernel = isinstance(self.k1, Kernel)
         k2_is_kernel = isinstance(self.k2, Kernel)
         if k1_is_kernel and k2_is_kernel:
diff --git a/ptgp/kernels/nonstationary.py b/ptgp/kernels/nonstationary.py
index 94d443b..f9526f2 100644
--- a/ptgp/kernels/nonstationary.py
+++ b/ptgp/kernels/nonstationary.py
@@ -136,7 +136,7 @@ class WarpedInput(Kernel):
     ----------
     input_dim : int
         Number of columns of ``X`` this kernel expects. May differ from
-        ``kernel_func.input_dim`` — the warp may change dimensionality.
+        ``kernel_func.input_dim``; the warp may change dimensionality.
     kernel_func : Kernel
         Inner kernel applied to the warped inputs.
     warp_func : callable
diff --git a/ptgp/likelihoods/base.py b/ptgp/likelihoods/base.py
index 24771d6..4842018 100644
--- a/ptgp/likelihoods/base.py
+++ b/ptgp/likelihoods/base.py
@@ -60,7 +60,7 @@ def predict_mean_and_var(self, mu, var):
         return E_mean, E_var + E_mean_sq - E_mean**2
 
     def predict_log_density(self, y, mu, var):
-        """log E_{q(f)}[p(y|f)] — predictive log-density at test points.
+        """log E_{q(f)}[p(y|f)], the predictive log-density at test points.
 
         Default: Gauss-Hermite quadrature in log-space for numerical stability.
         """
diff --git a/ptgp/likelihoods/gaussian.py b/ptgp/likelihoods/gaussian.py
index ef2cb8a..dbdcca7 100644
--- a/ptgp/likelihoods/gaussian.py
+++ b/ptgp/likelihoods/gaussian.py
@@ -33,7 +33,7 @@ def sigma_at(self, X_train, X_new):
 
         Sigma is treated as data-dependent iff it has any free symbolic
         tensor input in its graph (covers both raw ``pt.matrix`` X and
-        ``pm.Data`` X — the latter is a SharedVariable but its type is
+        ``pm.Data`` X, the latter is a SharedVariable but its type is
         ``TensorType``). For data-independent sigma the call is a no-op.
         For data-dependent sigma, strict-mode ``graph_replace`` raises a
         clear error if ``X_train`` is not in sigma's graph.
diff --git a/ptgp/linalg/operator.py b/ptgp/linalg/operator.py
index ce97521..635d20a 100644
--- a/ptgp/linalg/operator.py
+++ b/ptgp/linalg/operator.py
@@ -1,4 +1,4 @@
-# LinearOperatorType — placeholder for future PyTensor Op implementation.
+# LinearOperatorType: placeholder for a future PyTensor Op.
 # Will define a new PyTensor type representing a matrix implicitly via
 # its matrix-vector product, enabling CG solves and Lanczos log-det
 # without materialising the full N x N kernel matrix.
diff --git a/ptgp/linalg/ops.py b/ptgp/linalg/ops.py
index d0f3154..ffc1283 100644
--- a/ptgp/linalg/ops.py
+++ b/ptgp/linalg/ops.py
@@ -1,5 +1,5 @@
 # Custom PyTensor Ops for lazy linear algebra.
-# KernelLinearOp  — represents K(X, X) as a LinearOperatorType
-# LinearOpMatvec  — computes A @ v without materialising A
-# LinearOpSolve   — CG-based solve: A x = b
-# LinearOpLogdet  — Lanczos-based log-determinant estimator
+# KernelLinearOp: represents K(X, X) as a LinearOperatorType
+# LinearOpMatvec: computes A @ v without materialising A
+# LinearOpSolve: CG-based solve, A x = b
+# LinearOpLogdet: Lanczos-based log-determinant estimator
diff --git a/ptgp/objectives.py b/ptgp/objectives.py
index 32c0f4f..fffece9 100644
--- a/ptgp/objectives.py
+++ b/ptgp/objectives.py
@@ -153,7 +153,7 @@ def collapsed_elbo(vfe, X, y):
 def fitc_log_marginal_likelihood(vfe, X, y):
     """FITC (Fully Independent Training Conditional) approximate log marginal likelihood.
 
-    Unlike ``collapsed_elbo``, FITC is not a lower bound — it approximates the
+    Unlike ``collapsed_elbo``, FITC is not a lower bound; it approximates the
     log marginal likelihood using the true per-point diagonal rather than the
     Nystrom diagonal throughout. The FITC covariance is::
 
@@ -185,9 +185,9 @@ def fitc_log_marginal_likelihood(vfe, X, y):
     Returns
     -------
     FITCTerms
-        ``fitc`` — FITC approximate log marginal likelihood (fit + logdet).
-        ``fit`` — quadratic term: ``-0.5 * (y^T K_fitc^{-1} y + N log 2π)``.
-        ``logdet`` — log-determinant term: ``-0.5 log|K_fitc|``.
+        ``fitc``: FITC approximate log marginal likelihood (fit + logdet).
+        ``fit``: quadratic term: ``-0.5 * (y^T K_fitc^{-1} y + N log 2π)``.
+        ``logdet``: log-determinant term: ``-0.5 log|K_fitc|``.
     """
     sigma2 = vfe.likelihood.sigma**2
     N = X.shape[0]
@@ -294,20 +294,20 @@ def vfe_diagnostics(vfe, X, y):
     """Collapsed ELBO terms plus sigma and two normalised fit metrics.
 
     Returns a ``VFEDiagnostics`` namedtuple of symbolic TensorVariables,
-    suitable for use with :func:`ptgp.optim.compile_scipy_diagnostics`.
+    for use with :func:`ptgp.optim.compile_scipy_diagnostics`.
 
     Fields
     ------
     elbo, fit, trace_penalty
         Direct from :func:`collapsed_elbo`.
     nystrom_residual
-        ``tr(Kff - Qff) / N`` — per-point Nyström approximation error.
+        ``tr(Kff - Qff) / N``, the per-point Nyström approximation error.
     sigma
         Likelihood noise (constrained space).
     fit_per_n
-        ``fit / N`` — scale-invariant data fit.
+        ``fit / N``, the per-point data fit.
     excess_fit_per_n
-        ``fit_per_n + 0.5 * log(2π σ²)`` — how much better than noise floor.
+        ``fit_per_n + 0.5 * log(2π σ²)``, how much better than the noise floor.
         Goes to zero when the model fits at the noise level only.
     """
     terms = collapsed_elbo(vfe, X, y)
diff --git a/ptgp/optim/api.py b/ptgp/optim/api.py
index e82de13..3ccaf33 100644
--- a/ptgp/optim/api.py
+++ b/ptgp/optim/api.py
@@ -49,7 +49,7 @@ def fit(
     compile_kwargs=None,
     **scipy_kwargs,
 ):
-    """Batteries-included estimation method for GP models.
+    """Estimation method for GP models.
 
     Estimation uses reasonable defaults based on the provided approximation. For
     fine-grained control, use the low-level API.
@@ -70,18 +70,18 @@ def fit(
         ``gp_model.default_objective``.
     method : str
         Optimization method passed to :func:`scipy.optimize.minimize`. Must
-        accept a Jacobian — ``fit`` always supplies one (``jac=True``).
+        accept a Jacobian; ``fit`` always supplies one (``jac=True``).
         Default ``"L-BFGS-B"``.
     init : str
         Strategy for the initial parameter vector. One of:
 
-        - ``"prior_median"`` (default) — median of each prior via
+        - ``"prior_median"`` (default): median of each prior via
           ``pm.icdf(rv, 0.5)``, with a 500-sample fallback when icdf is
           unimplemented and a per-RV fallback to PyMC's initial point for
           improper priors.
-        - ``"prior_draw"`` — one draw from each prior. Stochastic; pin with
+        - ``"prior_draw"``: one draw from each prior. Stochastic; pin with
           ``init_rng`` for reproducibility.
-        - ``"unconstrained_zero"`` — PyMC's ``model.initial_point()`` (0 in
+        - ``"unconstrained_zero"``: PyMC's ``model.initial_point()`` (0 in
           unconstrained space unless ``initval=`` was set on the RV).
     init_rng : int or numpy.random.Generator, optional
         Seed for the sampling-based portions of ``"prior_median"`` and
@@ -93,7 +93,7 @@ def fit(
     **scipy_kwargs
         Additional keyword arguments forwarded to
         :func:`scipy.optimize.minimize` (e.g. ``tol``, ``options``,
-        ``bounds``). Do not pass ``jac`` — it is set internally.
+        ``bounds``). Do not pass ``jac``; it is set internally.
 
     Returns
     -------
diff --git a/ptgp/optim/training.py b/ptgp/optim/training.py
index c85401d..2724113 100644
--- a/ptgp/optim/training.py
+++ b/ptgp/optim/training.py
@@ -1,7 +1,7 @@
 """Training and prediction compilation for PTGP models with PyMC priors.
 
 Uses ``pytensor.shared`` variables so that training automatically updates
-the parameters used by prediction — no model reconstruction needed.
+the parameters used by prediction; no model reconstruction needed.
 """
 
 import logging
@@ -96,7 +96,7 @@ def _make_initial_point(model, init="prior_median", rng=None, n_median_samples=5
     Some priors cannot be sampled or quantiled. Improper priors like
     ``pm.HalfFlat`` and ``pm.Flat`` raise ``NotImplementedError`` from both
     ``pm.icdf`` and ``pm.draw`` because they have no proper measure. When
-    every method fails for a given RV — or returns non-finite values — that
+    every method fails for a given RV, or returns non-finite values, that
     RV's value is taken from ``model.initial_point()`` instead. PyMC's
     initial point is 0 in unconstrained space for every parameter unless
     ``initval`` was set explicitly, which corresponds to:
@@ -185,7 +185,7 @@ def _make_shared_params(
     If ``frozen_vars`` is given, value vars appearing as keys are initialized
     to their freeze values rather than the PyMC initial point. Without this,
     a frozen var's shared slot (and its scipy theta slice) would carry the
-    initial point — wrong for diagnostics and for ``unpack_to_shared``.
+    initial point, wrong for diagnostics and for ``unpack_to_shared``.
 
     Returns
     -------
@@ -275,7 +275,7 @@ def compile_training_step(
     model : pm.Model, optional
         PyMC model. Uses the enclosing ``with pm.Model()`` context if
         None. Every continuous free RV in the model is automatically
-        made into a trainable shared variable — you do not need to list
+        made into a trainable shared variable, so you do not need to list
         them.
     optimizer_fn : callable, optional
         Optimizer function (default: ``adam``). Must have signature
@@ -320,7 +320,7 @@ def compile_training_step(
         ``(X_batch, y_batch) -> loss_value``. Updates shared parameters
         in place.
     shared_params : dict
-        ``{value_var: shared_var}`` — the shared variables holding the
+        ``{value_var: shared_var}``, the shared variables holding the
         unconstrained parameter values. Needed by ``compile_predict``.
     shared_extras : list
         Shared variables for ``extra_vars``. Needed by ``compile_predict``.
@@ -440,7 +440,7 @@ def compile_scipy_objective(
         PTGP model whose hyperparameters are PyMC RVs.
     X_var : TensorVariable
         Symbolic input placeholder. ``X`` is passed to the compiled
-        function on each scipy iteration — typically the full training
+        function on each scipy iteration, typically the full training
         inputs for GP/VFE (batching is not used with quasi-Newton methods).
     y_var : TensorVariable
         Symbolic target placeholder, handled like ``X_var``.
@@ -516,7 +516,7 @@ def compile_scipy_objective(
     shared_params : dict
         ``{value_var: shared_var}`` for every continuous PyMC value var.
         Needed by :func:`compile_predict` and :func:`get_trained_params`.
-        Not read by ``fun`` — present only for the predict handoff.
+        Not read by ``fun``; present only for the predict handoff.
     shared_extras : list
         Shared variables for ``extra_vars``, in the same order. Needed
         by :func:`compile_predict`. Not read by ``fun``.
@@ -749,7 +749,7 @@ def tracked_minimize(fun, theta0, args, diag_fn=None, print_every=None, **scipy_
     theta0 : ndarray
         Initial parameter vector.
     args : tuple
-        Extra arguments forwarded to both ``fun`` and ``diag_fn`` — typically
+        Extra arguments forwarded to both ``fun`` and ``diag_fn``, typically
         ``(X, y)``.
     diag_fn : callable, optional
         ``(theta, *args) -> namedtuple`` from
@@ -924,17 +924,17 @@ def minimize_staged_vfe(
     the trace penalty instead of improving ``Z``. This function blocks that
     failure mode with a structured four-phase schedule:
 
-    1. **Phase 1** — freeze ``sigma`` at ``sigma_init``; train all other
+    1. **Phase 1**: freeze ``sigma`` at ``sigma_init``; train all other
        hyperparameters and (by default) ``Z`` together.
-    2. **Phase 2a** — freeze all hyperparameters; train ``Z`` only.
-    3. **Phase 2b** — freeze ``Z``; train all hyperparameters (sigma now free).
-    4. **Phase 3** — joint fine-tuning of everything (nothing frozen).
+    2. **Phase 2a**: freeze all hyperparameters; train ``Z`` only.
+    3. **Phase 2b**: freeze ``Z``; train all hyperparameters (sigma now free).
+    4. **Phase 3**: joint fine-tuning of everything (nothing frozen).
 
     Phases 2a/2b repeat ``phase2_cycles`` times. Pass ``phase2_cycles=0`` to
     skip straight from phase 1 to phase 3 when the ``Z`` initialisation is
     already good (e.g. from ``greedy_variance_init``).
 
-    ``sigma`` is auto-detected from ``gp_model.likelihood.sigma`` — no need to
+    ``sigma`` is auto-detected from ``gp_model.likelihood.sigma``; no need to
     pass a separate ``sigma_rv`` argument.
 
     All diagnostic history entries are :class:`ptgp.objectives.VFEDiagnostics`
@@ -977,7 +977,7 @@ def minimize_staged_vfe(
     init : str
         Initialisation strategy for phase 1 hyperparameters. Same options as
         :func:`compile_scipy_objective`: ``"prior_median"`` (default),
-        ``"prior_draw"``, or ``"unconstrained_zero"`` (PyMC initial point —
+        ``"prior_draw"``, or ``"unconstrained_zero"`` (PyMC initial point,
         0 in unconstrained space). Only affects phase 1; later phases
         inherit converged values from the previous phase.
     init_rng : int or numpy Generator, optional
diff --git a/ptgp/utils.py b/ptgp/utils.py
index d0dfc24..e62aede 100644
--- a/ptgp/utils.py
+++ b/ptgp/utils.py
@@ -12,7 +12,7 @@ def get_initial_params(model, init="prior_median", rng=None, n_median_samples=50
     """Return constrained-space initial values for all free RVs in a PyMC model.
 
     Uses the same initialization strategies as ``compile_scipy_objective``.
-    Useful for building proxy kernels with concrete float values to pass to
+    Use it to build proxy kernels with concrete float values to pass to
     ``greedy_variance_init``.
 
     Parameters