UPDE — Bayesian Uncertainty

Why Bayesian uncertainty exists in this project

Deterministic UPDE runs provide a single trajectory. Bayesian uncertainty surfaces how sensitive that trajectory is to uncertainty in estimated frequencies and couplings. For production planning and safety review, that difference is material.

The module exposes uncertainty in a structured way so teams can make explicit risk decisions rather than relying only on point estimates.

scpn_phase_orchestrator.upde.bayesian propagates uncertainty in natural frequencies omega and coupling matrices K_nm through the existing UPDE integrator. It reports posterior-predictive order parameter summaries as R ± sigma plus a configurable credible interval.

The production backend is deterministic NumPy Monte Carlo over explicit array distributions. numpyro and blackjax are reserved backend names and raise NotImplementedError until their samplers are implemented, benchmarked, and validated against the NumPy propagation baseline.

fit_gaussian_upde_posterior() provides the deterministic production baseline for posterior fitting from observed Kuramoto phase trajectories. It uses a finite-difference regression against the same UPDE coupling surface, enforces non-negative zero-diagonal coupling, emits JSON-safe diagnostics, and feeds the resulting Gaussian distributions directly into bayesian_upde_run().

audit_bayesian_backend_status() probes backend names through the same execution path. The NumPy backend must execute, while reserved sampler names such as numpyro and blackjax must fail closed with audit records until they have validated implementations and benchmark evidence.

Minimal Example

import numpy as np

from scpn_phase_orchestrator.upde import (
    BayesianUPDEConfig,
    GaussianArrayDistribution,
    bayesian_upde_run,
)

phases = np.array([0.0, 0.4, 1.1, 1.9])
omega_mean = np.array([0.9, 1.0, 1.08, 1.16])
knm_mean = np.full((4, 4), 0.18)
np.fill_diagonal(knm_mean, 0.0)
alpha = np.zeros((4, 4))

result = bayesian_upde_run(
    phases,
    omega=GaussianArrayDistribution(omega_mean, np.full(4, 0.015)),
    knm=GaussianArrayDistribution(
        knm_mean,
        np.full((4, 4), 0.01),
        non_negative=True,
        zero_diagonal=True,
    ),
    alpha=alpha,
    zeta=0.02,
    psi=0.1,
    config=BayesianUPDEConfig(n_samples=256, seed=7, n_steps=25),
)

r_mean, r_sigma = result.r_plus_minus

Semantics

• omega and knm may be deterministic arrays or distribution objects.
• GaussianArrayDistribution samples independent normal uncertainty per array entry and can enforce non-negative coupling and zero self-coupling.
• fit_gaussian_upde_posterior() estimates Gaussian omega and K_nm distributions from finite phase trajectories with explicit ridge and uncertainty floors.
• audit_bayesian_backend_status() records executable and fail-closed backend state for release and safety-review evidence.
• The existing UPDE kernel performs every rollout, so deterministic engine validation, phase wrapping, and backend dispatch semantics are preserved.
• BayesianUPDEResult.to_audit_record() emits JSON-safe uncertainty diagnostics suitable for safety review and replay logs.

How this is used operationally

Use Bayesian runs when estimates are sparse, noisy, or manually inferred. In that setting the output carries both central tendency and uncertainty so decision logic can enforce explicit safety thresholds before proposing changes.

Production interpretation

• Use Bayesian UPDE when uncertainty is itself a policy input, not an afterthought.
• The fail-closed backend status check is a compliance boundary: unavailable advanced samplers must not silently become “best effort” paths.
• r_plus_minus is intended for risk-aware controller policy: a narrower sigma can justify higher coupling, while a wider sigma should force conservative action bounds.

Practical overview

Bayesian UPDE is the uncertainty surface for domains where a single trajectory is not enough to support operational action.

The module keeps one strict boundary: uncertainty must remain explicit in the audit record. r_plus_minus is not a display-only value; it is intended to feed risk-aware decision logic and conservative policy envelopes.

That is why backend-gated execution is important here. The code path refuses to promote incomplete uncertainty backends into a production claim while preserving the deterministic NumPy baseline as an auditable anchor.

How teams usually use this surface

