quadsv.statistics ================= .. py:module:: quadsv.statistics Functions --------- .. autoapisummary:: quadsv.statistics.auto_chunk_size quadsv.statistics.compute_null_params quadsv.statistics.effective_rank quadsv.statistics.gene_pattern_diversity quadsv.statistics.liu_sf quadsv.statistics.spatial_q_test quadsv.statistics.spatial_r_test quadsv.statistics.within_group_pattern_diversity Module Contents --------------- .. py:function:: auto_chunk_size(kernel, n_jobs = 1, budget_bytes = _DEFAULT_CHUNK_BUDGET) Pick a per-backend-optimal ``chunk_size`` for the Q / R test. The returned value is used by :func:`spatial_q_test` / :func:`spatial_r_test` (and by :meth:`DetectorGrid.compute_qstat` / :meth:`DetectorIrregular.compute_qstat`) to split a multi-feature batch into chunks. It is the smaller of two caps: 1. **Cache sweet-spot cap** — empirical sweep of per-feature time vs ``chunk`` at ``n ∈ {30k, 100k, 300k, 1M}``: .. list-table:: :header-rows: 1 :widths: 50 50 * - Backend - chunk cap * - :class:`~quadsv.FFTKernel` - 32 * - :class:`~quadsv.NUFFTKernel` - 64 * - MatrixKernel (any sub-type) - 16 (``n < 200k``); 8 (``n ≥ 200k``) Matrix backends don't vectorise over RHS columns (scipy CSR SpMV, SuperLU triangular solve), and the chunk size cap is determined empirically for best per-feature speed under the given memory constraints. FFT / NUFFT *do* benefit from BLAS / ``n_transf`` batching, but their complex workspace spills L3 past the listed cap (15× slowdown for FFT at chunk=512, 1.9× for NUFFT at chunk=256). 2. **Memory cap** — ``budget_bytes / n_jobs // per_feat``, where ``per_feat`` is the backend-specific transient bytes per feature: - MatrixKernel dense / sparse: ``16 · n`` - MatrixKernel precision-stored CAR: ``24 · n`` - FFTKernel: ``24 · n`` - NUFFTKernel: ``16 · ny·nx + 8 · n`` :param kernel: The backend kernel the chunk will operate on. :type kernel: Kernel :param n_jobs: Number of parallel workers the caller plans to use. The ``budget_bytes`` is divided by ``n_jobs`` so aggregate live memory stays bounded. :type n_jobs: int, default 1 :param budget_bytes: Aggregate live-memory cap across *all* workers. :type budget_bytes: int, default 2 GiB :returns: A ``chunk_size`` in ``[8, chunk_cap]``. The lower bound of 8 ensures measurements stay meaningful even when a single feature consumes most of the per-worker budget. :rtype: int .. py:function:: compute_null_params(kernel, method = 'welch', k_eigen = None, dirichlet_correction = True, liu_n_probes = None) Pre-compute null distribution parameters for spatial tests. Call this ONCE before running parallel tests on thousands of features. Caches the expensive computations (traces, cumulants, shifted-χ² fit) for reuse across both Q-tests and R-tests. :param kernel: The spatial kernel object (MatrixKernel, FFTKernel, NUFFTKernel, or compatible). :type kernel: Kernel :param method: Null approximation method for the **Q-test**. The R-test entry ``var_R = trace(K̃²)`` is always populated alongside, regardless of ``method`` — R-tests use a Normal approximation and only need this one moment. - 'clt': Central Limit Theorem (Z-score normal approximation) - 'welch': Welch-Satterthwaite moment matching (fast, uses traces) - 'liu': Liu eigenvalue-based approximation (accurate tail, slower) :type method: {'clt', 'welch', 'liu'}, default 'welch' :param k_eigen: Number of top eigenvalues to compute if method='liu' and kernel is sparse. If None, computes all available eigenvalues. Ignored when ``liu_n_probes`` is set (eigenvalues are bypassed entirely). :type k_eigen: int, optional :param liu_n_probes: If set, bypass the eigendecomposition for ``method='liu'`` and estimate the four spectral cumulants ``c_p = trace(K̃^p)``, ``p = 1..4``, directly from the kernel via Hutchinson probing (:func:`_hutchinson_cumulants`). Cost drops from :math:`O(n^3)` (dense eigensolve) or :math:`O(r^3)` (reduced Toeplitz-M) to :math:`2 \cdot n_\mathrm{probes}` kernel matvecs, at the cost of :math:`O(n_\mathrm{probes}^{-1/2})` Monte-Carlo error in each cumulant. Rule of thumb: ``n_probes = 60`` gives Liu p-values within ``~5\%`` of the eigenvalue-exact baseline; ``n_probes = 120`` within ``~0.2\%``. When ``None`` (default), the eigenvalue path is used. :type liu_n_probes: int, optional :param dirichlet_correction: When ``True``: use the finite-``n`` Dirichlet(1/2) ratio ``Var[Q] = 2 · (m · trace((HKH)²) - trace(HKH)²) / (m+2)`` with ``m = n-1``. When ``False``: drop the ``mean²`` term to the large-``n`` limit ``Var[Q] = 2 · trace((HKH)²)``, a monotonic upper bound that slightly overestimates ``Var[Q]`` at finite ``n``. See **Notes** for the derivation. :type dirichlet_correction: bool, default True :returns: * *dict[str, float or np.ndarray]* -- Always populated (regardless of ``method``): - ``'method'`` : str — the Q-test approximation selected. - ``'var_R'`` : float — ``trace(K̃²)``, the null variance of ``R`` (used by :func:`spatial_r_test`). Method-specific additions: - ``method='liu'`` (default for FFT / NUFFT kernels): * ``'cumulants'`` : ``dict {1: c_1, 2: c_2, 3: c_3, 4: c_4}`` with ``c_p = trace(K̃^p)``. Computed from the full eigendecomposition when available (``liu_n_probes is None``) or from :math:`2m` Hutchinson probes otherwise. * ``'liu_coef'`` : ``dict`` with cached Liu coefficients ``{'mu_Q', 'sigma_Q', 'mu_x', 'sigma_x', 'dof_x', 'delta_x'}`` derived from ``cumulants`` once; consumed by :func:`spatial_q_test` so per-feature p-values reduce to a pure :math:`t`-broadcast. - ``method='welch'`` (default for MatrixKernel Q-tests): * ``'mean_Q'`` : ``trace(K)`` * ``'var_Q'`` : ``2 · trace(K²)`` * ``'scale_g'`` : Welch scale parameter ``var_Q / (2 · mean_Q)`` * ``'df_h'`` : Welch df ``2 · mean_Q² / var_Q`` - ``method='clt'``: ``'mean_Q'``, ``'var_Q'`` only. * Consumers (``spatial_q_test`` / ``spatial_r_test``) read only the keys * *their approximation needs; the dict is safe to reuse across calls.* :raises AssertionError: If method is not one of 'clt', 'welch', 'liu'. .. rubric:: Examples >>> kernel = MatrixKernel.from_coordinates(coords, method='gaussian') >>> params = compute_null_params(kernel, method='welch') >>> Q, pval = spatial_q_test(data, kernel, null_params=params) >>> R, r_pval = spatial_r_test(x, y, kernel, null_params=params) .. rubric:: Notes :func:`spatial_q_test` standardizes its input as :math:`Z = (X - \bar{X}\,\mathbf{1}) / \sigma`, so the realized quadratic form is .. math:: Q \;=\; Z^{\top} K Z \;=\; X^{\top}\, H K H\, X / \sigma^{2}, \qquad H = I - \mathbf{1}\mathbf{1}^{\top} / n. Null moments are therefore for the double-centered operator :math:`\tilde{K} = H K H`, not raw :math:`K`. :math:`Q` is additionally a *ratio* of quadratic forms — the denominator :math:`\sigma^{2} = X^{\top} H X / (n-1)` is a random variable correlated with the numerator. ``dirichlet_correction=True`` applies the finite-:math:`n` correction derived from the Dirichlet(1/2) distribution of :math:`Y_i = X_i^{2} / \sum_j X_j^{2}`: .. math:: \mathrm{Var}[Q] \;=\; \frac{2 \bigl[\, m \cdot \mathrm{tr}(\tilde K^{2}) - \mathrm{tr}(\tilde K)^{2} \,\bigr]}{m + 2}, \qquad m = n - 1. With ``dirichlet_correction=False`` the finite term drops out and :math:`\mathrm{Var}[Q] = 2\,\mathrm{tr}(\tilde K^{2})` (large-:math:`n` limit). .. py:function:: effective_rank(cov, weights = None) Effective rank (participation ratio) of a covariance matrix. Computes .. math:: K_\mathrm{eff} \;=\; \frac{\big(\sum_k \lambda_k\big)^2}{\sum_k \lambda_k^2} where :math:`\lambda_k` are the eigenvalues of ``cov`` (or, when ``weights`` is given, of :math:`W^{1/2} \mathrm{cov}\, W^{1/2}` with :math:`W = \mathrm{diag}(w)`). The result is bounded by 1 (rank-1, all variance on a single direction) and ``K = cov.shape[0]`` (uniformly spread, all eigenvalues equal). It coincides with the standard inverse Herfindahl index of the (normalised) eigenvalue distribution and quantifies the "effective number of independent components" of a quadratic-form statistic :math:`T^2 = X^\top \mathrm{cov} X` where :math:`X \sim \mathcal{N}(0,I)`. :param cov: Symmetric ``(K, K)`` covariance matrix. Negative eigenvalues from numerical noise are clipped to 0. :type cov: np.ndarray :param weights: Non-negative weights of length ``K``. When provided, returns the effective rank of the weighted form :math:`W^{1/2} \mathrm{cov}\,W^{1/2}`, useful for analysing how a frequency-weighted L2 statistic actually distributes its sensitivity across eigen-directions. :type weights: np.ndarray, optional :returns: Effective rank in ``[1, K]``. Returns ``nan`` when the trace is non-positive (degenerate covariance). :rtype: float .. rubric:: Examples >>> import numpy as np >>> effective_rank(np.eye(10)) 10.0 >>> # Rank-1 outer product >>> v = np.zeros(10); v[0] = 1.0 >>> effective_rank(np.outer(v, v)) 1.0 .. py:function:: gene_pattern_diversity(spectra, weights = None, *, eps = 1e-12) Spatial-pattern diversity across genes within a single sample. Quantifies how heterogeneous the per-gene spatial-frequency profiles are within one sample. Computes the cross-gene covariance of the log-spectra, .. math:: \hat\Sigma_\mathrm{genes} \;=\; \frac{1}{G - 1} \sum_g \big(\ell_g - \bar\ell\big)\big(\ell_g - \bar\ell\big)^\top where :math:`\ell_g \in \mathbb{R}^K` is gene :math:`g`'s radially-binned log-spectrum and :math:`\bar\ell` the mean across genes, then returns ``effective_rank(Σ_genes, weights)``. Interpretation: - **Low diversity** (:math:`K_\mathrm{eff} \approx 1`): most genes share the same spatial-frequency profile — the sample's spatial patterns collapse onto a single dominant mode (e.g. all "smooth" or all "punctate"). - **High diversity** (:math:`K_\mathrm{eff} \approx K`): genes vary widely in their spatial structure — the sample carries a rich mix of spatial scales. :param spectra: ``(n_genes, K)`` raw spectrum matrix (typically a single sample's radially-binned spectrum). Log is taken internally with an ``eps`` floor. :type spectra: np.ndarray :param weights: Per-bin weights, same semantics as :func:`effective_rank`. :type weights: np.ndarray, optional :param eps: Floor for ``log(spectra)`` to handle exact-zero bins. :type eps: float, default 1e-12 .. py:function:: liu_sf(t, lambs, dofs = None, deltas = None, kurtosis = False, n = None) Liu approximation to a linear combination of non-central chi-squared variables. One-shot convenience wrapper equivalent to ``_liu_apply(t, _liu_prepare(lambs, ...))``. For multi-feature workloads, prefer the split form: call :func:`_liu_prepare` once on the spectrum and :func:`_liu_apply` for each Q-batch (:func:`compute_null_params` already caches the coefficients under ``null_params['liu_coef']``, so :func:`spatial_q_test` does this automatically). :param t: Test statistic value(s). Array input is broadcast efficiently through a single :func:`scipy.stats.ncx2.sf` call. :type t: float or np.ndarray :param lambs: Eigenvalues of the kernel matrix, shape ``(n_evals,)``. :type lambs: np.ndarray :param dofs: Per-eigenvalue degrees of freedom. Default: ones (chi-squared). :type dofs: np.ndarray, optional :param deltas: Non-centrality parameters. Default: zeros (central). :type deltas: np.ndarray, optional :param kurtosis: If True, use the kurtosis-based edge-case approximation. :type kurtosis: bool, default False :param n: Sample size. When provided, applies the Dirichlet(1/2) variance correction ``Var[Q] = 2·(m·c_2 - c_1²)/(m+2)`` with ``m = n-1`` for the z-scored ratio ``Q = XᵀK̃X/σ̂²``. Essential for broad-spectrum PSD kernels (CAR, graph_laplacian) on dense grids, where the large-:math:`n` limit ``2·c_2`` overestimates ``Var[Q]`` and collapses the tail to zero. Default ``None`` keeps the original large-:math:`n` behavior for back-compat with callers supplying a raw eigenvalue mixture. :type n: int, optional :returns: Tail probability ``Pr(Q > t)`` with the same shape as ``t``. :rtype: float or np.ndarray .. py:function:: spatial_q_test(Xn, kernel, null_params = None, return_pval = True, is_standardized = False, chunk_size = 'auto', show_progress = False) Univariate spatial Q-test for detecting spatial variability. Top-level chunking wrapper — splits the feature batch along the trailing axis into blocks of ``chunk_size`` features, dispatches each block to the backend-specific per-chunk helper (:func:`quadsv.kernels.fft._q_test_fft`, :func:`quadsv.kernels.nufft._q_test_nufft`, or :func:`_q_test_matrix`), and concatenates the results. The per-chunk helpers do **not** handle chunking themselves. :param Xn: Input data of shape ``(n,)`` / ``(n, M)`` for MatrixKernel and NUFFTKernel, or ``(ny, nx)`` / ``(ny, nx, M)`` for FFTKernel. Can be dense numpy array or sparse matrix (CSC/CSR recommended) for MatrixKernel; FFT/NUFFT paths require dense input. :type Xn: np.ndarray or scipy.sparse matrix :param kernel: Pre-constructed :class:`~quadsv.kernels.Kernel` (``MatrixKernel`` / ``FFTKernel`` / ``NUFFTKernel``) or a raw dense / sparse kernel matrix. :type kernel: Kernel :param null_params: Pre-computed null distribution parameters from :func:`compute_null_params`. Resolved once at the top level if ``None`` and shared across chunks (no redundant recomputation). :type null_params: dict, optional :param return_pval: If True, returns ``(Q, pval)``; else returns ``Q`` only. :type return_pval: bool, default True :param is_standardized: If True, skips Z-score standardization internally. :type is_standardized: bool, default False :param chunk_size: Number of features processed per per-chunk dispatch call. ``'auto'`` defers to :func:`auto_chunk_size` (backend-specific cache sweet spot under a 2 GiB live-memory budget); ``-1`` processes the full batch in a single call. For a cross-backend cost model see :doc:`/guides/scaling`. :type chunk_size: int or ``'auto'``, default ``'auto'`` :param show_progress: If True, displays a tqdm bar over chunks (only when ``M > chunk_size``). :type show_progress: bool, default False :returns: * **Q** (*float or np.ndarray*) -- Test statistic value(s). Shape ``(M,)`` for 2-D / 3-D inputs, scalar for 1-D. * **pval** (*float or np.ndarray, optional*) -- Tail probability under null hypothesis; returned only if ``return_pval=True``. Same shape as Q. .. rubric:: Notes Under H₀: data is spatially independent. Under H₁: mean shift present. The test statistic ``Q = xᵀ K x`` approximates a chi-squared mixture under the null; see :doc:`/guides/theory` and :doc:`/guides/scaling`. .. rubric:: Examples >>> coords = np.random.randn(100, 2) >>> kernel = MatrixKernel.from_coordinates(coords, method='gaussian') >>> data = np.random.randn(100) >>> Q, pval = spatial_q_test(data, kernel) >>> # Sparse-matrix batch of features (auto-chunked): >>> from scipy.sparse import csr_matrix >>> sparse_data = csr_matrix(np.random.randn(100, 1000)) >>> Q, pval = spatial_q_test(sparse_data, kernel, show_progress=True) .. py:function:: spatial_r_test(Xn, Yn, kernel, null_params = None, return_pval = True, is_standardized = False, chunk_size = 'auto', show_progress = False) Bivariate spatial R-test for correlation between two spatial variables. Top-level chunking wrapper — splits the paired feature batch along the trailing axis into blocks of ``chunk_size`` features, dispatches each block to the backend-specific per-chunk helper (:func:`quadsv.kernels.fft._r_test_fft`, :func:`quadsv.kernels.nufft._r_test_nufft`, or :func:`_r_test_matrix`), and concatenates the results. The per-chunk helpers do **not** handle chunking themselves. :param Xn: First input. Shape ``(n,)`` / ``(n, M)`` for MatrixKernel and NUFFTKernel, ``(ny, nx)`` / ``(ny, nx, M)`` for FFTKernel. :type Xn: np.ndarray or scipy.sparse matrix :param Yn: Second input, same shape as ``Xn`` (paired R-test). For ``NUFFTKernel`` a bipartite mode with ``M_x != M_y`` is passed through without chunking. :type Yn: np.ndarray or scipy.sparse matrix :param kernel: Pre-constructed :class:`~quadsv.kernels.Kernel`. :type kernel: Kernel :param null_params: Pre-computed null parameters; only ``'var_R'`` is consumed. Resolved once at the top level if ``None`` and shared across chunks. :type null_params: dict, optional :param return_pval: If True, returns ``(R, pval)``; else returns ``R`` only. :type return_pval: bool, default True :param is_standardized: If True, skips Z-score standardization internally. :type is_standardized: bool, default False :param chunk_size: Number of feature pairs processed per per-chunk dispatch call. ``'auto'`` defers to :func:`auto_chunk_size`; ``-1`` processes the full batch in a single call. For a cross-backend cost model see :doc:`/guides/scaling`. :type chunk_size: int or ``'auto'``, default ``'auto'`` :param show_progress: If True, displays a tqdm bar over chunks (only when ``M > chunk_size``). :type show_progress: bool, default False :returns: * **R** (*float or np.ndarray*) -- Test statistic value(s). Shape ``(M,)`` for 2-D / 3-D inputs, scalar for 1-D. For bipartite NUFFT input (``M_x != M_y``), shape ``(M_x, M_y)``. * **pval** (*float or np.ndarray, optional*) -- Two-tailed tail probability under null hypothesis; returned only if ``return_pval=True``. .. rubric:: Notes Under H₀, ``R = xᵀ K y`` is approximated as :math:`\mathcal{N}(0, \mathrm{trace}((HKH)^2))` on z-scored inputs. See :doc:`/guides/theory` and :doc:`/guides/scaling`. .. rubric:: Examples >>> coords = np.random.randn(100, 2) >>> kernel = MatrixKernel.from_coordinates(coords, method='gaussian') >>> x_data = np.random.randn(100) >>> y_data = np.random.randn(100) >>> R, pval = spatial_r_test(x_data, y_data, kernel) .. py:function:: within_group_pattern_diversity(spectra, groups, weights = None, *, eps = 1e-12) Spatial-pattern diversity of the within-group residual covariance. For a cohort of samples partitioned into two groups, computes the pooled-across-genes within-group covariance of log-spectra (the same estimator used by ``log_l2 + null='wald'`` in the comparator), then returns its effective rank. Interpretation: - **Low diversity** (:math:`K_\mathrm{eff} \approx 1`): within-group sample-to-sample variation aligns with one spatial-frequency direction. Wald-type tests on this cohort effectively reduce to a 1-DoF test → high power per direction but very sensitive to estimation noise in that single eigenvalue. - **High diversity** (:math:`K_\mathrm{eff} \approx K`): noise spreads over many directions; Wald's analytic null is more accurate (CLT smoothing). :param spectra: ``(n_samples, n_genes, K)`` spectrum tensor (raw, not logged). :type spectra: np.ndarray :param groups: ``(n_samples,)`` with exactly two distinct labels. :type groups: np.ndarray :param weights: Per-bin weights, same semantics as :func:`effective_rank`. :type weights: np.ndarray, optional :param eps: Floor for ``log(spectra)``. :type eps: float, default 1e-12 :returns: Effective rank of the pooled within-group covariance. :rtype: float