Significance module

gp_methods

Gaussian Process-based significance testing methods.

estimate_correlation_length(chain, nbins=20, k_start=0.001, k_end=0.34, min_correlation=0.1, binning=None)

Estimate correlation length from empirical bin-to-bin correlations.

Fits exponential decay model to the empirical correlation matrix: \(\mathrm{corr}(\Delta \log k) \sim \exp(-|\Delta \log k| / \ell)\)

This provides a data-driven estimate of the correlation length \(\ell\) in \(\log(k)\) space by analyzing how correlation decays with separation between bins.

PARAMETER	DESCRIPTION
chain	MCMC chain object containing delta_i parameters.
nbins	Number of bins in power spectrum reconstruction (default: 20) TYPE: int DEFAULT: 20
k_start	Minimum wavenumber in \(\mathrm{Mpc}^{-1}\) (default: 1e-3) TYPE: float DEFAULT: 0.001
k_end	Maximum wavenumber in \(\mathrm{Mpc}^{-1}\) (default: 0.34) TYPE: float DEFAULT: 0.34
min_correlation	Minimum correlation threshold to include in exponential fit. Only bin pairs with correlation > min_correlation are used. TYPE: float DEFAULT: 0.1
binning	optional BinningScheme object (supersedes nbins, k_start, k_end if provided) TYPE: Optional[BinningScheme] DEFAULT: None

RETURNS	DESCRIPTION
float	Estimated correlation length \(\ell\) in \(\log(k)\) units. Returns 0.5 as
float	default if fewer than 5 bin pairs meet the correlation threshold.

Source code in src/primefeat/significance/gp_methods.py

def estimate_correlation_length(
    chain,
    nbins: int = 20,
    k_start: float = 1e-3,
    k_end: float = 0.34,
    min_correlation: float = 0.1,
    binning: Optional[BinningScheme] = None,
) -> float:
    """
    Estimate correlation length from empirical bin-to-bin correlations.

    Fits exponential decay model to the empirical correlation matrix:
    $\\mathrm{corr}(\\Delta \\log k) \\sim \\exp(-|\\Delta \\log k| / \\ell)$

    This provides a data-driven estimate of the correlation length $\\ell$ in
    $\\log(k)$ space by analyzing how correlation decays with separation between bins.

    Args:
        chain: MCMC chain object containing delta_i parameters.
        nbins: Number of bins in power spectrum reconstruction (default: 20)
        k_start: Minimum wavenumber in $\\mathrm{Mpc}^{-1}$ (default: 1e-3)
        k_end: Maximum wavenumber in $\\mathrm{Mpc}^{-1}$ (default: 0.34)
        min_correlation: Minimum correlation threshold to include in exponential fit.
            Only bin pairs with correlation > min_correlation are used.
        binning: optional BinningScheme object (supersedes nbins, k_start, k_end if provided)

    Returns:
        Estimated correlation length $\\ell$ in $\\log(k)$ units. Returns 0.5 as
        default if fewer than 5 bin pairs meet the correlation threshold.
    """
    # Resolve binning scheme
    b = _resolve_binning(binning, k_start, k_end, nbins)

    log_k = b.log_bin_centers

    # Extract delta samples and compute correlation matrix
    delta_samples = extract_delta_samples(chain, b)
    corr_matrix = compute_delta_correlation(chain, b)

    # Extract distance-correlation pairs
    distances = []
    correlations = []
    for i in range(b.nbins):
        for j in range(i + 1, b.nbins):
            dist = abs(log_k[j] - log_k[i])
            corr = corr_matrix[i, j]
            if corr > min_correlation:  # Only use positive correlations
                distances.append(dist)
                correlations.append(corr)

    if len(distances) < 5:
        print(
            f"Warning: Only {len(distances)} bin pairs with correlation > {min_correlation}"
        )
        print("Using default length scale = 0.5")
        return 0.5

    distances = np.array(distances)
    correlations = np.array(correlations)

    # Fit exponential decay: corr = exp(-dist / length_scale)
    # Taking log: log(corr) = -dist / length_scale
    # So: length_scale = -mean(dist) / mean(log(corr))

    log_corr = np.log(correlations)
    length_scale = -np.mean(distances) / np.mean(log_corr)

    print(f"Estimated correlation length scale: {length_scale:.3f} in log(k)")
    print(f"  Based on {len(distances)} bin pairs with corr > {min_correlation}")

    return length_scale

gp_significance_test(chain, nbins=20, k_start=0.001, k_end=0.34, alpha=0.05, method='null', length_scale=None, sigma=1.0, use_full_covariance=True, noise_level=0.01, compute_landscape=False, landscape_kwargs=None, kernel_type=KernelType.RBF, kernel_params=None, sigma_range=None, length_scale_range=None, n_sigma=100, n_length=100, backend='numpy', max_samples=None, sample_selection='stride', random_state=None, n_bootstrap=200, pre_optimize_extra_params=True, grid_config=None, binning=None)

Test significance using correlation-aware Gaussian Process methods.

This function implements three complementary significance testing approaches for primordial power spectrum features, each with different null hypotheses and use cases.

Method: 'null' (default) - Statistically rigorous feature detection

Null hypothesis: $H_0: \delta \sim \mathcal{N}(0, \Sigma_{\mathrm{post}})$

Test statistic: $\chi^2 = \delta^T \Sigma_{\mathrm{post}}^{-1} \delta \sim \chi^2(\mathrm{nbins})$

Use for: Initial feature detection. No circular reasoning from estimating
correlation lengths. Tests only against posterior covariance.

Method: 'lml' - Bayesian model comparison for characterization

Null: $H_0: \delta \sim \mathcal{N}(0, \Sigma_{\mathrm{post}})$

Alternative: $H_1: \delta \sim \mathcal{N}(0, \sigma^2 K(\ell) + \Sigma_{\mathrm{post}})$

Use for: Characterizing detected features. Optimizes signal amplitude
$\sigma$ and correlation length $\ell$ via log-marginal likelihood.
Computes Bayes factors $\mathrm{BF} = p(D|H_1) / p(D|H_0)$ for evidence
quantification.