• Start from deterministic fitting (fit_gaussian_upde_posterior) on observed trajectories.
• Pass distributions through bayesian_upde_run.
• Compare posterior spread against control limits before accepting aggressive knob proposals.

bayesian

Uncertainty propagation for Kuramoto UPDE rollouts.

The shipped backend is deterministic NumPy Monte Carlo over explicit distributions for omega and K_nm. Probabilistic-programming backends are reserved as fail-closed names until their samplers are implemented and benchmarked against this reproducible baseline.

Classes

GaussianArrayDistribution dataclass

GaussianArrayDistribution(
    mean: object,
    std: object,
    non_negative: bool = False,
    zero_diagonal: bool = False,
)

Independent Gaussian array distribution with optional matrix guards.

Attributes

shape property

shape: tuple[int, ...]

Return the event shape sampled by this Gaussian distribution.

Returns

tuple[int, ...] Return the event shape sampled by this Gaussian distribution.

Methods:

sample

sample(rng: Generator, n_samples: int) -> FloatArray

Draw finite Gaussian samples with optional support guards applied.

Parameters

rng : np.random.Generator NumPy random generator used for sampling. n_samples : int Number of samples to draw.

Returns

FloatArray Finite Gaussian samples, shape (n_samples, *event_shape).

Raises

TypeError If rng is not a NumPy random generator.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def sample(self, rng: np.random.Generator, n_samples: int) -> FloatArray:
    """Draw finite Gaussian samples with optional support guards applied.

    Parameters
    ----------
    rng : np.random.Generator
        NumPy random generator used for sampling.
    n_samples : int
        Number of samples to draw.

    Returns
    -------
    FloatArray
        Finite Gaussian samples, shape ``(n_samples, *event_shape)``.

    Raises
    ------
    TypeError
        If ``rng`` is not a NumPy random generator.
    """
    if not isinstance(rng, np.random.Generator):
        raise TypeError("rng must be a numpy.random.Generator")
    sample_count = _validate_positive_integer(
        n_samples,
        name="n_samples",
        minimum=1,
    )
    draws = rng.normal(
        loc=np.asarray(self.mean, dtype=np.float64),
        scale=np.asarray(self.std, dtype=np.float64),
        size=(sample_count, *self.shape),
    )
    samples = np.asarray(draws, dtype=np.float64)
    if self.non_negative:
        samples = np.maximum(samples, 0.0)
    if self.zero_diagonal:
        diag = np.arange(self.shape[0])
        samples[:, diag, diag] = 0.0
    return np.ascontiguousarray(samples, dtype=np.float64)

BayesianUPDEConfig dataclass

BayesianUPDEConfig(
    n_samples: int = 128,
    seed: int | None = None,
    dt: float = 0.01,
    n_steps: int = 1,
    method: MethodName = "rk4",
    credible_interval: float = 0.95,
    backend: BackendName = "numpy",
    n_substeps: int = 1,
    atol: float = 1e-06,
    rtol: float = 0.001,
)

Configuration for Bayesian UPDE uncertainty propagation.

BayesianUPDEResult dataclass

BayesianUPDEResult(
    r_samples: FloatArray,
    final_phase_samples: FloatArray,
    omega_mean: FloatArray,
    knm_mean: FloatArray,
    r_mean: float,
    r_sigma: float,
    r_lower: float,
    r_upper: float,
    psi_mean: float,
    sample_count: int,
    credible_interval: float,
    backend: str,
    method: str,
)

Posterior predictive order-parameter summary.

Attributes

r_plus_minus property

r_plus_minus: tuple[float, float]

Return the compact R ± sigma pair.

Returns

tuple[float, float] Return the compact R ± sigma pair.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe uncertainty diagnostics.

Returns

dict[str, object] Return JSON-safe uncertainty diagnostics.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe uncertainty diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe uncertainty diagnostics.
    """
    return {
        "kind": "bayesian_upde",
        "backend": self.backend,
        "method": self.method,
        "sample_count": self.sample_count,
        "credible_interval": self.credible_interval,
        "r_summary": {
            "mean": self.r_mean,
            "sigma": self.r_sigma,
            "lower": self.r_lower,
            "upper": self.r_upper,
            "plus_minus": [self.r_mean, self.r_sigma],
        },
        "psi_mean": self.psi_mean,
        "omega_mean": self.omega_mean.tolist(),
        "knm_mean": self.knm_mean.tolist(),
        "final_phase_mean": np.mean(self.final_phase_samples, axis=0).tolist(),
        "diagnostics": {
            "finite_samples": bool(
                np.all(np.isfinite(self.r_samples))
                and np.all(np.isfinite(self.final_phase_samples))
            ),
            "r_min": float(np.min(self.r_samples)),
            "r_max": float(np.max(self.r_samples)),
        },
    }

BayesianBackendStatus dataclass

BayesianBackendStatus(
    backend: str,
    available: bool,
    fail_closed: bool,
    reason: str,
    sample_count: int,
)

Execution status for one Bayesian UPDE backend name.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe backend availability diagnostics.

Returns

dict[str, object] Return JSON-safe backend availability diagnostics.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe backend availability diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe backend availability diagnostics.
    """
    return {
        "kind": "bayesian_backend_status",
        "backend": self.backend,
        "available": self.available,
        "fail_closed": self.fail_closed,
        "reason": self.reason,
        "sample_count": self.sample_count,
    }

GaussianUPDEPosteriorFit dataclass

GaussianUPDEPosteriorFit(
    omega: GaussianArrayDistribution,
    knm: GaussianArrayDistribution,
    residual_rmse: float,
    sample_count: int,
    dt: float,
    ridge: float,
    backend: str = "numpy_lstsq",
)

Gaussian posterior approximation fitted from observed phase trajectories.

Methods:

to_audit_record

to_audit_record() -> dict[str, object]

Return JSON-safe posterior-fit diagnostics.

Returns

dict[str, object] Return JSON-safe posterior-fit diagnostics.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def to_audit_record(self) -> dict[str, object]:
    """Return JSON-safe posterior-fit diagnostics.

    Returns
    -------
    dict[str, object]
        Return JSON-safe posterior-fit diagnostics.
    """
    return {
        "kind": "gaussian_upde_posterior_fit",
        "backend": self.backend,
        "sample_count": self.sample_count,
        "dt": self.dt,
        "ridge": self.ridge,
        "residual_rmse": self.residual_rmse,
        "omega": {
            "mean": np.asarray(self.omega.mean, dtype=np.float64).tolist(),
            "std": np.asarray(self.omega.std, dtype=np.float64).tolist(),
        },
        "knm": {
            "mean": np.asarray(self.knm.mean, dtype=np.float64).tolist(),
            "std": np.asarray(self.knm.std, dtype=np.float64).tolist(),
        },
        "diagnostics": {
            "finite": bool(
                np.all(np.isfinite(np.asarray(self.omega.mean)))
                and np.all(np.isfinite(np.asarray(self.omega.std)))
                and np.all(np.isfinite(np.asarray(self.knm.mean)))
                and np.all(np.isfinite(np.asarray(self.knm.std)))
                and np.isfinite(self.residual_rmse)
            ),
            "zero_diagonal": bool(
                np.allclose(np.diag(np.asarray(self.knm.mean)), 0.0)
                and np.allclose(np.diag(np.asarray(self.knm.std)), 0.0)
            ),
        },
    }

Functions:

fit_gaussian_upde_posterior

fit_gaussian_upde_posterior(
    phase_trajectory: object,
    *,
    dt: float,
    alpha: object | None = None,
    ridge: float = 1e-06,
    coupling_std_floor: float = 1e-06,
    omega_std_floor: float = 1e-06,
) -> GaussianUPDEPosteriorFit

Fit Gaussian omega and K_nm priors from observed phases.

The estimator is a deterministic NumPy ridge least-squares baseline. It fits the Kuramoto right-hand side independently per target oscillator:

d theta_i / dt = omega_i + sum_j K_ij sin(theta_j - theta_i - alpha_ij).

The result is intentionally review-only: it produces distributions that can feed :func:bayesian_upde_run, but it does not apply control actions.

Parameters

phase_trajectory : object Observed phase trajectory, shape (T, N). dt : float Integration step size. alpha : object | None Phase-lag matrix in radians, shape (N, N), or None for no lag. ridge : float Ridge-regularisation strength for the prior fit. coupling_std_floor : float Lower bound on the fitted coupling standard deviation. omega_std_floor : float Lower bound on the fitted natural-frequency standard deviation.