Method: 'null+GP' - Null test with GP-fitted covariance

Null hypothesis: $H_0: \delta \sim \mathcal{N}(0, \sigma^2 K(\ell) + \Sigma_{\mathrm{post}})$

Test statistic: $\chi^2 = \delta^T K_{\mathrm{total}}^{-1} \delta \sim \chi^2(\mathrm{nbins})$

Workflow:
    1. Fit GP hyperparameters $(\sigma^*, \ell^*)$ on $\bar{\delta}$ via LML grid search
    2. Construct $K_{\mathrm{total}} = \sigma^{*2} K_{\mathrm{GP}}(\ell^*) + \Sigma_{\mathrm{post}}$
    3. Run null test using $K_{\mathrm{total}}$ as covariance
    4. Return same structure as 'null' for compatibility

Use for: Testing with signal covariance included in null hypothesis.

PARAMETER	DESCRIPTION
chain	MCMC chain object containing delta_i parameters.
nbins	Number of bins in power spectrum reconstruction (default: 20) TYPE: int DEFAULT: 20
k_start	Minimum wavenumber in \(\mathrm{Mpc}^{-1}\) (default: 1e-3) TYPE: float DEFAULT: 0.001
k_end	Maximum wavenumber in \(\mathrm{Mpc}^{-1}\) (default: 0.34) TYPE: float DEFAULT: 0.34
alpha	Significance level \(\alpha\) for hypothesis testing. TYPE: float DEFAULT: 0.05
binning	optional BinningScheme object (supersedes nbins, k_start, k_end if provided) TYPE: Optional[BinningScheme] DEFAULT: None
method	Testing method - 'null', 'lml', or 'null+GP'. TYPE: str DEFAULT: 'null'
kernel_type	GP kernel type. Can be KernelType enum or string: - KernelType.RBF or 'rbf': Squared exponential kernel (default) - KernelType.MATERN or 'matern': Matern kernel - KernelType.RATIONAL_QUADRATIC or 'rq': Multi-scale features - KernelType.PERIODIC or 'periodic': Oscillatory features - KernelType.LOCALLY_PERIODIC or 'locally_periodic': Decaying periodic TYPE: Union[str, KernelType] DEFAULT: RBF
kernel_params	Kernel-specific parameters: - Matern: {'nu': 0.5, 1.5, or 2.5} - smoothness parameter (default: 1.5) - Rational quadratic: {'alpha': float} - scale mixture parameter - Periodic: {'period': float} - oscillation period in \(\log(k)\) space - Locally periodic: {'period': float, 'length_scale_rbf': float} TYPE: Optional[Dict] DEFAULT: None
sigma_range	Range for signal amplitude grid \((\sigma_{\mathrm{min}}, \sigma_{\mathrm{max}})\). Used in 'lml' and 'null+GP' methods. Default: auto-scaled from data. TYPE: Optional[Tuple[float, float]] DEFAULT: None
length_scale_range	Range for correlation length grid \((\ell_{\mathrm{min}}, \ell_{\mathrm{max}})\). Used in 'lml' and 'null+GP' methods. Default: based on \(\log(k)\) range. TYPE: Optional[Tuple[float, float]] DEFAULT: None
n_sigma	Number of grid points for \(\sigma\) in LML optimization. TYPE: int DEFAULT: 100
n_length	Number of grid points for \(\ell\) in LML optimization. TYPE: int DEFAULT: 100
backend	Computation backend: - 'numpy': Standard NumPy implementation - 'jax': GPU-accelerated implementation (50-100x faster for large grids) TYPE: str DEFAULT: 'numpy'
max_samples	Maximum number of posterior samples to use. If provided and smaller than chain size, downselects samples before covariance/test computations to reduce memory usage. TYPE: Optional[int] DEFAULT: None
sample_selection	Downselection strategy when max_samples is set: - 'stride': Evenly spaced deterministic subsample (default) - 'random': Random subsample without replacement TYPE: str DEFAULT: 'stride'
random_state	Random seed used when sample_selection='random'. TYPE: Optional[int] DEFAULT: None
length_scale	Fixed correlation length \(\ell\) (deprecated, for backward compatibility). TYPE: Optional[float] DEFAULT: None
sigma	Fixed signal amplitude \(\sigma\) (deprecated, for backward compatibility). TYPE: float DEFAULT: 1.0
use_full_covariance	Whether to use \(\Sigma_{\mathrm{post}}\) (deprecated, always True). TYPE: bool DEFAULT: True
noise_level	Diagonal noise regularization (deprecated). TYPE: float DEFAULT: 0.01
compute_landscape	If True, force method='lml' (backward compatibility). TYPE: bool DEFAULT: False
landscape_kwargs	Additional kwargs for landscape computation (deprecated). TYPE: Optional[Dict] DEFAULT: None

RETURNS

DESCRIPTION

GPSignificanceResult

GPSignificanceResult containing: - test_statistics: \(\chi^2\) values for each posterior sample - p_values: P-values under \(H_0\) for each sample - global_p_value: Median p-value across samples - fraction_significant: Fraction of samples with p < \(\alpha\) - bin_contributions: Per-bin contribution to \(\chi^2\) statistic - length_scale: Fitted \(\ell\) (None for 'null' method) - covariance_matrix: Covariance \(K\) used in null hypothesis - lml_landscape: LML grid and Bayes factors (for 'lml' and 'null+GP' methods)

Examples:

Standard null test (recommended for initial detection):

>>> result = gp_significance_test(chain, method='null')
>>> print(f"Global p-value: {result.global_p_value:.4f}")
>>> print(f"Significant: {result.fraction_significant:.1%}")

Bayesian model comparison for characterization:

>>> result = gp_significance_test(chain, method='lml')
>>> print(f"Optimal sigma: {result.lml_landscape['optimal_sigma']:.4f}")
>>> print(f"Log Bayes factor: {result.lml_landscape['log_bayes_factor']:.2f}")

Null test with GP-fitted covariance:

>>> result = gp_significance_test(chain, method='null+GP')
>>> print(f"Fitted length scale: {result.length_scale:.3f}")
>>> print(f"Global p-value: {result.global_p_value:.4f}")

GPU-accelerated computation with JAX backend:

>>> result = gp_significance_test(chain, method='lml', backend='jax')

Source code in src/primefeat/significance/gp_methods.py

def gp_significance_test(
    chain,
    nbins: int = 20,
    k_start: float = 1e-3,
    k_end: float = 0.34,
    alpha: float = 0.05,
    method: str = "null",
    length_scale: Optional[float] = None,
    sigma: float = 1.0,
    use_full_covariance: bool = True,
    noise_level: float = 0.01,
    compute_landscape: bool = False,
    landscape_kwargs: Optional[Dict] = None,
    kernel_type: Union[str, KernelType] = KernelType.RBF,
    kernel_params: Optional[Dict] = None,
    sigma_range: Optional[Tuple[float, float]] = None,
    length_scale_range: Optional[Tuple[float, float]] = None,
    n_sigma: int = 100,
    n_length: int = 100,
    backend: str = "numpy",
    max_samples: Optional[int] = None,
    sample_selection: str = "stride",
    random_state: Optional[int] = None,
    n_bootstrap: int = 200,
    pre_optimize_extra_params: bool = True,
    grid_config: Optional[GridSearchConfig] = None,
    binning: Optional[BinningScheme] = None,
) -> GPSignificanceResult:
    """
    Test significance using correlation-aware Gaussian Process methods.

    This function implements three complementary significance testing approaches
    for primordial power spectrum features, each with different null hypotheses
    and use cases.

    **Method: 'null' (default)** - Statistically rigorous feature detection

        Null hypothesis: $H_0: \\delta \\sim \\mathcal{N}(0, \\Sigma_{\\mathrm{post}})$

        Test statistic: $\\chi^2 = \\delta^T \\Sigma_{\\mathrm{post}}^{-1} \\delta \\sim \\chi^2(\\mathrm{nbins})$

        Use for: Initial feature detection. No circular reasoning from estimating
        correlation lengths. Tests only against posterior covariance.

    **Method: 'lml'** - Bayesian model comparison for characterization

        Null: $H_0: \\delta \\sim \\mathcal{N}(0, \\Sigma_{\\mathrm{post}})$

        Alternative: $H_1: \\delta \\sim \\mathcal{N}(0, \\sigma^2 K(\\ell) + \\Sigma_{\\mathrm{post}})$

        Use for: Characterizing detected features. Optimizes signal amplitude
        $\\sigma$ and correlation length $\\ell$ via log-marginal likelihood.
        Computes Bayes factors $\\mathrm{BF} = p(D|H_1) / p(D|H_0)$ for evidence
        quantification.

    **Method: 'null+GP'** - Null test with GP-fitted covariance

        Null hypothesis: $H_0: \\delta \\sim \\mathcal{N}(0, \\sigma^2 K(\\ell) + \\Sigma_{\\mathrm{post}})$

        Test statistic: $\\chi^2 = \\delta^T K_{\\mathrm{total}}^{-1} \\delta \\sim \\chi^2(\\mathrm{nbins})$

        Workflow:
            1. Fit GP hyperparameters $(\\sigma^*, \\ell^*)$ on $\\bar{\\delta}$ via LML grid search
            2. Construct $K_{\\mathrm{total}} = \\sigma^{*2} K_{\\mathrm{GP}}(\\ell^*) + \\Sigma_{\\mathrm{post}}$
            3. Run null test using $K_{\\mathrm{total}}$ as covariance
            4. Return same structure as 'null' for compatibility

        Use for: Testing with signal covariance included in null hypothesis.

    Args:
        chain: MCMC chain object containing delta_i parameters.
        nbins: Number of bins in power spectrum reconstruction (default: 20)
        k_start: Minimum wavenumber in $\\mathrm{Mpc}^{-1}$ (default: 1e-3)
        k_end: Maximum wavenumber in $\\mathrm{Mpc}^{-1}$ (default: 0.34)
        alpha: Significance level $\\alpha$ for hypothesis testing.
        binning: optional BinningScheme object (supersedes nbins, k_start, k_end if provided)
        method: Testing method - 'null', 'lml', or 'null+GP'.
        kernel_type: GP kernel type. Can be KernelType enum or string:
            - KernelType.RBF or 'rbf': Squared exponential kernel (default)
            - KernelType.MATERN or 'matern': Matern kernel
            - KernelType.RATIONAL_QUADRATIC or 'rq': Multi-scale features
            - KernelType.PERIODIC or 'periodic': Oscillatory features
            - KernelType.LOCALLY_PERIODIC or 'locally_periodic': Decaying periodic
        kernel_params: Kernel-specific parameters:
            - Matern: {'nu': 0.5, 1.5, or 2.5} - smoothness parameter (default: 1.5)
            - Rational quadratic: {'alpha': float} - scale mixture parameter
            - Periodic: {'period': float} - oscillation period in $\\log(k)$ space
            - Locally periodic: {'period': float, 'length_scale_rbf': float}
        sigma_range: Range for signal amplitude grid $(\\sigma_{\\mathrm{min}}, \\sigma_{\\mathrm{max}})$.
            Used in 'lml' and 'null+GP' methods. Default: auto-scaled from data.
        length_scale_range: Range for correlation length grid $(\\ell_{\\mathrm{min}}, \\ell_{\\mathrm{max}})$.
            Used in 'lml' and 'null+GP' methods. Default: based on $\\log(k)$ range.
        n_sigma: Number of grid points for $\\sigma$ in LML optimization.
        n_length: Number of grid points for $\\ell$ in LML optimization.
        backend: Computation backend:
            - 'numpy': Standard NumPy implementation
            - 'jax': GPU-accelerated implementation (50-100x faster for large grids)
        max_samples: Maximum number of posterior samples to use. If provided and
            smaller than chain size, downselects samples before covariance/test
            computations to reduce memory usage.
        sample_selection: Downselection strategy when max_samples is set:
            - 'stride': Evenly spaced deterministic subsample (default)
            - 'random': Random subsample without replacement
        random_state: Random seed used when sample_selection='random'.
        length_scale: Fixed correlation length $\\ell$ (deprecated, for backward compatibility).
        sigma: Fixed signal amplitude $\\sigma$ (deprecated, for backward compatibility).
        use_full_covariance: Whether to use $\\Sigma_{\\mathrm{post}}$ (deprecated, always True).
        noise_level: Diagonal noise regularization (deprecated).
        compute_landscape: If True, force method='lml' (backward compatibility).
        landscape_kwargs: Additional kwargs for landscape computation (deprecated).

    Returns:
        GPSignificanceResult containing:
            - test_statistics: $\\chi^2$ values for each posterior sample
            - p_values: P-values under $H_0$ for each sample
            - global_p_value: Median p-value across samples
            - fraction_significant: Fraction of samples with p < $\\alpha$
            - bin_contributions: Per-bin contribution to $\\chi^2$ statistic
            - length_scale: Fitted $\\ell$ (None for 'null' method)
            - covariance_matrix: Covariance $K$ used in null hypothesis
            - lml_landscape: LML grid and Bayes factors (for 'lml' and 'null+GP' methods)

    Examples:
        Standard null test (recommended for initial detection):

        >>> result = gp_significance_test(chain, method='null')
        >>> print(f"Global p-value: {result.global_p_value:.4f}")
        >>> print(f"Significant: {result.fraction_significant:.1%}")

        Bayesian model comparison for characterization:

        >>> result = gp_significance_test(chain, method='lml')
        >>> print(f"Optimal sigma: {result.lml_landscape['optimal_sigma']:.4f}")
        >>> print(f"Log Bayes factor: {result.lml_landscape['log_bayes_factor']:.2f}")

        Null test with GP-fitted covariance:

        >>> result = gp_significance_test(chain, method='null+GP')
        >>> print(f"Fitted length scale: {result.length_scale:.3f}")
        >>> print(f"Global p-value: {result.global_p_value:.4f}")

        GPU-accelerated computation with JAX backend:

        >>> result = gp_significance_test(chain, method='lml', backend='jax')
    """
    # Resolve binning scheme
    b = _resolve_binning(binning, k_start, k_end, nbins)

    if grid_config is not None:
        if sigma_range is None:
            sigma_range = grid_config.sigma_range
        if length_scale_range is None:
            length_scale_range = grid_config.length_scale_range
        if n_sigma == 100:
            n_sigma = grid_config.n_sigma
        if n_length == 100:
            n_length = grid_config.n_length

    # Handle backward compatibility
    if compute_landscape:
        method = "lml"

    log_k = b.log_bin_centers.reshape(-1, 1)

    if sample_selection not in {"stride", "random"}:
        raise ValueError(
            f"sample_selection must be 'stride' or 'random', got '{sample_selection}'"
        )

    n_chain_samples = len(chain[b.bin_param_names[0]])
    sample_indices = None

    if max_samples is not None:
        if max_samples <= 0:
            raise ValueError(f"max_samples must be positive, got {max_samples}")

        if max_samples < n_chain_samples:
            if sample_selection == "stride":
                sample_indices = np.linspace(
                    0, n_chain_samples - 1, num=max_samples, dtype=int
                )
            else:
                rng = np.random.default_rng(random_state)
                sample_indices = np.sort(
                    rng.choice(n_chain_samples, size=max_samples, replace=False)
                )

    # Extract posterior samples (optionally downselected)
    if sample_indices is None:
        delta_samples = extract_delta_samples(chain, b)
    else:
        delta_samples = np.array(
            [chain[name][sample_indices] for name in b.bin_param_names]
        ).T

    n_samples = len(delta_samples)

    # Compute posterior covariance (always needed)
    Sigma_post = np.cov(delta_samples.T)

    _console.print(
        f"[bold]GP Significance Test[/bold] | method=[cyan]{method}[/cyan] | "
        f"samples=[green]{n_samples}[/green] | bins=[green]{b.nbins}[/green] | "
        f"k=[{b.k_start}, {b.k_end}] Mpc⁻¹ | backend=[cyan]{backend}[/cyan]"
    )

    if sample_indices is not None:
        _console.print(
            f"[yellow]Downselected posterior samples:[/yellow] "
            f"{n_samples}/{n_chain_samples} using [cyan]{sample_selection}[/cyan]"
        )

    if method == "null":
        # =================================================================
        # null TEST: Test against Σ_post only (statistically correct)
        # =================================================================
        _console.print("[dim]H₀: δ ~ N(0, Σ_post)[/dim]")

        if backend == "jax":
            # JAX backend: vectorized computation (50-100x faster)
            try:
                from ..backends.jax import (
                    run_null_test_jax,
                    is_jax_significance_available,
                )

                if not is_jax_significance_available():
                    raise ImportError("JAX not available")

                test_statistics, contributions_matrix = run_null_test_jax(
                    delta_samples, Sigma_post, nbins
                )
                # contributions_matrix is (nbins, n_samples)

                # Build bin_contributions list from JAX output
                bin_contributions = []
                for i in range(nbins):
                    contributions = contributions_matrix[i]
                    bin_contributions.append(
                        {
                            "bin_index": i + 1,
                            "k_center": b.bin_centers[i],
                            "mean_contribution": float(np.mean(contributions)),
                            "std_contribution": float(np.std(contributions)),
                            "percentile_95": float(np.percentile(contributions, 95)),
                        }
                    )

            except ImportError as e:
                _console.print(
                    f"[yellow]Warning:[/yellow] JAX unavailable ({e}), using NumPy. "
                    "Please check that jax/tinygp are installed in the active environment."
                )
                backend = "numpy"

        if backend == "numpy":
            # NumPy backend: sequential computation
            # Regularize if needed
            try:
                Sigma_inv_factor = cho_factor(Sigma_post)
            except np.linalg.LinAlgError:
                _console.print(
                    "[yellow]Warning:[/yellow] Σ_post singular, adding regularization"
                )
                Sigma_post = Sigma_post + 1e-6 * np.eye(b.nbins)
                Sigma_inv_factor = cho_factor(Sigma_post)

            # Compute test statistic: χ² = δᵀ Σ_post⁻¹ δ
            test_statistics = []
            for delta in delta_samples:
                chi2_val = delta @ cho_solve(Sigma_inv_factor, delta)
                test_statistics.append(chi2_val)

            test_statistics = np.array(test_statistics)

            # Compute bin contributions: c_i = delta_i * [Sigma_post^{-1} delta]_i
            # Exact decomposition: sum_i(c_i) = T^2 for every sample
            Sigma_inv_deltas = cho_solve(
                Sigma_inv_factor, delta_samples.T
            ).T  # (n_samples, nbins)
            per_sample_contributions = (
                delta_samples * Sigma_inv_deltas
            )  # (n_samples, nbins)

            bin_contributions = []
            for i in range(b.nbins):
                contributions = per_sample_contributions[:, i]
                bin_contributions.append(
                    {
                        "bin_index": i + 1,
                        "k_center": b.bin_centers[i],
                        "mean_contribution": np.mean(contributions),
                        "std_contribution": np.std(contributions),
                        "percentile_95": np.percentile(contributions, 95),
                    }
                )

        # Compute p-values from test statistics
        p_values = 1 - chi2.cdf(test_statistics, df=nbins)
        global_p_value = np.median(p_values)
        fraction_significant = np.mean(p_values < alpha)

        # Build summary panel
        if global_p_value < alpha:
            verdict = " Features detected " + f"at {alpha}"
            verdict_icon = "✓"
            verdict_style = "bold green"
        else:
            verdict = "[dim]No significant features[/dim]" + f" at {alpha}"
            verdict_icon = "✗"
            verdict_style = "bold red"

        # Color p-value based on significance
        if global_p_value < 0.01:
            p_color = "green"
        elif global_p_value < 0.05:
            p_color = "yellow"
        else:
            p_color = "red"

        summary_text = Text()
        summary_text.append("Global p-value: ", style="bold")
        summary_text.append(f"{global_p_value:.4f}", style=p_color)
        summary_text.append("  |  Fraction significant: ", style="bold")
        summary_text.append(f"{fraction_significant:.1%}", style="cyan")
        summary_text.append("\n" + f"{verdict_icon} {verdict}", style=verdict_style)

        _console.print(Panel(summary_text, title="Results", border_style="blue"))

        _print_bin_contributions_table(bin_contributions)

        return GPSignificanceResult(
            test_statistics=test_statistics,
            p_values=p_values,
            fraction_significant=fraction_significant,
            bin_contributions=bin_contributions,
            length_scale=None,
            covariance_matrix=Sigma_post,
            global_p_value=global_p_value,
            log_k=log_k,
            delta_values=delta_samples.mean(axis=0),
            lml_landscape=None,
        )

    elif method == "lml":
        # =================================================================
        # LML LANDSCAPE: Bayesian model comparison for characterization
        # =================================================================
        _console.print("[dim]H₀ vs H₁: σ² K(ℓ) + Σ_post[/dim]")

        delta_mean = delta_samples.mean(axis=0)

        # Set defaults for ranges
        if sigma_range is None:
            data_std = np.std(delta_mean)
            sigma_range = (0.01 * data_std, 5.0 * data_std)

        if length_scale_range is None:
            log_k_range = np.log(b.k_end / b.k_start)
            length_scale_range = (log_k_range / 4, log_k_range)

        sigma_grid = np.linspace(sigma_range[0], sigma_range[1], n_sigma)
        length_scale_grid = np.linspace(
            length_scale_range[0], length_scale_range[1], n_length
        )

        # Normalize kernel type to enum
        kernel_type_enum = _normalize_kernel_type(kernel_type)
        _console.print(
            f"[dim]Grid: {n_sigma}×{n_length} | kernel={kernel_type_enum.value}[/dim]"
        )

        # Pre-optimize extra hyperparameters for complex kernels
        if pre_optimize_extra_params and _kernel_has_extra_hyperparams(
            kernel_type_enum
        ):
            _console.print("[dim]Pre-optimising extra hyperparameters...[/dim]")
            kernel_params = _pre_optimize_extra_params(
                kernel_type_enum,
                kernel_params or {},
                delta_mean,
                log_k.ravel(),
                Sigma_post,
                sigma_init=float(np.mean(sigma_range)),
                length_scale_init=float(np.mean(length_scale_range)),
                verbose=False,
            )
            _console.print(f"[dim]Extra params → {kernel_params}[/dim]")

        if backend == "jax":
            # JAX backend: vectorized grid evaluation (50-100x faster)
            try:
                from ..backends.jax import (
                    run_lml_grid_jax,
                    is_jax_significance_available,
                )

                if not is_jax_significance_available():
                    raise ImportError("JAX not available")

                lml_grid = run_lml_grid_jax(
                    delta_mean,
                    log_k.ravel(),
                    Sigma_post,
                    sigma_grid,
                    length_scale_grid,
                    kernel_type=kernel_type_enum,
                    kernel_params=kernel_params,
                )
            except ImportError as e:
                _console.print(
                    f"[yellow]Warning:[/yellow] JAX unavailable ({e}), using NumPy. "
                    "Please check that jax/tinygp are installed in the active environment."
                )
                backend = "numpy"

        null_lml, K_null_factor, log_det_K_null = _compute_null_lml(
            Sigma_post, delta_mean, b.nbins
        )

        if backend == "numpy":
            lml_grid = _compute_lml_grid_numpy(
                delta_mean,
                log_k,
                Sigma_post,
                sigma_grid,
                length_scale_grid,
                kernel_type_enum,
                kernel_params,
                b.nbins,
            )

        optimal_sigma, optimal_length_scale, optimal_lml = _find_grid_optimum(
            lml_grid, sigma_grid, length_scale_grid
        )
        log_BF = optimal_lml - null_lml
        BF = np.exp(log_BF)

        # Profile likelihood CIs for sigma and ell (free — grid already computed)
        # Profile: maximize LML over the other parameter
        # Thresholds from Wilks' theorem for 1 parameter:
        #   1sigma (68%): delta_LML > -0.5
        #   2sigma (95%): delta_LML > -2.0
        profile_sigma = lml_grid.max(axis=1)  # max over ell for each sigma
        profile_ell = lml_grid.max(axis=0)  # max over sigma for each ell

        sigma_1sigma = sigma_grid[profile_sigma > optimal_lml - 0.5]
        sigma_2sigma = sigma_grid[profile_sigma > optimal_lml - 2.0]
        ell_1sigma = length_scale_grid[profile_ell > optimal_lml - 0.5]
        ell_2sigma = length_scale_grid[profile_ell > optimal_lml - 2.0]

        sigma_ci_1sigma = (
            (float(sigma_1sigma.min()), float(sigma_1sigma.max()))
            if len(sigma_1sigma)
            else (optimal_sigma, optimal_sigma)
        )
        sigma_ci_2sigma = (
            (float(sigma_2sigma.min()), float(sigma_2sigma.max()))
            if len(sigma_2sigma)
            else (optimal_sigma, optimal_sigma)
        )
        ell_ci_1sigma = (
            (float(ell_1sigma.min()), float(ell_1sigma.max()))
            if len(ell_1sigma)
            else (optimal_length_scale, optimal_length_scale)
        )
        ell_ci_2sigma = (
            (float(ell_2sigma.min()), float(ell_2sigma.max()))
            if len(ell_2sigma)
            else (optimal_length_scale, optimal_length_scale)
        )

        # BIC and AIC penalty-corrected model comparison
        # n = nbins (dimension of data vector delta), k = 2 free parameters (sigma, ell)
        # Delta = metric(H1) - metric(H0); negative values favor H1
        delta_BIC = -2 * log_BF + 2 * np.log(nbins)
        delta_AIC = -2 * log_BF + 4

        # Interpret Bayes factor (Kass & Raftery 1995 scale)
        if log_BF < 0:
            bf_interp = "Favors null (no signal)"
            bf_color = "red"
        elif log_BF < 1:
            bf_interp = "Weak evidence"
            bf_color = "dim"
        elif log_BF < 3:
            bf_interp = "Positive evidence"
            bf_color = "yellow"
        elif log_BF < 5:
            bf_interp = "Strong evidence"
            bf_color = "green"
        else:
            bf_interp = "Very strong evidence"
            bf_color = "bold green"

        # Build results table
        table = Table(
            title="LML Landscape Results",
            box=box.ROUNDED,
            show_header=True,
            header_style="bold cyan",
        )
        table.add_column("Parameter", style="bold")
        table.add_column("Value", justify="right")
        table.add_column("Description", style="dim")

        table.add_row("σ (amplitude)", f"{optimal_sigma:.4f}", "Signal std deviation")
        table.add_row(
            "σ 1σ CI",
            f"[{sigma_ci_1sigma[0]:.4f}, {sigma_ci_1sigma[1]:.4f}]",
            "Profile likelihood",
        )
        table.add_row(
            "ℓ (length scale)", f"{optimal_length_scale:.3f}", "In log(k) space"
        )
        table.add_row(
            "ℓ 1σ CI",
            f"[{ell_ci_1sigma[0]:.3f}, {ell_ci_1sigma[1]:.3f}]",
            "Profile likelihood",
        )
        table.add_row("Max LML", f"{optimal_lml:.2f}", "Log-marginal likelihood")
        table.add_row("Null LML", f"{null_lml:.2f}", "σ → 0 limit")
        table.add_row(
            "log(BF)",
            f"[{bf_color}]{log_BF:.2f}[/{bf_color}]",
            f"[{bf_color}]{bf_interp}[/{bf_color}]",
        )
        table.add_row("Bayes Factor", f"[{bf_color}]{BF:.2e}[/{bf_color}]", "H1 / H0")
        bic_color = "green" if delta_BIC < 0 else "red"
        aic_color = "green" if delta_AIC < 0 else "red"
        table.add_row(
            r"$\Delta$BIC",
            f"[{bic_color}]{delta_BIC:.2f}[/{bic_color}]",
            "< 0 favors H1",
        )
        table.add_row(
            r"$\Delta$AIC",
            f"[{aic_color}]{delta_AIC:.2f}[/{aic_color}]",
            "< 0 favors H1",
        )

        # Build optimal K
        K_signal_opt = _build_signal_kernel_matrix(
            log_k, optimal_sigma, optimal_length_scale, kernel_type_enum, kernel_params
        )
        K_opt = K_signal_opt + Sigma_post

        # Bootstrap credible interval on ln B
        # Resample posterior mean to capture uncertainty from finite chain length.
        # Hyperparameters (sigma*, ell*) are fixed at the values found above.
        rng = np.random.default_rng(random_state)
        K_opt_factor = cho_factor(K_opt)
        log_det_K_opt = 2 * np.sum(np.log(np.diag(K_opt_factor[0])))
        log_BF_bootstrap = []
        for _ in range(n_bootstrap):
            idx = rng.choice(len(delta_samples), size=len(delta_samples), replace=True)
            delta_boot = delta_samples[idx].mean(axis=0)
            null_lml_boot = -0.5 * (
                delta_boot @ cho_solve(K_null_factor, delta_boot)
                + log_det_K_null
                + nbins * np.log(2 * np.pi)
            )
            alt_lml_boot = -0.5 * (
                delta_boot @ cho_solve(K_opt_factor, delta_boot)
                + log_det_K_opt
                + nbins * np.log(2 * np.pi)
            )
            log_BF_bootstrap.append(alt_lml_boot - null_lml_boot)
        log_BF_bootstrap = np.array(log_BF_bootstrap)
        log_BF_ci = (
            float(np.percentile(log_BF_bootstrap, 16)),
            float(np.percentile(log_BF_bootstrap, 84)),
        )

        table.add_row(
            "ln(B) 68% CI",
            f"[{bf_color}][{log_BF_ci[0]:.2f}, {log_BF_ci[1]:.2f}][/{bf_color}]",
            f"Bootstrap ({n_bootstrap} resamples)",
        )
        _console.print(table)

        # Per-bin decomposition of the data-fit part of ln B
        # ln B = data_fit_term + complexity_penalty
        # data_fit_term = 0.5 * [delta^T K_null^{-1} delta - delta^T K_opt^{-1} delta]
        #               = sum_i c_i^fit,  where c_i^fit = 0.5 * delta_i * [(K_null^{-1} - K_opt^{-1}) delta]_i
        # complexity_penalty = 0.5 * [log|K_null| - log|K_opt|]  (scalar, Occam factor)
        K_null_inv_delta_mean = cho_solve(K_null_factor, delta_mean)
        K_opt_inv_delta_mean = cho_solve(K_opt_factor, delta_mean)
        bin_lml_contributions = (
            0.5 * delta_mean * (K_null_inv_delta_mean - K_opt_inv_delta_mean)
        )
        complexity_penalty = 0.5 * (log_det_K_null - log_det_K_opt)

        landscape = {
            "sigma_grid": sigma_grid,
            "length_scale_grid": length_scale_grid,
            "lml_grid": lml_grid,
            "optimal_sigma": optimal_sigma,
            "optimal_length_scale": optimal_length_scale,
            "optimal_lml": optimal_lml,
            "null_lml": null_lml,
            "log_bayes_factor": log_BF,
            "bayes_factor": BF,
            "delta_BIC": delta_BIC,
            "delta_AIC": delta_AIC,
            "log_bayes_factor_bootstrap": log_BF_bootstrap,
            "log_bayes_factor_ci": log_BF_ci,
            "bin_lml_contributions": bin_lml_contributions,
            "complexity_penalty": complexity_penalty,
            "profile_sigma": profile_sigma,
            "profile_ell": profile_ell,
            "sigma_ci_1sigma": sigma_ci_1sigma,
            "sigma_ci_2sigma": sigma_ci_2sigma,
            "ell_ci_1sigma": ell_ci_1sigma,
            "ell_ci_2sigma": ell_ci_2sigma,
            "log_k": log_k,
            "delta_values": delta_mean,
            "noise_level": np.sqrt(np.mean(np.diag(Sigma_post))),
            "K": K_opt,
            "K_signal": K_signal_opt,
            "K_noise": Sigma_post,
            "optimized_kernel_params": kernel_params,
        }

        return GPSignificanceResult(
            test_statistics=None,
            p_values=None,
            fraction_significant=None,
            bin_contributions=[],
            length_scale=optimal_length_scale,
            covariance_matrix=K_opt,
            global_p_value=None,
            log_k=log_k,
            delta_values=delta_mean,
            lml_landscape=landscape,
            log_bayes_factor_ci=log_BF_ci,
        )

    elif method == "null+GP":
        # =================================================================
        # null+GP: Null test with GP-fitted covariance
        # =================================================================
        # Like 'null' but uses K_total = σ*² K_GP(ℓ*) + Σ_post instead of Σ_post
        # where (σ*, ℓ*) are fitted via LML optimization on the posterior mean.
        # =================================================================
        _console.print("[dim]H₀: δ ~ N(0, σ²K(ℓ) + Σ_post)[/dim]")

        delta_mean = delta_samples.mean(axis=0)

        # -----------------------------------------------------------------
        # Step 1: Run LML grid to find optimal (σ*, ℓ*)
        # -----------------------------------------------------------------
        if sigma_range is None:
            data_std = np.std(delta_mean)
            sigma_range = (0.01 * data_std, 5 * data_std)

        if length_scale_range is None:
            log_k_range = np.log(k_end / k_start)
            length_scale_range = (log_k_range / 4, 3 * log_k_range)

        # sigma_grid = np.linspace(sigma_range[0], sigma_range[1], n_sigma)
        sigma_grid = np.geomspace(sigma_range[0], sigma_range[1], n_sigma)
        # length_scale_grid = np.linspace(
        length_scale_grid = np.geomspace(
            length_scale_range[0], length_scale_range[1], n_length
        )

        # Normalize kernel type to enum
        kernel_type_enum = _normalize_kernel_type(kernel_type)
        _console.print(
            f"[dim]Grid: {n_sigma}×{n_length} | kernel={kernel_type_enum.value}[/dim]"
        )

        # Pre-optimize extra hyperparameters for complex kernels
        if pre_optimize_extra_params and _kernel_has_extra_hyperparams(
            kernel_type_enum
        ):
            _console.print("[dim]Pre-optimising extra hyperparameters...[/dim]")
            kernel_params = _pre_optimize_extra_params(
                kernel_type_enum,
                kernel_params or {},
                delta_mean,
                log_k.ravel(),
                Sigma_post,
                sigma_init=float(np.mean(sigma_range)),
                length_scale_init=float(np.mean(length_scale_range)),
                verbose=False,
            )
            _console.print(f"[dim]Extra params → {kernel_params}[/dim]")

        if backend == "jax":
            try:
                from ..backends.jax import (
                    run_lml_grid_jax,
                    is_jax_significance_available,
                )

                if not is_jax_significance_available():
                    raise ImportError("JAX not available")

                lml_grid = run_lml_grid_jax(
                    delta_mean,
                    log_k.ravel(),
                    Sigma_post,
                    sigma_grid,
                    length_scale_grid,
                    kernel_type=kernel_type_enum,
                    kernel_params=kernel_params,
                )
            except ImportError as e:
                _console.print(
                    f"[yellow]Warning:[/yellow] JAX unavailable ({e}), using NumPy. "
                    "Please check that jax/tinygp are installed in the active environment."
                )
                backend = "numpy"

        null_lml, _, _ = _compute_null_lml(Sigma_post, delta_mean, nbins)

        if backend == "numpy":
            lml_grid = _compute_lml_grid_numpy(
                delta_mean,
                log_k,
                Sigma_post,
                sigma_grid,
                length_scale_grid,
                kernel_type_enum,
                kernel_params,
                b.nbins,
            )

        optimal_sigma, optimal_length_scale, optimal_lml = _find_grid_optimum(
            lml_grid, sigma_grid, length_scale_grid
        )
        log_BF = optimal_lml - null_lml

        # -----------------------------------------------------------------
        # Step 2: Build K_total = σ*² K_GP(ℓ*) + Σ_post
        # -----------------------------------------------------------------
        K_signal_opt = _build_signal_kernel_matrix(
            log_k, optimal_sigma, optimal_length_scale, kernel_type_enum, kernel_params
        )
        K_total = K_signal_opt + Sigma_post

        # -----------------------------------------------------------------
        # Step 3: Run null test using K_total as covariance
        # -----------------------------------------------------------------
        if backend == "jax":
            try:
                from ..backends.jax import (
                    run_null_test_jax,
                    is_jax_significance_available,
                )

                if not is_jax_significance_available():
                    raise ImportError("JAX not available")

                test_statistics, contributions_matrix = run_null_test_jax(
                    delta_samples, K_total, b.nbins
                )

                # Build bin_contributions list from JAX output
                bin_contributions = []
                for i in range(b.nbins):
                    contributions = contributions_matrix[i]
                    bin_contributions.append(
                        {
                            "bin_index": i + 1,
                            "k_center": b.bin_centers[i],
                            "mean_contribution": float(np.mean(contributions)),
                            "std_contribution": float(np.std(contributions)),
                            "percentile_95": float(np.percentile(contributions, 95)),
                        }
                    )

            except ImportError as e:
                _console.print(
                    f"[yellow]Warning:[/yellow] JAX unavailable ({e}), using NumPy. "
                    "Please check that jax/tinygp are installed in the active environment."
                )
                backend = "numpy"

        if backend == "numpy":
            # Regularize if needed
            try:
                K_total_factor = cho_factor(K_total)
            except np.linalg.LinAlgError:
                _console.print(
                    "[yellow]Warning:[/yellow] K_total singular, adding regularization"
                )
                K_total = K_total + 1e-6 * np.eye(b.nbins)
                K_total_factor = cho_factor(K_total)

            # Compute test statistic: χ² = δᵀ K_total⁻¹ δ
            test_statistics = []
            for delta in delta_samples:
                chi2_val = delta @ cho_solve(K_total_factor, delta)
                test_statistics.append(chi2_val)

            test_statistics = np.array(test_statistics)

            # Compute bin contributions: c_i = delta_i * [K_total^{-1} delta]_i
            # Exact decomposition: sum_i(c_i) = T^2 for every sample
            K_total_inv_deltas = cho_solve(
                K_total_factor, delta_samples.T
            ).T  # (n_samples, nbins)
            per_sample_contributions = (
                delta_samples * K_total_inv_deltas
            )  # (n_samples, nbins)

            bin_contributions = []
            for i in range(b.nbins):
                contributions = per_sample_contributions[:, i]
                bin_contributions.append(
                    {
                        "bin_index": i + 1,
                        "k_center": b.bin_centers[i],
                        "mean_contribution": np.mean(contributions),
                        "std_contribution": np.std(contributions),
                        "percentile_95": np.percentile(contributions, 95),
                    }
                )

        # Compute p-values from test statistics
        p_values = 1 - chi2.cdf(test_statistics, df=b.nbins)
        global_p_value = np.median(p_values)
        fraction_significant = np.mean(p_values < alpha)

        # Build summary panel
        if global_p_value < alpha:
            verdict = f"Features detected at α={alpha}"
            verdict_icon = "✓"
        else:
            verdict = f"No significant features at α={alpha}"
            verdict_icon = "✗"

        # Color p-value based on significance
        if global_p_value < 0.01:
            p_color = "green"
        elif global_p_value < 0.05:
            p_color = "yellow"
        else:
            p_color = "red"

        summary_text = Text()
        summary_text.append("Fitted GP: ", style="bold")
        summary_text.append(f"σ = {optimal_sigma:.4f}", style="cyan")
        summary_text.append(", ", style="dim")
        summary_text.append(f"ℓ = {optimal_length_scale:.3f}", style="cyan")
        summary_text.append("\nGlobal p-value: ", style="bold")
        summary_text.append(f"{global_p_value:.4f}", style=p_color)
        summary_text.append("  |  Fraction significant: ", style="bold")
        summary_text.append(f"{fraction_significant:.1%}", style="cyan")
        summary_text.append(
            f"\n{verdict_icon} {verdict}",
            style="green" if global_p_value < alpha else "red",
        )

        _console.print(Panel(summary_text, title="Results", border_style="blue"))

        _print_bin_contributions_table(bin_contributions)

        # Build landscape dict for compatibility
        landscape = {
            "sigma_grid": sigma_grid,
            "length_scale_grid": length_scale_grid,
            "lml_grid": lml_grid,
            "optimal_sigma": optimal_sigma,
            "optimal_length_scale": optimal_length_scale,
            "optimal_lml": optimal_lml,
            "null_lml": null_lml,
            "log_bayes_factor": log_BF,
            "bayes_factor": np.exp(log_BF),
            "log_k": log_k,
            "delta_values": delta_mean,
            "noise_level": np.sqrt(np.mean(np.diag(Sigma_post))),
            "K": K_total,
            "K_signal": K_signal_opt,
            "K_noise": Sigma_post,
            "optimized_kernel_params": kernel_params,
        }

        return GPSignificanceResult(
            test_statistics=test_statistics,
            p_values=p_values,
            fraction_significant=fraction_significant,
            bin_contributions=bin_contributions,
            length_scale=optimal_length_scale,
            covariance_matrix=K_total,
            global_p_value=global_p_value,
            log_k=log_k,
            delta_values=delta_mean,
            lml_landscape=landscape,
        )

    else:
        raise ValueError(f"Unknown method: {method}")