Returns

GaussianUPDEPosteriorFit The fitted Gaussian omega and K_nm prior.

Raises

ValueError If the phase trajectory is empty or non-finite.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def fit_gaussian_upde_posterior(
    phase_trajectory: object,
    *,
    dt: float,
    alpha: object | None = None,
    ridge: float = 1e-6,
    coupling_std_floor: float = 1e-6,
    omega_std_floor: float = 1e-6,
) -> GaussianUPDEPosteriorFit:
    """Fit Gaussian ``omega`` and ``K_nm`` priors from observed phases.

    The estimator is a deterministic NumPy ridge least-squares baseline. It
    fits the Kuramoto right-hand side independently per target oscillator:

    ``d theta_i / dt = omega_i + sum_j K_ij sin(theta_j - theta_i - alpha_ij)``.

    The result is intentionally review-only: it produces distributions that can
    feed :func:`bayesian_upde_run`, but it does not apply control actions.

    Parameters
    ----------
    phase_trajectory : object
        Observed phase trajectory, shape ``(T, N)``.
    dt : float
        Integration step size.
    alpha : object | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    ridge : float
        Ridge-regularisation strength for the prior fit.
    coupling_std_floor : float
        Lower bound on the fitted coupling standard deviation.
    omega_std_floor : float
        Lower bound on the fitted natural-frequency standard deviation.

    Returns
    -------
    GaussianUPDEPosteriorFit
        The fitted Gaussian ``omega`` and ``K_nm`` prior.

    Raises
    ------
    ValueError
        If the phase trajectory is empty or non-finite.
    """
    trajectory = _as_finite_array(phase_trajectory, name="phase_trajectory")
    if trajectory.ndim != 2:
        raise ValueError(
            f"phase_trajectory must be a 2-D array, got shape {trajectory.shape}"
        )
    if trajectory.shape[0] < 3:
        raise ValueError("phase_trajectory must contain at least three samples")
    dt_value = _validate_positive_finite(dt, name="dt")
    ridge_value = _validate_non_negative_finite(ridge, name="ridge")
    coupling_floor = _validate_non_negative_finite(
        coupling_std_floor,
        name="coupling_std_floor",
    )
    omega_floor = _validate_non_negative_finite(
        omega_std_floor,
        name="omega_std_floor",
    )
    n_samples, n_oscillators = trajectory.shape
    alpha_array = (
        np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
        if alpha is None
        else _as_finite_array(alpha, name="alpha")
    )
    if alpha_array.shape != (n_oscillators, n_oscillators):
        raise ValueError(
            f"alpha must have shape {(n_oscillators, n_oscillators)}, "
            f"got {alpha_array.shape}"
        )

    unwrapped = np.unwrap(trajectory, axis=0)
    derivatives = np.diff(unwrapped, axis=0) / dt_value
    theta = trajectory[:-1]
    omega_mean = np.empty(n_oscillators, dtype=np.float64)
    omega_std = np.empty(n_oscillators, dtype=np.float64)
    knm_mean = np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
    knm_std = np.zeros((n_oscillators, n_oscillators), dtype=np.float64)
    residuals: list[FloatArray] = []

    for target in range(n_oscillators):
        source_indices = [source for source in range(n_oscillators) if source != target]
        features = np.column_stack(
            [
                np.ones(theta.shape[0], dtype=np.float64),
                *[
                    np.sin(
                        theta[:, source]
                        - theta[:, target]
                        - alpha_array[target, source]
                    )
                    for source in source_indices
                ],
            ]
        )
        target_derivative = derivatives[:, target]
        gram = features.T @ features
        penalty = ridge_value * np.eye(gram.shape[0], dtype=np.float64)
        penalty[0, 0] = 0.0
        coeffs = np.linalg.solve(gram + penalty, features.T @ target_derivative)
        predicted = features @ coeffs
        residual = target_derivative - predicted
        residuals.append(residual)
        dof = max(1, features.shape[0] - features.shape[1])
        residual_sigma = float(np.sqrt(float(residual @ residual) / dof))
        covariance = residual_sigma**2 * np.linalg.pinv(gram + penalty)
        coefficient_std = np.sqrt(np.maximum(np.diag(covariance), 0.0))
        omega_mean[target] = coeffs[0]
        omega_std[target] = max(float(coefficient_std[0]), omega_floor)
        for offset, source in enumerate(source_indices, start=1):
            value = max(float(coeffs[offset]), 0.0)
            knm_mean[target, source] = value
            knm_std[target, source] = (
                max(float(coefficient_std[offset]), coupling_floor)
                if value > 0.0
                else 0.0
            )

    residual_vector = np.concatenate(residuals)
    residual_rmse = float(np.sqrt(float(np.mean(residual_vector**2))))
    return GaussianUPDEPosteriorFit(
        omega=GaussianArrayDistribution(omega_mean, omega_std),
        knm=GaussianArrayDistribution(
            knm_mean,
            knm_std,
            non_negative=True,
            zero_diagonal=True,
        ),
        residual_rmse=residual_rmse,
        sample_count=n_samples,
        dt=dt_value,
        ridge=ridge_value,
    )

audit_bayesian_backend_status

audit_bayesian_backend_status(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
    backends: tuple[BackendName, ...] = (
        "numpy",
        "numpyro",
        "blackjax",
    ),
) -> tuple[BayesianBackendStatus, ...]

Probe Bayesian backend names without silently accepting unsupported ones.

Parameters

phases : object Oscillator phases in radians, shape (N,). omega : object Natural-frequency distribution or array. knm : object Coupling matrix K_nm, shape (N, N). alpha : object Phase-lag matrix in radians, shape (N, N). Use a zero matrix for no lag. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. config : BayesianUPDEConfig | None Optional configuration object, or None for defaults. backends : tuple[BackendName, ...] Backend names to probe, in priority order.

Returns

tuple[BayesianBackendStatus, ...] Per-backend availability diagnostics.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def audit_bayesian_backend_status(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
    backends: tuple[BackendName, ...] = ("numpy", "numpyro", "blackjax"),
) -> tuple[BayesianBackendStatus, ...]:
    """Probe Bayesian backend names without silently accepting unsupported ones.

    Parameters
    ----------
    phases : object
        Oscillator phases in radians, shape ``(N,)``.
    omega : object
        Natural-frequency distribution or array.
    knm : object
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : object
        Phase-lag matrix in radians, shape ``(N, N)``. Use a zero matrix for no
        lag.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    config : BayesianUPDEConfig | None
        Optional configuration object, or ``None`` for defaults.
    backends : tuple[BackendName, ...]
        Backend names to probe, in priority order.

    Returns
    -------
    tuple[BayesianBackendStatus, ...]
        Per-backend availability diagnostics.
    """
    base_config = config or BayesianUPDEConfig(n_samples=8, seed=0, n_steps=1)
    statuses: list[BayesianBackendStatus] = []
    for backend in backends:
        backend_config = replace(base_config, backend=backend)
        try:
            result = bayesian_upde_run(
                phases,
                omega=omega,
                knm=knm,
                alpha=alpha,
                zeta=zeta,
                psi=psi,
                config=backend_config,
            )
        except NotImplementedError as exc:
            statuses.append(
                BayesianBackendStatus(
                    backend=backend,
                    available=False,
                    fail_closed=True,
                    reason=str(exc),
                    sample_count=0,
                )
            )
        else:
            statuses.append(
                BayesianBackendStatus(
                    backend=backend,
                    available=True,
                    fail_closed=False,
                    reason="executed",
                    sample_count=result.sample_count,
                )
            )
    return tuple(statuses)

bayesian_upde_run

bayesian_upde_run(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
) -> BayesianUPDEResult

Run UPDE over sampled omega and K_nm distributions.

Parameters

phases : object Oscillator phases in radians, shape (N,). omega : object Natural-frequency distribution or array. knm : object Coupling matrix K_nm, shape (N, N). alpha : object Phase-lag matrix in radians, shape (N, N). Use a zero matrix for no lag. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. config : BayesianUPDEConfig | None Optional configuration object, or None for defaults.

Returns

BayesianUPDEResult The Bayesian UPDE result with uncertainty diagnostics.

Raises

NotImplementedError If the requested backend is not implemented. ValueError If the sampled inputs are invalid.

Source code in src/scpn_phase_orchestrator/upde/bayesian.py

def bayesian_upde_run(
    phases: object,
    *,
    omega: object,
    knm: object,
    alpha: object,
    zeta: float,
    psi: float,
    config: BayesianUPDEConfig | None = None,
) -> BayesianUPDEResult:
    """Run UPDE over sampled ``omega`` and ``K_nm`` distributions.

    Parameters
    ----------
    phases : object
        Oscillator phases in radians, shape ``(N,)``.
    omega : object
        Natural-frequency distribution or array.
    knm : object
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    alpha : object
        Phase-lag matrix in radians, shape ``(N, N)``. Use a zero matrix for no
        lag.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    config : BayesianUPDEConfig | None
        Optional configuration object, or ``None`` for defaults.

    Returns
    -------
    BayesianUPDEResult
        The Bayesian UPDE result with uncertainty diagnostics.

    Raises
    ------
    NotImplementedError
        If the requested backend is not implemented.
    ValueError
        If the sampled inputs are invalid.
    """
    resolved = config or BayesianUPDEConfig()
    if resolved.backend != "numpy":
        raise NotImplementedError(
            f"{resolved.backend} Bayesian UPDE backend is not implemented; "
            "use backend='numpy' for reproducible Monte Carlo propagation"
        )

    phase_array = _as_finite_array(phases, name="phases")
    if phase_array.ndim != 1:
        raise ValueError(f"phases must be 1-D, got shape {phase_array.shape}")
    n_oscillators = int(phase_array.shape[0])
    alpha_array = _as_finite_array(alpha, name="alpha")
    if alpha_array.shape != (n_oscillators, n_oscillators):
        raise ValueError(
            f"alpha must have shape {(n_oscillators, n_oscillators)}, "
            f"got {alpha_array.shape}"
        )
    zeta_value = _validate_finite_real(zeta, name="zeta")
    psi_value = _validate_finite_real(psi, name="psi")

    rng = np.random.default_rng(resolved.seed)
    omega_samples = _sample_array(
        omega,
        expected_shape=(n_oscillators,),
        n_samples=resolved.n_samples,
        rng=rng,
        name="omega",
    )
    knm_samples = _sample_array(
        knm,
        expected_shape=(n_oscillators, n_oscillators),
        n_samples=resolved.n_samples,
        rng=rng,
        name="knm",
    )

    final_phase_samples: FloatArray = np.empty(
        (resolved.n_samples, n_oscillators), dtype=np.float64
    )
    r_samples: FloatArray = np.empty(resolved.n_samples, dtype=np.float64)
    psi_samples: FloatArray = np.empty(resolved.n_samples, dtype=np.float64)
    for idx in range(resolved.n_samples):
        final_phases = upde_run(
            phase_array,
            omega_samples[idx],
            knm_samples[idx],
            alpha_array,
            zeta_value,
            psi_value,
            resolved.dt,
            resolved.n_steps,
            method=resolved.method,
            n_substeps=resolved.n_substeps,
            atol=resolved.atol,
            rtol=resolved.rtol,
        )
        final_phase_samples[idx] = final_phases
        order_r, order_psi = compute_order_parameter(final_phases)
        r_samples[idx] = order_r
        psi_samples[idx] = order_psi

    tail = (1.0 - resolved.credible_interval) / 2.0
    lower, upper = np.quantile(r_samples, [tail, 1.0 - tail])
    return BayesianUPDEResult(
        r_samples=np.ascontiguousarray(r_samples, dtype=np.float64),
        final_phase_samples=np.ascontiguousarray(
            final_phase_samples,
            dtype=np.float64,
        ),
        omega_mean=np.mean(omega_samples, axis=0),
        knm_mean=np.mean(knm_samples, axis=0),
        r_mean=float(np.mean(r_samples)),
        r_sigma=float(np.std(r_samples, ddof=1)),
        r_lower=float(lower),
        r_upper=float(upper),
        psi_mean=float(np.mean(psi_samples)),
        sample_count=resolved.n_samples,
        credible_interval=resolved.credible_interval,
        backend=resolved.backend,
        method=resolved.method,
    )