`neurospatial.stats`¶

stats ¶

Statistical methods.

This module provides general-purpose statistical utilities used across encoding, decoding, and behavior analysis.

Submodules

circular : Circular statistics, Rayleigh test shuffle : Shuffle controls, permutation tests surrogates : Surrogate data generation

Imports

from neurospatial.stats import rayleigh_test, circular_mean from neurospatial.stats import shuffle_time_bins, compute_shuffle_pvalue from neurospatial.stats import generate_poisson_surrogates, generate_jittered_spikes from neurospatial.stats.circular import circular_basis from neurospatial.stats.shuffle import ShuffleTestResult from neurospatial.stats.surrogates import generate_inhomogeneous_poisson_surrogates

Classes¶

CircularBasisResult `dataclass` ¶

CircularBasisResult(sin_component: NDArray[float64], cos_component: NDArray[float64], angles: NDArray[float64])

Result from circular_basis() function.

Contains sin/cos components for use in GLM design matrices, plus the original angles for reference.

Attributes:

Name	Type	Description
`sin_component`	`(ndarray, shape(n_samples))`	Sine of angles: sin(angles).
`cos_component`	`(ndarray, shape(n_samples))`	Cosine of angles: cos(angles).
`angles`	`(ndarray, shape(n_samples))`	Original angles (in radians).

Properties

design_matrix : ndarray, shape (n_samples, 2) Design matrix with columns [sin_component, cos_component].

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_basis
>>> angles = np.linspace(0, 2 * np.pi, 100, endpoint=False)
>>> result = circular_basis(angles)
>>> result.design_matrix.shape
(100, 2)

Attributes¶

design_matrix `property` ¶

design_matrix: NDArray[float64]

Return design matrix with columns [sin_component, cos_component].

Returns:

Type	Description
`(ndarray, shape(n_samples, 2))`	Design matrix suitable for GLM.

ShuffleTestResult `dataclass` ¶

ShuffleTestResult(observed_score: float, null_scores: NDArray[float64], p_value: float, z_score: float, shuffle_type: str, n_shuffles: int)

Result container for shuffle-based significance testing.

A frozen dataclass containing the observed score, null distribution, and computed statistics (p-value, z-score) from a shuffle test.

Parameters:

Name	Type	Description	Default
`observed_score`	`float`	The score computed from the original (non-shuffled) data.	required
`null_scores`	`NDArray[float64]`	Array of scores computed from shuffled data, forming the null distribution.	required
`p_value`	`float`	Monte Carlo p-value with correction: (k + 1) / (n + 1) where k is the count of null scores at least as extreme as observed and n is the number of shuffles.	required
`z_score`	`float`	Standard score: (observed - mean(null)) / std(null). NaN if null has zero variance.	required
`shuffle_type`	`str`	Name of the shuffle method used (e.g., "time_bins", "cell_identity").	required
`n_shuffles`	`int`	Number of shuffles performed.	required

Attributes:

Name	Type	Description
`is_significant`	`bool`	True if p_value < 0.05 (convenience property).

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import ShuffleTestResult

>>> result = ShuffleTestResult(
...     observed_score=5.0,
...     null_scores=np.array([1.0, 2.0, 3.0, 4.0]),
...     p_value=0.2,
...     z_score=1.5,
...     shuffle_type="time_bins",
...     n_shuffles=4,
... )
>>> result.is_significant
False
>>> result.p_value
0.2

See Also

compute_shuffle_pvalue : Compute Monte Carlo p-value from null distribution. compute_shuffle_zscore : Compute z-score from null distribution.

Attributes¶

is_significant `property` ¶

is_significant: bool

Return True if p_value < 0.05.

Returns:

Type	Description
`bool`	True if the result is statistically significant at alpha=0.05.

Functions¶

plot ¶

plot(ax: Axes | None = None, **kwargs) -> Axes

Plot the null distribution with observed score.

Creates a histogram of the null distribution with a vertical line showing the observed score and annotations for p-value and z-score.

Parameters:

Name	Type	Description	Default
`ax`	`Axes \| None`	Axes to plot on. If None, creates a new figure.	`None`
`**kwargs`		Additional keyword arguments passed to matplotlib's hist().	`{}`

Returns:

Type	Description
`Axes`	The axes containing the plot.

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import ShuffleTestResult

>>> result = ShuffleTestResult(
...     observed_score=5.0,
...     null_scores=np.random.default_rng(42).normal(2.0, 1.0, 100),
...     p_value=0.01,
...     z_score=3.0,
...     shuffle_type="time_bins",
...     n_shuffles=100,
... )
>>> ax = result.plot()

Source code in src/neurospatial/stats/shuffle.py

def plot(self, ax: Axes | None = None, **kwargs) -> Axes:
    """Plot the null distribution with observed score.

    Creates a histogram of the null distribution with a vertical line
    showing the observed score and annotations for p-value and z-score.

    Parameters
    ----------
    ax : matplotlib.axes.Axes | None, default=None
        Axes to plot on. If None, creates a new figure.
    **kwargs
        Additional keyword arguments passed to matplotlib's hist().

    Returns
    -------
    matplotlib.axes.Axes
        The axes containing the plot.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import ShuffleTestResult

    >>> result = ShuffleTestResult(
    ...     observed_score=5.0,
    ...     null_scores=np.random.default_rng(42).normal(2.0, 1.0, 100),
    ...     p_value=0.01,
    ...     z_score=3.0,
    ...     shuffle_type="time_bins",
    ...     n_shuffles=100,
    ... )
    >>> ax = result.plot()  # doctest: +SKIP
    """
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots()

    # Default histogram settings
    hist_kwargs = {
        "bins": 30,
        "alpha": 0.7,
        "color": "steelblue",
        "edgecolor": "white",
    }
    hist_kwargs.update(kwargs)

    # Plot null distribution
    ax.hist(self.null_scores, **hist_kwargs)  # type: ignore[arg-type]

    # Add observed score line
    ax.axvline(
        self.observed_score,
        color="red",
        linestyle="--",
        linewidth=2,
        label=f"Observed: {self.observed_score:.3f}",
    )

    # Add annotations
    sig_marker = "*" if self.is_significant else ""
    ax.set_title(
        f"Shuffle Test ({self.shuffle_type})\n"
        f"p = {self.p_value:.4f}{sig_marker}, z = {self.z_score:.2f}"
    )
    ax.set_xlabel("Score")
    ax.set_ylabel("Count")
    ax.legend()

    return ax

Functions¶

circular_basis ¶

circular_basis(angles: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad') -> CircularBasisResult

Compute sine/cosine basis functions for circular variables.

Creates a design matrix for use in GLMs (Generalized Linear Models) where the sin and cos components capture circular modulation. This is the standard approach for including circular predictors in regression models.

Parameters:

Name	Type	Description	Default
`angles`	`(array, shape(n_samples))`	Circular variable. Common use cases: Head direction: Animal's facing direction (typically in degrees) Theta phase: LFP phase at spike times (0 to 2π radians) Running direction: Direction of movement Time of day: Circadian phase	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles. Use `'deg'` for head direction data that is typically recorded in degrees.	`'rad'`

Returns:

Type	Description
`CircularBasisResult`	Result object containing: sin_component: sin(angles) cos_component: cos(angles) angles: original angles (in radians) design_matrix property: (n_samples, 2) array for GLM

See Also

circular_basis_metrics : Compute amplitude/phase from GLM coefficients. is_modulated : Quick significance test for circular modulation. plot_circular_basis_tuning : Visualize fitted tuning curve.

Notes

Why sin/cos basis?

For a circular predictor θ, using sin(θ) and cos(θ) as separate predictors allows the model to capture any phase and amplitude of circular modulation. The fitted coefficients (β_sin, β_cos) can be converted to:

Amplitude: sqrt(β_sin² + β_cos²) — modulation strength
Phase: atan2(β_sin, β_cos) — preferred direction/phase

GLM vs. Binned Tuning Curves:

Binned: Non-parametric, intuitive, good for visualization. Requires many spikes per bin.
GLM (circular_basis): Parametric, statistically rigorous, can include covariates. Better for low spike counts.

Examples:

Head direction tuning (GLM approach):

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_basis
>>> # Head direction at each spike (in degrees)
>>> hd_at_spikes = np.array([45, 50, 48, 52, 180, 185, 170])  # degrees
>>> result = circular_basis(hd_at_spikes, angle_unit="deg")
>>> result.design_matrix.shape
(7, 2)

Theta phase modulation:

>>> # LFP phase at each spike (in radians)
>>> theta_phases = np.random.uniform(0, 2 * np.pi, 100)
>>> result = circular_basis(theta_phases)  # angle_unit='rad' is default
>>> X = result.design_matrix
>>> # Now fit GLM: model.fit(sm.add_constant(X), spike_counts)

Source code in src/neurospatial/stats/circular.py

def circular_basis(
    angles: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
) -> CircularBasisResult:
    """
    Compute sine/cosine basis functions for circular variables.

    Creates a design matrix for use in GLMs (Generalized Linear Models) where
    the sin and cos components capture circular modulation. This is the standard
    approach for including circular predictors in regression models.

    Parameters
    ----------
    angles : array, shape (n_samples,)
        Circular variable. Common use cases:

        - **Head direction**: Animal's facing direction (typically in degrees)
        - **Theta phase**: LFP phase at spike times (0 to 2π radians)
        - **Running direction**: Direction of movement
        - **Time of day**: Circadian phase

    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles. Use ``'deg'`` for head direction data that
        is typically recorded in degrees.

    Returns
    -------
    CircularBasisResult
        Result object containing:

        - sin_component: sin(angles)
        - cos_component: cos(angles)
        - angles: original angles (in radians)
        - design_matrix property: (n_samples, 2) array for GLM

    See Also
    --------
    circular_basis_metrics : Compute amplitude/phase from GLM coefficients.
    is_modulated : Quick significance test for circular modulation.
    plot_circular_basis_tuning : Visualize fitted tuning curve.

    Notes
    -----
    **Why sin/cos basis?**

    For a circular predictor θ, using sin(θ) and cos(θ) as separate
    predictors allows the model to capture any phase and amplitude of circular
    modulation. The fitted coefficients (β_sin, β_cos) can be converted to:

    - Amplitude: sqrt(β_sin² + β_cos²) — modulation strength
    - Phase: atan2(β_sin, β_cos) — preferred direction/phase

    **GLM vs. Binned Tuning Curves**:

    - **Binned**: Non-parametric, intuitive, good for visualization.
      Requires many spikes per bin.
    - **GLM** (``circular_basis``): Parametric, statistically rigorous, can
      include covariates. Better for low spike counts.

    Examples
    --------
    **Head direction tuning (GLM approach)**:

    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_basis
    >>> # Head direction at each spike (in degrees)
    >>> hd_at_spikes = np.array([45, 50, 48, 52, 180, 185, 170])  # degrees
    >>> result = circular_basis(hd_at_spikes, angle_unit="deg")
    >>> result.design_matrix.shape
    (7, 2)

    **Theta phase modulation**:

    >>> # LFP phase at each spike (in radians)
    >>> theta_phases = np.random.uniform(0, 2 * np.pi, 100)
    >>> result = circular_basis(theta_phases)  # angle_unit='rad' is default
    >>> X = result.design_matrix
    >>> # Now fit GLM: model.fit(sm.add_constant(X), spike_counts)
    """
    # Convert to radians if needed
    angles = np.asarray(angles, dtype=np.float64).ravel()
    angles_rad = _to_radians(angles, angle_unit)

    # Compute sin/cos components
    sin_component = np.sin(angles_rad)
    cos_component = np.cos(angles_rad)

    return CircularBasisResult(
        sin_component=sin_component,
        cos_component=cos_component,
        angles=angles_rad,
    )

circular_basis_metrics ¶

circular_basis_metrics(beta_sin: float, beta_cos: float, cov_matrix: NDArray[float64] | None = None) -> tuple[float, float, float | None]

Compute amplitude, phase, and p-value from GLM coefficients.

Given fitted coefficients for sin and cos basis functions, compute: - Amplitude (strength of modulation) - Phase (preferred angle) - P-value (statistical significance via Wald test)

Parameters:

Name	Type	Description	Default
`beta_sin`	`float`	Coefficient for sin(angle) from fitted GLM.	required
`beta_cos`	`float`	Coefficient for cos(angle) from fitted GLM.	required
`cov_matrix`	`(ndarray, shape(2, 2))`	Covariance matrix for [beta_sin, beta_cos] from GLM fit. If provided, computes p-value via Wald test.	`None`

Returns:

Name	Type	Description
`amplitude`	`float`	Modulation amplitude: sqrt(beta_sin^2 + beta_cos^2). Larger values indicate stronger circular modulation.
`phase`	`float`	Preferred phase angle in radians, range [-π, π]. Computed as atan2(beta_sin, beta_cos).
`pvalue`	`float or None`	P-value testing H0: amplitude = 0 (no modulation). None if cov_matrix not provided.

See Also

circular_basis : Create design matrix for GLM.

Notes

Interpretation:

amplitude > 0.3 typically indicates meaningful modulation
pvalue < 0.05 indicates statistically significant modulation
phase tells you the preferred angle (e.g., preferred head direction)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_basis_metrics
>>> # After fitting GLM with circular basis
>>> beta_sin, beta_cos = 0.5, 0.3
>>> cov = np.array([[0.01, 0.001], [0.001, 0.01]])
>>> amplitude, phase, pval = circular_basis_metrics(beta_sin, beta_cos, cov)
>>> print(f"Amplitude: {amplitude:.2f}, Phase: {np.degrees(phase):.1f}°")
Amplitude: 0.58, Phase: 59.0°

Source code in src/neurospatial/stats/circular.py

def circular_basis_metrics(
    beta_sin: float,
    beta_cos: float,
    cov_matrix: NDArray[np.float64] | None = None,
) -> tuple[float, float, float | None]:
    """
    Compute amplitude, phase, and p-value from GLM coefficients.

    Given fitted coefficients for sin and cos basis functions, compute:
    - Amplitude (strength of modulation)
    - Phase (preferred angle)
    - P-value (statistical significance via Wald test)

    Parameters
    ----------
    beta_sin : float
        Coefficient for sin(angle) from fitted GLM.
    beta_cos : float
        Coefficient for cos(angle) from fitted GLM.
    cov_matrix : ndarray, shape (2, 2), optional
        Covariance matrix for [beta_sin, beta_cos] from GLM fit.
        If provided, computes p-value via Wald test.

    Returns
    -------
    amplitude : float
        Modulation amplitude: sqrt(beta_sin^2 + beta_cos^2).
        Larger values indicate stronger circular modulation.
    phase : float
        Preferred phase angle in radians, range [-π, π].
        Computed as atan2(beta_sin, beta_cos).
    pvalue : float or None
        P-value testing H0: amplitude = 0 (no modulation).
        None if cov_matrix not provided.

    See Also
    --------
    circular_basis : Create design matrix for GLM.

    Notes
    -----
    **Interpretation**:

    - amplitude > 0.3 typically indicates meaningful modulation
    - pvalue < 0.05 indicates statistically significant modulation
    - phase tells you the preferred angle (e.g., preferred head direction)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_basis_metrics
    >>> # After fitting GLM with circular basis
    >>> beta_sin, beta_cos = 0.5, 0.3
    >>> cov = np.array([[0.01, 0.001], [0.001, 0.01]])
    >>> amplitude, phase, pval = circular_basis_metrics(beta_sin, beta_cos, cov)
    >>> print(f"Amplitude: {amplitude:.2f}, Phase: {np.degrees(phase):.1f}°")
    Amplitude: 0.58, Phase: 59.0°
    """
    # Compute amplitude: sqrt(beta_sin^2 + beta_cos^2)
    amplitude = float(np.sqrt(beta_sin**2 + beta_cos**2))

    # Compute phase: atan2(beta_sin, beta_cos)
    phase = float(np.arctan2(beta_sin, beta_cos))

    # Compute p-value if covariance provided
    pvalue: float | None = None
    if cov_matrix is not None:
        cov_matrix = np.asarray(cov_matrix, dtype=np.float64)
        if cov_matrix.shape != (2, 2):
            raise ValueError(f"cov_matrix must be shape (2, 2), got {cov_matrix.shape}")
        pvalue = _wald_test_magnitude(beta_sin, beta_cos, cov_matrix)

    return amplitude, phase, pvalue

circular_circular_correlation ¶

circular_circular_correlation(angles1: NDArray[float64], angles2: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad') -> tuple[float, float]

Circular-circular correlation coefficient.

Computes the Fisher & Lee (1983) circular correlation coefficient between two circular variables.

Parameters:

Name	Type	Description	Default
`angles1`	`ndarray of shape (n_samples,)`	First circular variable in radians or degrees.	required
`angles2`	`ndarray of shape (n_samples,)`	Second circular variable in radians or degrees. Must have same length as `angles1`.	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles (applies to both arrays).	`'rad'`

Returns:

Name	Type	Description
`r`	`float`	Circular correlation coefficient in [-1, 1]. - r > 0: positive association (angles increase together) - r < 0: negative association (angles move in opposite directions) - r = 0: no circular correlation
`pval`	`float`	P-value from normal approximation.

Raises:

Type	Description
`ValueError`	If arrays have different lengths or insufficient samples.

See Also

circular_linear_correlation : For circular-linear correlation.

Notes

Formula (Fisher & Lee, 1983; Jammalamadaka & SenGupta, 2001, p. 176):

rho = sum(sin(a1 - mean1) * sin(a2 - mean2)) /
      sqrt(sum(sin(a1 - mean1)^2) * sum(sin(a2 - mean2)^2))

Interpretation:

|r| > 0.3: Moderate association
|r| > 0.5: Strong association

Important Properties:

Symmetric: r(a1, a2) == r(a2, a1)
Invariant to constant offsets: r(a, a + c) = 1.0 for any constant c (because deviations from circular means remain identical)
Anticorrelation requires reflection: r(a, -a) approaches -1.0

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_circular_correlation
>>> angles1 = np.random.default_rng(42).uniform(0, 2 * np.pi, 100)
>>> angles2 = angles1 + 0.1 * np.random.default_rng(42).standard_normal(100)
>>> r, p = circular_circular_correlation(angles1, angles2)
>>> r > 0.8
True

References

Fisher, N.I. & Lee, A.J. (1983). A correlation coefficient for circular data. Biometrika, 70(2), 327-332. Jammalamadaka, S.R. & SenGupta, A. (2001). Topics in Circular Statistics. World Scientific, p. 176.

Source code in src/neurospatial/stats/circular.py

def circular_circular_correlation(
    angles1: NDArray[np.float64],
    angles2: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
) -> tuple[float, float]:
    """
    Circular-circular correlation coefficient.

    Computes the Fisher & Lee (1983) circular correlation coefficient
    between two circular variables.

    Parameters
    ----------
    angles1 : ndarray of shape (n_samples,)
        First circular variable in radians or degrees.
    angles2 : ndarray of shape (n_samples,)
        Second circular variable in radians or degrees. Must have same
        length as ``angles1``.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles (applies to both arrays).

    Returns
    -------
    r : float
        Circular correlation coefficient in [-1, 1].
        - r > 0: positive association (angles increase together)
        - r < 0: negative association (angles move in opposite directions)
        - r = 0: no circular correlation
    pval : float
        P-value from normal approximation.

    Raises
    ------
    ValueError
        If arrays have different lengths or insufficient samples.

    See Also
    --------
    circular_linear_correlation : For circular-linear correlation.

    Notes
    -----
    **Formula** (Fisher & Lee, 1983; Jammalamadaka & SenGupta, 2001, p. 176):

        rho = sum(sin(a1 - mean1) * sin(a2 - mean2)) /
              sqrt(sum(sin(a1 - mean1)^2) * sum(sin(a2 - mean2)^2))

    **Interpretation**:

    - |r| > 0.3: Moderate association
    - |r| > 0.5: Strong association

    **Important Properties**:

    - Symmetric: r(a1, a2) == r(a2, a1)
    - Invariant to constant offsets: r(a, a + c) = 1.0 for any constant c
      (because deviations from circular means remain identical)
    - Anticorrelation requires reflection: r(a, -a) approaches -1.0

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_circular_correlation
    >>> angles1 = np.random.default_rng(42).uniform(0, 2 * np.pi, 100)
    >>> angles2 = angles1 + 0.1 * np.random.default_rng(42).standard_normal(100)
    >>> r, p = circular_circular_correlation(angles1, angles2)
    >>> r > 0.8
    True

    References
    ----------
    Fisher, N.I. & Lee, A.J. (1983). A correlation coefficient for circular data.
        Biometrika, 70(2), 327-332.
    Jammalamadaka, S.R. & SenGupta, A. (2001). Topics in Circular Statistics.
        World Scientific, p. 176.
    """
    from scipy import stats

    # Convert to radians if needed
    angles1 = np.asarray(angles1, dtype=np.float64)
    angles2 = np.asarray(angles2, dtype=np.float64)
    angles1 = _to_radians(angles1, angle_unit)
    angles2 = _to_radians(angles2, angle_unit)

    # Validate paired inputs (handles NaN removal and length check)
    angles1, angles2 = _validate_paired_input(
        angles1, angles2, "angles1", "angles2", min_samples=5
    )

    n = len(angles1)

    # Compute circular means using scipy
    alpha_bar = float(stats.circmean(angles1, high=2 * np.pi, low=0))
    beta_bar = float(stats.circmean(angles2, high=2 * np.pi, low=0))

    # Compute deviations from circular mean
    sin_dev1 = np.sin(angles1 - alpha_bar)
    sin_dev2 = np.sin(angles2 - beta_bar)

    # Compute correlation coefficient (Fisher & Lee formula)
    num = np.sum(sin_dev1 * sin_dev2)
    den = np.sqrt(np.sum(sin_dev1**2) * np.sum(sin_dev2**2))

    if den == 0:
        # Degenerate case: no variation in one or both variables
        warnings.warn(
            "No variation in one or both circular variables. "
            "This may indicate constant angles. Returning r=0.",
            category=RuntimeWarning,
            stacklevel=2,
        )
        return 0.0, 1.0

    rho = num / den

    # Clip to valid range (numerical stability)
    rho = float(np.clip(rho, -1.0, 1.0))

    # Significance test (normal approximation)
    # From Jammalamadaka & SenGupta (2001), p. 177
    l20 = np.mean(sin_dev1**2)
    l02 = np.mean(sin_dev2**2)
    l22 = np.mean(sin_dev1**2 * sin_dev2**2)

    if l22 == 0:
        warnings.warn(
            "Cannot compute p-value: variance product is zero. Returning p=1.0.",
            category=RuntimeWarning,
            stacklevel=2,
        )
        return rho, 1.0

    # Test statistic follows standard normal under null hypothesis
    ts = np.sqrt((n * l20 * l02) / l22) * rho
    pval = float(2 * (1 - stats.norm.cdf(np.abs(ts))))

    # Ensure p-value is in valid range
    pval = float(np.clip(pval, 0.0, 1.0))

    return rho, pval

circular_linear_correlation ¶

circular_linear_correlation(angles: NDArray[float64], linear_values: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad') -> tuple[float, float]

Circular-linear correlation coefficient.

Computes the correlation between a circular and a linear variable using the Mardia & Jupp formula.

Parameters:

Name	Type	Description	Default
`angles`	`ndarray of shape (n_samples,)`	Circular variable (e.g., spike phases) in radians or degrees.	required
`linear_values`	`ndarray of shape (n_samples,)`	Linear variable (e.g., positions, time).	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles.	`'rad'`

Returns:

Name	Type	Description
`r`	`float`	Correlation coefficient in [0, 1]. Note: always non-negative.
`pval`	`float`	P-value from chi-squared distribution with 2 df.

Raises:

Type	Description
`ValueError`	If arrays have different lengths or insufficient samples.

See Also

phase_position_correlation : Alias with neuroscience naming. circular_circular_correlation : For two circular variables.

Notes

The circular-linear correlation is computed using the formula from Mardia & Jupp (2000, Section 11.3):

r^2 = (r_xs^2 + r_xc^2 - 2*r_xs*r_xc*r_cs) / (1 - r_cs^2)

where: - r_xs: Pearson correlation between linear variable x and sin(angles) - r_xc: Pearson correlation between linear variable x and cos(angles) - r_cs: Pearson correlation between cos(angles) and sin(angles)

The result is always non-negative because r is the square root of r^2.

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_linear_correlation
>>> positions = np.linspace(0, 100, 50)
>>> phases = np.linspace(0, 2 * np.pi, 50)  # Linear relationship
>>> r, p = circular_linear_correlation(phases, positions)
>>> r > 0.5 and p < 0.05  # Significant correlation
True

Source code in src/neurospatial/stats/circular.py

def circular_linear_correlation(
    angles: NDArray[np.float64],
    linear_values: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
) -> tuple[float, float]:
    """
    Circular-linear correlation coefficient.

    Computes the correlation between a circular and a linear variable
    using the Mardia & Jupp formula.

    Parameters
    ----------
    angles : ndarray of shape (n_samples,)
        Circular variable (e.g., spike phases) in radians or degrees.
    linear_values : ndarray of shape (n_samples,)
        Linear variable (e.g., positions, time).
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles.

    Returns
    -------
    r : float
        Correlation coefficient in [0, 1]. Note: always non-negative.
    pval : float
        P-value from chi-squared distribution with 2 df.

    Raises
    ------
    ValueError
        If arrays have different lengths or insufficient samples.

    See Also
    --------
    phase_position_correlation : Alias with neuroscience naming.
    circular_circular_correlation : For two circular variables.

    Notes
    -----
    The circular-linear correlation is computed using the formula from
    Mardia & Jupp (2000, Section 11.3):

        r^2 = (r_xs^2 + r_xc^2 - 2*r_xs*r_xc*r_cs) / (1 - r_cs^2)

    where:
    - r_xs: Pearson correlation between linear variable x and sin(angles)
    - r_xc: Pearson correlation between linear variable x and cos(angles)
    - r_cs: Pearson correlation between cos(angles) and sin(angles)

    The result is always non-negative because r is the square root of r^2.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_linear_correlation
    >>> positions = np.linspace(0, 100, 50)
    >>> phases = np.linspace(0, 2 * np.pi, 50)  # Linear relationship
    >>> r, p = circular_linear_correlation(phases, positions)
    >>> r > 0.5 and p < 0.05  # Significant correlation
    True
    """
    # chi2 is imported at module level; only pearsonr is needed locally.
    from scipy.stats import pearsonr

    # Convert to arrays and radians if needed
    angles = np.asarray(angles, dtype=np.float64)
    linear_values = np.asarray(linear_values, dtype=np.float64)
    angles = _to_radians(angles, angle_unit)

    # Validate paired inputs (handles length mismatch, NaN removal)
    angles, linear_values = _validate_paired_input(
        angles, linear_values, "angles", "linear_values", min_samples=3
    )

    n = len(angles)

    # Compute sin and cos of angles
    sin_angles = np.sin(angles)
    cos_angles = np.cos(angles)

    # Check for degenerate case: no variation in linear variable
    linear_std = np.std(linear_values)
    if linear_std < 1e-10:
        warnings.warn(
            "Linear variable has no variation (constant values). "
            "Circular-linear correlation is undefined. Returning r=0.",
            category=RuntimeWarning,
            stacklevel=2,
        )
        return 0.0, 1.0

    # Compute Pearson correlations
    # r_xs: correlation between x (linear) and sin(angles)
    # r_xc: correlation between x (linear) and cos(angles)
    # r_cs: correlation between cos(angles) and sin(angles)
    r_xs, _ = pearsonr(linear_values, sin_angles)
    r_xc, _ = pearsonr(linear_values, cos_angles)
    r_cs, _ = pearsonr(cos_angles, sin_angles)

    # Handle degenerate case: r_cs near 1 (cos and sin perfectly correlated)
    # This shouldn't happen for real circular data but handle it anyway
    if np.abs(r_cs) > 0.9999:
        warnings.warn(
            "Degenerate case: cos(angles) and sin(angles) are nearly perfectly "
            "correlated. This suggests angles have very limited range. "
            "Returning r=0.",
            category=RuntimeWarning,
            stacklevel=2,
        )
        return 0.0, 1.0

    # Mardia & Jupp formula for circular-linear correlation
    # r^2 = (r_xs^2 + r_xc^2 - 2*r_xs*r_xc*r_cs) / (1 - r_cs^2)
    r_squared = (r_xs**2 + r_xc**2 - 2 * r_xs * r_xc * r_cs) / (1 - r_cs**2)

    # Ensure r_squared is non-negative (numerical issues can make it slightly negative)
    r_squared = max(0.0, r_squared)

    # r is always non-negative by definition
    r = float(np.sqrt(r_squared))

    # P-value from chi-squared distribution with 2 degrees of freedom
    # Test statistic: n * r^2 follows chi-squared(2) under null hypothesis
    chi2_stat = n * r_squared
    pval = float(1.0 - chi2.cdf(chi2_stat, df=2))

    # Ensure p-value is in valid range
    pval = float(np.clip(pval, 0.0, 1.0))

    return r, pval

circular_mean ¶

circular_mean(angles: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad', weights: NDArray[float64] | None = None) -> float

Compute the circular mean (mean direction) of angles.

Parameters:

Name	Type	Description	Default
`angles`	`ndarray of shape (n_samples,)`	Sample of angles.	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles.	`'rad'`
`weights`	`ndarray of shape (n_samples,)`	Weights for each angle. If None, uniform weights are used. Must be the same length as `angles` and non-negative.	`None`

Returns:

Type	Description
`float`	Circular mean in radians, range [-π, π].

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_mean
>>> angles = np.array([0, np.pi / 4, np.pi / 2])
>>> mean = circular_mean(angles)
>>> bool(np.abs(mean - np.pi / 4) < 0.1)
True

Source code in src/neurospatial/stats/circular.py

def circular_mean(
    angles: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
    weights: NDArray[np.float64] | None = None,
) -> float:
    """
    Compute the circular mean (mean direction) of angles.

    Parameters
    ----------
    angles : ndarray of shape (n_samples,)
        Sample of angles.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles.
    weights : ndarray of shape (n_samples,), optional
        Weights for each angle. If None, uniform weights are used.
        Must be the same length as ``angles`` and non-negative.

    Returns
    -------
    float
        Circular mean in radians, range [-π, π].

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_mean
    >>> angles = np.array([0, np.pi / 4, np.pi / 2])
    >>> mean = circular_mean(angles)
    >>> bool(np.abs(mean - np.pi / 4) < 0.1)
    True
    """
    angles = np.asarray(angles, dtype=np.float64)
    angles = _to_radians(angles, angle_unit)

    if len(angles) == 0:
        return np.nan

    if weights is not None:
        # Co-filter NaN angles together with their paired weights (same rule
        # as the weighted Rayleigh path), then compute on the remainder.
        # Non-finite WEIGHTS still raise inside the validator -- only NaN
        # ANGLES are legitimately dropped here.
        angles, weights = _validate_circular_input(
            angles, "angles", min_samples=1, check_range=False, weights=weights
        )
        assert weights is not None  # validator returns weights when given them
        cos_angles = np.cos(angles)
        sin_angles = np.sin(angles)
        weight_sum = np.sum(weights)
        if weight_sum == 0:
            return np.nan
        weights_norm = weights / weight_sum
        mean_cos = np.sum(weights_norm * cos_angles)
        mean_sin = np.sum(weights_norm * sin_angles)
    else:
        cos_angles = np.cos(angles)
        sin_angles = np.sin(angles)
        mean_cos = np.mean(cos_angles)
        mean_sin = np.mean(sin_angles)

    return float(np.arctan2(mean_sin, mean_cos))

circular_variance ¶

circular_variance(angles: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad', weights: NDArray[float64] | None = None) -> float

Compute the circular variance of angles.

The circular variance is defined as V = 1 - R, where R is the mean resultant length.

Parameters:

Name	Type	Description	Default
`angles`	`ndarray of shape (n_samples,)`	Sample of angles.	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles.	`'rad'`
`weights`	`ndarray of shape (n_samples,)`	Weights for each angle. If None, uniform weights are used. Must be the same length as `angles` and non-negative.	`None`

Returns:

Type	Description
`float`	Circular variance in [0, 1]. 0 indicates all angles are identical, 1 indicates maximum dispersion (uniform distribution).

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import circular_variance
>>> # Concentrated distribution - low variance
>>> concentrated = np.array([0.1, 0.15, 0.05, 0.12])
>>> var = circular_variance(concentrated)
>>> var < 0.1
True

Source code in src/neurospatial/stats/circular.py

def circular_variance(
    angles: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
    weights: NDArray[np.float64] | None = None,
) -> float:
    """
    Compute the circular variance of angles.

    The circular variance is defined as V = 1 - R, where R is the mean
    resultant length.

    Parameters
    ----------
    angles : ndarray of shape (n_samples,)
        Sample of angles.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles.
    weights : ndarray of shape (n_samples,), optional
        Weights for each angle. If None, uniform weights are used.
        Must be the same length as ``angles`` and non-negative.

    Returns
    -------
    float
        Circular variance in [0, 1]. 0 indicates all angles are identical,
        1 indicates maximum dispersion (uniform distribution).

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import circular_variance
    >>> # Concentrated distribution - low variance
    >>> concentrated = np.array([0.1, 0.15, 0.05, 0.12])
    >>> var = circular_variance(concentrated)
    >>> var < 0.1
    True
    """
    angles = np.asarray(angles, dtype=np.float64)
    angles = _to_radians(angles, angle_unit)

    if weights is not None:
        # Co-filter NaN angles with their paired weights (same rule as the
        # weighted Rayleigh path). Non-finite WEIGHTS still raise.
        angles, weights = _validate_circular_input(
            angles, "angles", min_samples=1, check_range=False, weights=weights
        )

    r = _mean_resultant_length(angles, weights=weights)
    return float(1.0 - r)

is_modulated ¶

is_modulated(beta_sin: float, beta_cos: float, cov_matrix: NDArray[float64], *, alpha: float = 0.05, min_magnitude: float = 0.2) -> bool

Quick check: Do GLM coefficients show significant circular modulation?

This is a convenience function that combines statistical significance (Wald test p-value) with a practical magnitude threshold.

Parameters:

Name	Type	Description	Default
`beta_sin`	`float`	Coefficient for sin(angle) from fitted GLM.	required
`beta_cos`	`float`	Coefficient for cos(angle) from fitted GLM.	required
`cov_matrix`	`(ndarray, shape(2, 2))`	Covariance matrix for [beta_sin, beta_cos] from GLM fit. Required for significance testing.	required
`alpha`	`float`	Significance level for Wald test.	`0.05`
`min_magnitude`	`float`	Minimum modulation magnitude (amplitude) required. Amplitude = sqrt(beta_sin^2 + beta_cos^2). Default justification (0.2): This threshold provides a balance between sensitivity and specificity for detecting biologically meaningful circular modulation in neural data. Based on typical GLM coefficient scales for head direction cells and phase precession: Weak modulation: 0.1-0.3 (may be noise) Moderate modulation: 0.3-0.6 (reliable signal) Strong modulation: >0.6 (clear tuning) The 0.2 threshold filters out very weak effects while maintaining sensitivity to moderate modulation. Adjust based on your data scale and scientific question.	`0.2`

Returns:

Type	Description
`bool`	True if BOTH: - p-value < alpha (statistically significant) - magnitude >= min_magnitude (practically meaningful) False otherwise.

See Also

circular_basis_metrics : For full analysis with all metrics.

Notes

Why both thresholds?

Statistical significance alone can occur with very weak modulation when sample sizes are large. Requiring both significance AND a minimum magnitude ensures the modulation is both reliable and biologically meaningful.

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import is_modulated
>>> # Strong significant modulation
>>> is_modulated(2.0, 2.0, np.array([[0.01, 0], [0, 0.01]]))
True
>>> # Weak (non-significant) modulation
>>> is_modulated(0.1, 0.1, np.array([[1.0, 0], [0, 1.0]]))
False

Source code in src/neurospatial/stats/circular.py

def is_modulated(
    beta_sin: float,
    beta_cos: float,
    cov_matrix: NDArray[np.float64],
    *,
    alpha: float = 0.05,
    min_magnitude: float = 0.2,
) -> bool:
    """
    Quick check: Do GLM coefficients show significant circular modulation?

    This is a convenience function that combines statistical significance
    (Wald test p-value) with a practical magnitude threshold.

    Parameters
    ----------
    beta_sin : float
        Coefficient for sin(angle) from fitted GLM.
    beta_cos : float
        Coefficient for cos(angle) from fitted GLM.
    cov_matrix : ndarray, shape (2, 2)
        Covariance matrix for [beta_sin, beta_cos] from GLM fit.
        Required for significance testing.
    alpha : float, default=0.05
        Significance level for Wald test.
    min_magnitude : float, default=0.2
        Minimum modulation magnitude (amplitude) required.
        Amplitude = sqrt(beta_sin^2 + beta_cos^2).

        **Default justification (0.2)**: This threshold provides a balance
        between sensitivity and specificity for detecting biologically
        meaningful circular modulation in neural data. Based on typical
        GLM coefficient scales for head direction cells and phase precession:

        - Weak modulation: 0.1-0.3 (may be noise)
        - Moderate modulation: 0.3-0.6 (reliable signal)
        - Strong modulation: >0.6 (clear tuning)

        The 0.2 threshold filters out very weak effects while maintaining
        sensitivity to moderate modulation. Adjust based on your data scale
        and scientific question.

    Returns
    -------
    bool
        True if BOTH:
        - p-value < alpha (statistically significant)
        - magnitude >= min_magnitude (practically meaningful)
        False otherwise.

    See Also
    --------
    circular_basis_metrics : For full analysis with all metrics.

    Notes
    -----
    **Why both thresholds?**

    Statistical significance alone can occur with very weak modulation
    when sample sizes are large. Requiring both significance AND a
    minimum magnitude ensures the modulation is both reliable and
    biologically meaningful.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import is_modulated
    >>> # Strong significant modulation
    >>> is_modulated(2.0, 2.0, np.array([[0.01, 0], [0, 0.01]]))
    True
    >>> # Weak (non-significant) modulation
    >>> is_modulated(0.1, 0.1, np.array([[1.0, 0], [0, 1.0]]))
    False
    """
    amplitude, _, pvalue = circular_basis_metrics(beta_sin, beta_cos, cov_matrix)

    # pvalue should never be None since cov_matrix is required, but handle it
    if pvalue is None:
        return False

    return pvalue < alpha and amplitude >= min_magnitude

mean_resultant_length ¶

mean_resultant_length(angles: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad', weights: NDArray[float64] | None = None) -> float

Compute the mean resultant length (Rayleigh vector length).

The mean resultant length R measures concentration of circular data. It equals 1 if all angles are identical and approaches 0 for uniform distributions.

Parameters:

Name	Type	Description	Default
`angles`	`ndarray of shape (n_samples,)`	Sample of angles.	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles.	`'rad'`
`weights`	`ndarray of shape (n_samples,)`	Weights for each angle. If None, uniform weights are used. Must be the same length as `angles` and non-negative.	`None`

Returns:

Type	Description
`float`	Mean resultant length R in [0, 1].

Notes

The mean resultant length is defined as:

R = |mean(exp(i * angles))| = sqrt(mean(cos)^2 + mean(sin)^2)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import mean_resultant_length
>>> # All angles identical - R should be 1
>>> r = mean_resultant_length(np.array([0, 0, 0, 0]))
>>> bool(np.abs(r - 1.0) < 0.01)
True

Source code in src/neurospatial/stats/circular.py

def mean_resultant_length(
    angles: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
    weights: NDArray[np.float64] | None = None,
) -> float:
    """
    Compute the mean resultant length (Rayleigh vector length).

    The mean resultant length R measures concentration of circular data.
    It equals 1 if all angles are identical and approaches 0 for uniform
    distributions.

    Parameters
    ----------
    angles : ndarray of shape (n_samples,)
        Sample of angles.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles.
    weights : ndarray of shape (n_samples,), optional
        Weights for each angle. If None, uniform weights are used.
        Must be the same length as ``angles`` and non-negative.

    Returns
    -------
    float
        Mean resultant length R in [0, 1].

    Notes
    -----
    The mean resultant length is defined as:

        R = |mean(exp(i * angles))| = sqrt(mean(cos)^2 + mean(sin)^2)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import mean_resultant_length
    >>> # All angles identical - R should be 1
    >>> r = mean_resultant_length(np.array([0, 0, 0, 0]))
    >>> bool(np.abs(r - 1.0) < 0.01)
    True
    """
    angles = np.asarray(angles, dtype=np.float64)
    angles = _to_radians(angles, angle_unit)

    if weights is not None:
        # Co-filter NaN angles with their paired weights (same rule as the
        # weighted Rayleigh path). Non-finite WEIGHTS still raise.
        angles, weights = _validate_circular_input(
            angles, "angles", min_samples=1, check_range=False, weights=weights
        )

    return _mean_resultant_length(angles, weights=weights)

phase_position_correlation ¶

phase_position_correlation(phases: NDArray[float64], positions: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad') -> tuple[float, float]

Correlation between spike phases and positions.

This is an alias for circular_linear_correlation with neuroscience naming conventions.

Parameters:

Name	Type	Description	Default
`phases`	`ndarray of shape (n_spikes,)`	Spike phases relative to LFP theta in radians or degrees.	required
`positions`	`ndarray of shape (n_spikes,)`	Position at each spike.	required
`angle_unit`	`('rad', 'deg')`	Unit of input phases.	`'rad'`

Returns:

Name	Type	Description
`r`	`float`	Correlation coefficient in [0, 1].
`pval`	`float`	P-value from chi-squared distribution.

See Also

circular_linear_correlation : Generic version. phase_precession : Full phase precession analysis with slope.

Source code in src/neurospatial/stats/circular.py

def phase_position_correlation(
    phases: NDArray[np.float64],
    positions: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
) -> tuple[float, float]:
    """
    Correlation between spike phases and positions.

    This is an alias for ``circular_linear_correlation`` with neuroscience
    naming conventions.

    Parameters
    ----------
    phases : ndarray of shape (n_spikes,)
        Spike phases relative to LFP theta in radians or degrees.
    positions : ndarray of shape (n_spikes,)
        Position at each spike.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input phases.

    Returns
    -------
    r : float
        Correlation coefficient in [0, 1].
    pval : float
        P-value from chi-squared distribution.

    See Also
    --------
    circular_linear_correlation : Generic version.
    phase_precession : Full phase precession analysis with slope.
    """
    return circular_linear_correlation(phases, positions, angle_unit=angle_unit)

plot_circular_basis_tuning ¶

plot_circular_basis_tuning(beta_sin: float, beta_cos: float, angles: NDArray[float64] | None = None, rates: NDArray[float64] | None = None, ax: Axes | PolarAxes | None = None, *, intercept: float = 0.0, cov_matrix: NDArray[float64] | None = None, projection: Literal['polar', 'linear'] = 'polar', n_points: int = 100, show_data: bool = False, show_fit: bool = True, show_ci: bool = False, ci: float = 0.95, color: str = 'C0', data_color: str = 'gray', data_alpha: float = 0.5, ci_alpha: float = 0.3, line_kwargs: dict[str, Any] | None = None, scatter_kwargs: dict[str, Any] | None = None) -> Axes | PolarAxes

Plot tuning curve from GLM circular basis coefficients.

Visualizes the fitted circular tuning curve from a GLM with sin/cos basis functions, optionally overlaying raw binned data.

Parameters:

Name	Type	Description	Default
`beta_sin`	`float`	Coefficient for sin(angle) from fitted GLM.	required
`beta_cos`	`float`	Coefficient for cos(angle) from fitted GLM.	required
`angles`	`(array, shape(n_bins))`	Center of each angular bin (radians) for raw data overlay. Required if `show_data=True`.	`None`
`rates`	`(array, shape(n_bins))`	Firing rate or response in each bin. Required if `show_data=True`.	`None`
`ax`	`Axes`	Axes to plot on. If None, creates new figure with appropriate projection.	`None`
`intercept`	`float`	Intercept (baseline) coefficient from GLM. For Poisson GLM, this controls the baseline firing rate: exp(intercept).	`0.0`
`cov_matrix`	`(ndarray, shape(2, 2))`	Covariance matrix for [beta_sin, beta_cos] from GLM fit. Required if `show_ci=True` to compute confidence bands.	`None`
`projection`	`('polar', 'linear')`	Plot projection type.	`'polar'`
`n_points`	`int`	Number of points for smooth fitted curve.	`100`
`show_data`	`bool`	If True, overlay raw data points. Requires `angles` and `rates`.	`False`
`show_fit`	`bool`	If True, show smooth curve from GLM fit.	`True`
`show_ci`	`bool`	If True, show confidence band around fitted curve. Requires `cov_matrix`.	`False`
`ci`	`float`	Confidence level for band (e.g., 0.95 for 95% CI).	`0.95`
`color`	`str`	Color for fitted curve.	`'C0'`
`data_color`	`str`	Color for data points.	`'gray'`
`data_alpha`	`float`	Alpha (transparency) for data points.	`0.5`
`ci_alpha`	`float`	Alpha (transparency) for confidence band fill.	`0.3`
`line_kwargs`	`dict`	Additional kwargs for line plot (fitted curve).	`None`
`scatter_kwargs`	`dict`	Additional kwargs for scatter plot (data points).	`None`

Returns:

Type	Description
`Axes`	The axes object with the plot.

Raises:

Type	Description
`ValueError`	If `show_data=True` but `angles` or `rates` not provided. If `show_ci=True` but `cov_matrix` not provided.

See Also

circular_basis_metrics : Compute amplitude/phase from coefficients. circular_basis : Create design matrix for GLM.

Examples:

Head direction GLM tuning curve:

>>> import numpy as np
>>> from neurospatial.stats.circular import plot_circular_basis_tuning
>>> # After fitting GLM to head direction data
>>> # Coefficients: preferred direction ~45° with strong modulation
>>> beta_sin, beta_cos = 0.7, 0.7  # atan2(0.7, 0.7) ≈ 45°
>>> ax = plot_circular_basis_tuning(
...     beta_sin, beta_cos, intercept=2.0
... )

Source code in src/neurospatial/stats/circular.py

def plot_circular_basis_tuning(
    beta_sin: float,
    beta_cos: float,
    angles: NDArray[np.float64] | None = None,
    rates: NDArray[np.float64] | None = None,
    ax: Axes | PolarAxes | None = None,
    *,
    intercept: float = 0.0,
    cov_matrix: NDArray[np.float64] | None = None,
    projection: Literal["polar", "linear"] = "polar",
    n_points: int = 100,
    show_data: bool = False,
    show_fit: bool = True,
    show_ci: bool = False,
    ci: float = 0.95,
    color: str = "C0",
    data_color: str = "gray",
    data_alpha: float = 0.5,
    ci_alpha: float = 0.3,
    line_kwargs: dict[str, Any] | None = None,
    scatter_kwargs: dict[str, Any] | None = None,
) -> Axes | PolarAxes:
    """
    Plot tuning curve from GLM circular basis coefficients.

    Visualizes the fitted circular tuning curve from a GLM with sin/cos
    basis functions, optionally overlaying raw binned data.

    Parameters
    ----------
    beta_sin : float
        Coefficient for sin(angle) from fitted GLM.
    beta_cos : float
        Coefficient for cos(angle) from fitted GLM.
    angles : array, shape (n_bins,), optional
        Center of each angular bin (radians) for raw data overlay.
        Required if ``show_data=True``.
    rates : array, shape (n_bins,), optional
        Firing rate or response in each bin.
        Required if ``show_data=True``.
    ax : matplotlib.axes.Axes, optional
        Axes to plot on. If None, creates new figure with appropriate projection.
    intercept : float, default=0.0
        Intercept (baseline) coefficient from GLM. For Poisson GLM, this
        controls the baseline firing rate: exp(intercept).
    cov_matrix : ndarray, shape (2, 2), optional
        Covariance matrix for [beta_sin, beta_cos] from GLM fit.
        Required if ``show_ci=True`` to compute confidence bands.
    projection : {'polar', 'linear'}, default='polar'
        Plot projection type.
    n_points : int, default=100
        Number of points for smooth fitted curve.
    show_data : bool, default=False
        If True, overlay raw data points. Requires ``angles`` and ``rates``.
    show_fit : bool, default=True
        If True, show smooth curve from GLM fit.
    show_ci : bool, default=False
        If True, show confidence band around fitted curve. Requires ``cov_matrix``.
    ci : float, default=0.95
        Confidence level for band (e.g., 0.95 for 95% CI).
    color : str, default='C0'
        Color for fitted curve.
    data_color : str, default='gray'
        Color for data points.
    data_alpha : float, default=0.5
        Alpha (transparency) for data points.
    ci_alpha : float, default=0.3
        Alpha (transparency) for confidence band fill.
    line_kwargs : dict, optional
        Additional kwargs for line plot (fitted curve).
    scatter_kwargs : dict, optional
        Additional kwargs for scatter plot (data points).

    Returns
    -------
    matplotlib.axes.Axes
        The axes object with the plot.

    Raises
    ------
    ValueError
        If ``show_data=True`` but ``angles`` or ``rates`` not provided.
        If ``show_ci=True`` but ``cov_matrix`` not provided.

    See Also
    --------
    circular_basis_metrics : Compute amplitude/phase from coefficients.
    circular_basis : Create design matrix for GLM.

    Examples
    --------
    **Head direction GLM tuning curve**:

    >>> import numpy as np
    >>> from neurospatial.stats.circular import plot_circular_basis_tuning
    >>> # After fitting GLM to head direction data
    >>> # Coefficients: preferred direction ~45° with strong modulation
    >>> beta_sin, beta_cos = 0.7, 0.7  # atan2(0.7, 0.7) ≈ 45°
    >>> ax = plot_circular_basis_tuning(
    ...     beta_sin, beta_cos, intercept=2.0
    ... )  # doctest: +SKIP
    """
    import matplotlib.pyplot as plt
    from scipy import stats

    # Validate show_data requirements
    if show_data and (angles is None or rates is None):
        raise ValueError(
            "show_data=True requires both angles and rates arguments.\n\n"
            "To show model fit over raw data:\n"
            "  plot_circular_basis_tuning(beta_sin, beta_cos, "
            "angles=bins, rates=rates, show_data=True)\n\n"
            "To show only the fitted curve:\n"
            "  plot_circular_basis_tuning(beta_sin, beta_cos, show_data=False)"
        )

    # Validate show_ci requirements
    if show_ci and cov_matrix is None:
        raise ValueError(
            "show_ci=True requires cov_matrix argument.\n\n"
            "To show confidence bands:\n"
            "  cov = model.cov_params()[1:3, 1:3]  # Extract sin/cos covariance\n"
            "  plot_circular_basis_tuning(beta_sin, beta_cos, cov_matrix=cov, show_ci=True)"
        )

    # Validate cov_matrix shape if provided
    if cov_matrix is not None:
        cov_matrix = np.asarray(cov_matrix, dtype=np.float64)
        if cov_matrix.shape != (2, 2):
            raise ValueError(
                f"cov_matrix must be a 2x2 matrix for [beta_sin, beta_cos].\n"
                f"Got shape: {cov_matrix.shape}\n\n"
                f"Extract from GLM fit:\n"
                f"  cov = model.cov_params()[1:3, 1:3]  # Rows/cols for sin, cos"
            )

    # Create figure if needed
    if ax is None:
        if projection == "polar":
            _, ax = plt.subplots(subplot_kw={"projection": "polar"})
        else:
            _, ax = plt.subplots()

    # Default kwargs
    line_defaults: dict = {"color": color, "linewidth": 2, "zorder": 2}
    scatter_defaults: dict = {
        "color": data_color,
        "alpha": data_alpha,
        "s": 30,
        "zorder": 1,
    }

    # Merge with user kwargs
    line_kw = {**line_defaults, **(line_kwargs or {})}
    scatter_kw = {**scatter_defaults, **(scatter_kwargs or {})}

    # Generate smooth curve for fit
    theta_smooth = np.linspace(0, 2 * np.pi, n_points, endpoint=False)

    # GLM prediction: lambda(theta) = exp(intercept + beta_cos*cos + beta_sin*sin)
    linear_pred = (
        intercept + beta_cos * np.cos(theta_smooth) + beta_sin * np.sin(theta_smooth)
    )
    # Transform to response scale (Poisson GLM uses exp link)
    fitted_response = np.exp(linear_pred)

    # Compute confidence bands if requested
    ci_lower: NDArray[np.float64] | None = None
    ci_upper: NDArray[np.float64] | None = None

    if show_ci and cov_matrix is not None:
        # Design matrix for each angle: [sin(θ), cos(θ)]
        design = np.column_stack([np.sin(theta_smooth), np.cos(theta_smooth)])

        # Variance of linear predictor at each angle
        var_linear = np.sum((design @ cov_matrix) * design, axis=1)
        se_linear = np.sqrt(np.maximum(var_linear, 0))

        # Z-score for confidence level
        z = stats.norm.ppf((1 + ci) / 2)

        # CI on linear scale, then transform
        linear_lower = linear_pred - z * se_linear
        linear_upper = linear_pred + z * se_linear

        # Transform to response scale
        ci_lower = np.exp(linear_lower)
        ci_upper = np.exp(linear_upper)

    # Close the curve (append first point)
    theta_closed = np.concatenate([theta_smooth, [theta_smooth[0] + 2 * np.pi]])
    fitted_closed = np.concatenate([fitted_response, [fitted_response[0]]])

    # Close CI bands if computed
    if ci_lower is not None and ci_upper is not None:
        ci_lower_closed = np.concatenate([ci_lower, [ci_lower[0]]])
        ci_upper_closed = np.concatenate([ci_upper, [ci_upper[0]]])
    else:
        ci_lower_closed = None
        ci_upper_closed = None

    if projection == "polar":
        polar_ax = cast("PolarAxes", ax)

        # Configure polar plot: 0° at top (North), clockwise direction
        polar_ax.set_theta_zero_location("N")
        polar_ax.set_theta_direction(-1)

        # Plot confidence band (behind curve)
        if show_ci and ci_lower_closed is not None and ci_upper_closed is not None:
            polar_ax.fill_between(
                theta_closed,
                ci_lower_closed,
                ci_upper_closed,
                color=color,
                alpha=ci_alpha,
                zorder=0,
            )

        # Plot fitted curve
        if show_fit:
            polar_ax.plot(theta_closed, fitted_closed, **line_kw)

        # Plot data points
        if show_data and angles is not None and rates is not None:
            angles_arr = np.asarray(angles, dtype=np.float64).ravel()
            rates_arr = np.asarray(rates, dtype=np.float64).ravel()
            polar_ax.scatter(angles_arr, rates_arr, **scatter_kw)

    else:
        # Linear projection
        # Plot confidence band (behind curve)
        if show_ci and ci_lower_closed is not None and ci_upper_closed is not None:
            ax.fill_between(
                theta_closed,
                ci_lower_closed,
                ci_upper_closed,
                color=color,
                alpha=ci_alpha,
                zorder=0,
            )

        # Plot fitted curve
        if show_fit:
            ax.plot(theta_closed, fitted_closed, **line_kw)

        # Plot data points
        if show_data and angles is not None and rates is not None:
            angles_arr = np.asarray(angles, dtype=np.float64).ravel()
            rates_arr = np.asarray(rates, dtype=np.float64).ravel()
            ax.scatter(angles_arr, rates_arr, **scatter_kw)

        ax.set_xlabel("Angle (rad)")
        ax.set_ylabel("Response")
        ax.set_xlim(0, 2 * np.pi)

    return ax

rayleigh_test ¶

rayleigh_test(angles: NDArray[float64], *, angle_unit: Literal['rad', 'deg'] = 'rad', weights: NDArray[float64] | None = None) -> tuple[float, float]

Rayleigh test for non-uniformity of circular data.

Tests H0: angles are uniformly distributed on the circle. Rejection indicates a preferred direction exists.

Parameters:

Name	Type	Description	Default
`angles`	`ndarray of shape (n_samples,)`	Sample of angles. Values outside [0, 2pi] (or [0, 360] for degrees) will be wrapped.	required
`angle_unit`	`('rad', 'deg')`	Unit of input angles.	`'rad'`
`weights`	`ndarray of shape (n_samples,)`	Weights for each angle, interpreted as counts/frequencies (e.g., spike counts per phase bin). If None, uniform unit weights are used. Under this count convention, weighting an angle by an integer `k` is equivalent to observing that angle `k` times; the test statistic is `z = sum(weights) * R**2` (see Notes).	`None`

Returns:

Name	Type	Description
`z`	`float`	Rayleigh z-statistic (n * R^2 where R is mean resultant length). Range: [0, n]. Higher values indicate stronger directionality.
`pval`	`float`	P-value from Rayleigh approximation with finite-sample correction. Small p-values (< 0.05) indicate non-uniform distribution.

Raises:

Type	Description
`ValueError`	If angles array is empty, all NaN, or too short; if `weights` is not the same length as `angles` (a length-1 array is rejected, not broadcast) or contains negative values; or if the total weight (sum of weights) is zero.

See Also

circular_linear_correlation : Test correlation with linear variable.

Notes

The unweighted statistic is z = n * R**2 where R is the mean resultant length and n the number of angles. A finite-sample correction from Mardia & Jupp (2000, Section 5.3.2, p. 94) is applied for accurate p-values when the (effective) sample size is < 50.

Weighted variant. When weights are supplied this implements the count-based (grouped-data) weighted Rayleigh test of Mardia & Jupp (2000, Section 5.3.5): weights are treated as frequencies, the resultant length is the weighted mean resultant length R = |sum_i w_i exp(i theta_i)| / sum_i w_i, and the statistic is z = (sum_i w_i) * R**2 with the finite-sample correction evaluated at n = sum_i w_i. With unit weights this reduces exactly to the unweighted statistic, and integer weights match physically replicating each angle by its weight. This is intentionally not the effective-sample-size form n_eff = (sum w)**2 / sum(w**2).

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import rayleigh_test
>>> # Uniform distribution - expect high p-value
>>> uniform_angles = np.linspace(0, 2 * np.pi, 100, endpoint=False)
>>> z, p = rayleigh_test(uniform_angles)
>>> p > 0.5
True

Source code in src/neurospatial/stats/circular.py

def rayleigh_test(
    angles: NDArray[np.float64],
    *,
    angle_unit: Literal["rad", "deg"] = "rad",
    weights: NDArray[np.float64] | None = None,
) -> tuple[float, float]:
    """
    Rayleigh test for non-uniformity of circular data.

    Tests H0: angles are uniformly distributed on the circle.
    Rejection indicates a preferred direction exists.

    Parameters
    ----------
    angles : ndarray of shape (n_samples,)
        Sample of angles. Values outside [0, 2pi] (or [0, 360] for degrees)
        will be wrapped.
    angle_unit : {'rad', 'deg'}, default='rad'
        Unit of input angles.
    weights : ndarray of shape (n_samples,), optional
        Weights for each angle, interpreted as **counts/frequencies**
        (e.g., spike counts per phase bin). If None, uniform unit weights
        are used. Under this count convention, weighting an angle by an
        integer ``k`` is equivalent to observing that angle ``k`` times;
        the test statistic is ``z = sum(weights) * R**2`` (see Notes).

    Returns
    -------
    z : float
        Rayleigh z-statistic (n * R^2 where R is mean resultant length).
        Range: [0, n]. Higher values indicate stronger directionality.
    pval : float
        P-value from Rayleigh approximation with finite-sample correction.
        Small p-values (< 0.05) indicate non-uniform distribution.

    Raises
    ------
    ValueError
        If angles array is empty, all NaN, or too short; if ``weights`` is
        not the same length as ``angles`` (a length-1 array is rejected,
        not broadcast) or contains negative values; or if the total weight
        (sum of weights) is zero.

    See Also
    --------
    circular_linear_correlation : Test correlation with linear variable.

    Notes
    -----
    The unweighted statistic is ``z = n * R**2`` where ``R`` is the mean
    resultant length and ``n`` the number of angles. A finite-sample
    correction from Mardia & Jupp (2000, Section 5.3.2, p. 94) is applied
    for accurate p-values when the (effective) sample size is < 50.

    **Weighted variant.** When ``weights`` are supplied this implements the
    count-based (grouped-data) weighted Rayleigh test of Mardia & Jupp
    (2000, Section 5.3.5): weights are treated as frequencies, the resultant
    length is the weighted mean resultant length
    ``R = |sum_i w_i exp(i theta_i)| / sum_i w_i``, and the statistic is
    ``z = (sum_i w_i) * R**2`` with the finite-sample correction evaluated
    at ``n = sum_i w_i``. With unit weights this reduces exactly to the
    unweighted statistic, and integer weights match physically replicating
    each angle by its weight. This is intentionally *not* the
    effective-sample-size form ``n_eff = (sum w)**2 / sum(w**2)``.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import rayleigh_test
    >>> # Uniform distribution - expect high p-value
    >>> uniform_angles = np.linspace(0, 2 * np.pi, 100, endpoint=False)
    >>> z, p = rayleigh_test(uniform_angles)
    >>> p > 0.5
    True
    """
    # Convert to radians if needed
    angles = np.asarray(angles, dtype=np.float64)
    angles = _to_radians(angles, angle_unit)

    # Validate input (handles NaN drop + Inf + min samples) and co-filter
    # weights so angle/weight alignment survives NaN removal. Weight
    # non-negativity and equal-length are enforced inside the validator.
    #
    # The validator's ``min_samples`` gate counts the number of *angles*
    # (distinct samples). That is correct for UNWEIGHTED input but wrong for
    # WEIGHTED (count) input: a genuinely well-tuned cell can put 100 spikes
    # into just 1-2 angular bins, so the effective sample size is the total
    # weight (sum of counts), not the number of occupied bins. For the
    # weighted path we therefore only require >= 1 angle here and enforce the
    # real minimum-sample-size check on ``sum(weights)`` below (via n_eff).
    validator_min_samples = 1 if weights is not None else 3
    angles, weights = _validate_circular_input(
        angles,
        "angles",
        min_samples=validator_min_samples,
        check_range=False,
        weights=weights,
    )

    n = len(angles)

    # Weighted mean resultant length (weights already validated above).
    r_mean = _mean_resultant_length(angles, weights=weights)

    # Effective count for the z-statistic and finite-sample correction.
    # Count-based (grouped-data) weighted Rayleigh (Mardia & Jupp 2000,
    # Section 5.3.5): n_eff = sum(weights); for unit weights, n_eff = n.
    if weights is not None:
        n_eff: float = float(np.sum(weights))
    else:
        n_eff = float(n)

    # Guard: all-zero (or empty) weights give n_eff == 0. The weighted R
    # is already NaN in that case; surface an actionable error instead of a
    # bare ZeroDivisionError in the finite-sample correction below.
    if n_eff <= 0:
        raise ValueError(
            "Total weight (sum of weights) is zero; cannot run the weighted "
            "Rayleigh test. At least one angle must carry positive weight.\n"
            "Fix: check that your per-angle counts/occupancy are not all zero."
        )

    # Minimum effective sample size. For weighted (count) data the effective
    # sample size is the total weight, so we enforce the same minimum (3) on
    # ``n_eff`` that the validator enforces on the angle count for unweighted
    # data. This rejects genuinely-insufficient data (e.g. 2 total spikes)
    # while accepting a strongly-tuned cell whose many spikes land in only
    # 1-2 angular bins.
    if weights is not None and n_eff < 3:
        raise ValueError(
            f"Need at least 3 total weight (sum of counts) for the weighted "
            f"Rayleigh test. Got sum(weights) = {n_eff:g}.\n"
            f"Fix: provide more spikes/counts or use a different analysis method."
        )

    # Rayleigh z-statistic: z = n_eff * R^2.
    z = n_eff * r_mean**2

    # P-value with finite-sample correction (Mardia & Jupp, p. 94)
    # For large n, p = exp(-z) is a good approximation
    # For small n, we use the correction formula
    pval = float(np.exp(-z))

    # Apply finite-sample correction for more accurate p-values
    # From Mardia & Jupp (2000), Section 5.3.2, equation 5.3.6
    if n_eff < 50:
        # Correction terms
        term1 = (2 * z - z**2) / (4 * n_eff)
        term2 = (24 * z - 132 * z**2 + 76 * z**3 - 9 * z**4) / (288 * n_eff**2)
        pval = pval * (1 + term1 - term2)

    # Ensure p-value is in valid range
    pval = float(np.clip(pval, 0.0, 1.0))

    return float(z), pval

wrap_angle ¶

wrap_angle(angle: NDArray[float64]) -> NDArray[np.float64]

Wrap angle to (-π, π].

Parameters:

Name	Type	Description	Default
`angle`	`NDArray[float64]`	Angles in radians.	required

Returns:

Type	Description
`NDArray[float64]`	Angles wrapped to (-π, π].

Examples:

>>> import numpy as np
>>> from neurospatial.stats.circular import wrap_angle
>>> wrap_angle(np.array([3 * np.pi / 2]))
array([-1.57079633])

Source code in src/neurospatial/stats/circular.py

def wrap_angle(angle: NDArray[np.float64]) -> NDArray[np.float64]:
    """
    Wrap angle to (-π, π].

    Parameters
    ----------
    angle : NDArray[np.float64]
        Angles in radians.

    Returns
    -------
    NDArray[np.float64]
        Angles wrapped to (-π, π].

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.circular import wrap_angle
    >>> wrap_angle(np.array([3 * np.pi / 2]))
    array([-1.57079633])
    """
    return (angle + np.pi) % (2 * np.pi) - np.pi

compute_shuffle_pvalue ¶

compute_shuffle_pvalue(observed: float, null_scores: NDArray[float64], *, tail: Literal['greater', 'less', 'two-sided'] = 'greater') -> float

Compute Monte Carlo p-value with Phipson-Smyth correction.

Computes the probability of observing a score at least as extreme as the observed value under the null distribution using the formula (k + 1) / (n + 1), which provides an unbiased estimate and avoids p-values of exactly zero.

Parameters:

Name	Type	Description	Default
`observed`	`float`	The observed score from the original (non-shuffled) data.	required
`null_scores`	`NDArray[float64]`	Array of scores from shuffled data, forming the null distribution.	required
`tail`	`('greater', 'less', 'two-sided')`	Direction of the test: "greater": Tests if observed is significantly greater than null (counts null >= observed). "less": Tests if observed is significantly less than null (counts null <= observed). "two-sided": Tests if observed is significantly different from null in either direction. Computed as 2 * min(p_greater, p_less), capped at 1.0.	`"greater"`

Returns:

Type	Description
`float`	The Monte Carlo p-value in the range (0, 1].

Notes

The Phipson-Smyth correction (k + 1) / (n + 1) ensures that:

The p-value is never exactly zero, which is important for downstream corrections (e.g., Bonferroni, FDR).
The estimator is unbiased for uniformly distributed test statistics.

For replay analysis, "greater" is typically used because we expect significant sequences to have higher scores than shuffled controls.

Non-finite null scores (NaN or Inf) are dropped with a UserWarning and excluded from n before the p-value is computed; otherwise they would silently bias the p-value toward significance. If no finite null scores remain, a ValueError is raised.

References

.. [1] Phipson, B., & Smyth, G. K. (2010). Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, 9(1).

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import compute_shuffle_pvalue

>>> # Observed score higher than all null values
>>> observed = 10.0
>>> null = np.array([1.0, 2.0, 3.0, 4.0])
>>> compute_shuffle_pvalue(observed, null, tail="greater")
0.2

>>> # Two-sided test
>>> compute_shuffle_pvalue(observed, null, tail="two-sided")
0.4

See Also

compute_shuffle_zscore : Compute z-score from null distribution. ShuffleTestResult : Container for shuffle test results.

Source code in src/neurospatial/stats/shuffle.py

def compute_shuffle_pvalue(
    observed: float,
    null_scores: NDArray[np.float64],
    *,
    tail: Literal["greater", "less", "two-sided"] = "greater",
) -> float:
    """Compute Monte Carlo p-value with Phipson-Smyth correction.

    Computes the probability of observing a score at least as extreme as
    the observed value under the null distribution using the formula
    (k + 1) / (n + 1), which provides an unbiased estimate and avoids
    p-values of exactly zero.

    Parameters
    ----------
    observed : float
        The observed score from the original (non-shuffled) data.
    null_scores : NDArray[np.float64]
        Array of scores from shuffled data, forming the null distribution.
    tail : {"greater", "less", "two-sided"}, default="greater"
        Direction of the test:

        - "greater": Tests if observed is significantly greater than null
          (counts null >= observed).
        - "less": Tests if observed is significantly less than null
          (counts null <= observed).
        - "two-sided": Tests if observed is significantly different from
          null in either direction. Computed as 2 * min(p_greater, p_less),
          capped at 1.0.

    Returns
    -------
    float
        The Monte Carlo p-value in the range (0, 1].

    Notes
    -----
    The Phipson-Smyth correction (k + 1) / (n + 1) ensures that:

    1. The p-value is never exactly zero, which is important for downstream
       corrections (e.g., Bonferroni, FDR).
    2. The estimator is unbiased for uniformly distributed test statistics.

    For replay analysis, "greater" is typically used because we expect
    significant sequences to have higher scores than shuffled controls.

    Non-finite null scores (NaN or Inf) are dropped with a ``UserWarning``
    and excluded from ``n`` before the p-value is computed; otherwise they
    would silently bias the p-value toward significance. If no finite null
    scores remain, a ``ValueError`` is raised.

    References
    ----------
    .. [1] Phipson, B., & Smyth, G. K. (2010). Permutation P-values should
           never be zero: calculating exact P-values when permutations are
           randomly drawn. Statistical Applications in Genetics and Molecular
           Biology, 9(1).

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import compute_shuffle_pvalue

    >>> # Observed score higher than all null values
    >>> observed = 10.0
    >>> null = np.array([1.0, 2.0, 3.0, 4.0])
    >>> compute_shuffle_pvalue(observed, null, tail="greater")
    0.2

    >>> # Two-sided test
    >>> compute_shuffle_pvalue(observed, null, tail="two-sided")
    0.4

    See Also
    --------
    compute_shuffle_zscore : Compute z-score from null distribution.
    ShuffleTestResult : Container for shuffle test results.
    """
    if not np.isfinite(observed):
        raise ValueError(
            f"observed must be a finite value, got {observed!r}. A non-finite "
            "observed score makes all comparisons against the null False, which "
            "would silently return the floor p-value 1/(n+1)."
        )

    null_scores = np.asarray(null_scores, dtype=np.float64).ravel()
    finite = np.isfinite(null_scores)
    n_total = null_scores.size
    if not finite.all():
        n_dropped = int(n_total - finite.sum())
        warnings.warn(
            f"Dropped {n_dropped} non-finite null score(s) out of {n_total} "
            f"before computing the p-value. Using {int(finite.sum())} finite "
            f"null scores; non-finite nulls would otherwise bias the p-value "
            f"toward significance.",
            category=UserWarning,
            stacklevel=2,
        )
        null_scores = null_scores[finite]

    n = len(null_scores)
    if n == 0:
        raise ValueError(
            "null_scores contains no finite values; cannot compute a "
            "shuffle p-value. Check that the shuffle produced valid scores."
        )

    if tail == "greater":
        # Count how many null values are >= observed
        k = int(np.sum(null_scores >= observed))
        return float((k + 1) / (n + 1))

    elif tail == "less":
        # Count how many null values are <= observed
        k = int(np.sum(null_scores <= observed))
        return float((k + 1) / (n + 1))

    elif tail == "two-sided":
        # Compute both tails and take 2 * min, capped at 1.0
        k_greater = int(np.sum(null_scores >= observed))
        k_less = int(np.sum(null_scores <= observed))
        p_greater = (k_greater + 1) / (n + 1)
        p_less = (k_less + 1) / (n + 1)
        return float(min(2 * min(p_greater, p_less), 1.0))

    else:
        raise ValueError(
            f"tail must be 'greater', 'less', or 'two-sided', got '{tail}'"
        )

compute_shuffle_zscore ¶

compute_shuffle_zscore(observed: float, null_scores: NDArray[float64]) -> float

Compute z-score of observed value relative to null distribution.

Calculates the standard score: (observed - mean(null)) / std(null). Returns NaN if the null distribution has zero variance.

Parameters:

Name	Type	Description	Default
`observed`	`float`	The observed score from the original (non-shuffled) data.	required
`null_scores`	`NDArray[float64]`	Array of scores from shuffled data, forming the null distribution.	required

Returns:

Type	Description
`float`	The z-score (standard score). NaN if null has zero variance or contains only one value.

Notes

The z-score indicates how many standard deviations the observed value is from the mean of the null distribution. Positive values indicate the observed is above the null mean; negative values indicate below.

Common interpretation thresholds:

|z| > 1.96: Approximately p < 0.05 (two-tailed) if null is normal
|z| > 2.58: Approximately p < 0.01 (two-tailed) if null is normal
|z| > 3.0: Strong evidence against null hypothesis

Unlike p-values, z-scores are useful for comparing effect sizes across different analyses.

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import compute_shuffle_zscore

>>> null = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
>>> observed = 5.0
>>> z = compute_shuffle_zscore(observed, null)
>>> bool(np.isclose(z, (5.0 - 3.0) / np.std(null)))
True

>>> # Zero variance returns NaN
>>> compute_shuffle_zscore(5.0, np.array([3.0, 3.0, 3.0]))
nan

See Also

compute_shuffle_pvalue : Compute Monte Carlo p-value from null distribution. ShuffleTestResult : Container for shuffle test results.

Source code in src/neurospatial/stats/shuffle.py

def compute_shuffle_zscore(
    observed: float,
    null_scores: NDArray[np.float64],
) -> float:
    """Compute z-score of observed value relative to null distribution.

    Calculates the standard score: (observed - mean(null)) / std(null).
    Returns NaN if the null distribution has zero variance.

    Parameters
    ----------
    observed : float
        The observed score from the original (non-shuffled) data.
    null_scores : NDArray[np.float64]
        Array of scores from shuffled data, forming the null distribution.

    Returns
    -------
    float
        The z-score (standard score). NaN if null has zero variance or
        contains only one value.

    Notes
    -----
    The z-score indicates how many standard deviations the observed value
    is from the mean of the null distribution. Positive values indicate
    the observed is above the null mean; negative values indicate below.

    Common interpretation thresholds:

    - |z| > 1.96: Approximately p < 0.05 (two-tailed) if null is normal
    - |z| > 2.58: Approximately p < 0.01 (two-tailed) if null is normal
    - |z| > 3.0: Strong evidence against null hypothesis

    Unlike p-values, z-scores are useful for comparing effect sizes across
    different analyses.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import compute_shuffle_zscore

    >>> null = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
    >>> observed = 5.0
    >>> z = compute_shuffle_zscore(observed, null)
    >>> bool(np.isclose(z, (5.0 - 3.0) / np.std(null)))
    True

    >>> # Zero variance returns NaN
    >>> compute_shuffle_zscore(5.0, np.array([3.0, 3.0, 3.0]))
    nan

    See Also
    --------
    compute_shuffle_pvalue : Compute Monte Carlo p-value from null distribution.
    ShuffleTestResult : Container for shuffle test results.
    """
    mean_null = float(np.mean(null_scores))
    std_null = float(np.std(null_scores))

    # Handle zero variance (all values equal) or single value
    if std_null == 0.0:
        return float("nan")

    return float((observed - mean_null) / std_null)

shuffle_cell_identity ¶

shuffle_cell_identity(spike_counts: NDArray[int64], encoding_models: NDArray[float64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[tuple[NDArray[np.int64], NDArray[np.float64]], None, None]

Shuffle mapping between spike trains and place fields.

Primary test for spatial code coherence. Disrupts the learned relationship between a neuron's activity and its encoded spatial location by randomly permuting which spike train is associated with which place field.

This shuffle randomly permutes columns (neuron axis) of spike_counts, effectively reassigning which place field goes with which spike train. The encoding models are returned unchanged.

Parameters:

Name	Type	Description	Default
`spike_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts per neuron per time bin.	required
`encoding_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Firing rate maps (place fields) for each neuron.	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts with neuron identities permuted (columns shuffled).
`encoding_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Original encoding models (unchanged, same object).

Raises:

Type	Description
`ValueError`	If the number of neurons in `spike_counts` (columns) does not match the number of rows in `encoding_models`. Row i of `encoding_models` is the place field for neuron i (column i of `spike_counts`); a mismatch would silently yield wrong decodes.

Notes

Randomly permutes columns of spike_counts (neuron axis)
Encoding models remain fixed (same object returned each iteration)
Equivalent to randomly reassigning which place field goes with which spike train
Preserves spike counts per neuron per time bin
Preserves total spikes per time bin
Caution: Can introduce noise if firing rates differ greatly between cells

Alternative implementation (equivalent): Instead of shuffling spike_counts columns, shuffle encoding_models rows. This yields the same decoded posteriors.

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_cell_identity

>>> spike_counts = np.array([[0, 1, 2], [2, 0, 1]], dtype=np.int64)
>>> encoding_models = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
>>> for i, (shuffled, models) in enumerate(
...     shuffle_cell_identity(spike_counts, encoding_models, n_shuffles=3, rng=42)
... ):
...     print(
...         f"Shuffle {i}: sum={shuffled.sum()}, models unchanged={models is encoding_models}"
...     )
Shuffle 0: sum=6, models unchanged=True
Shuffle 1: sum=6, models unchanged=True
Shuffle 2: sum=6, models unchanged=True

See Also

shuffle_place_fields_circular : Circular shift of place fields shuffle_time_bins : Temporal order shuffle

Source code in src/neurospatial/stats/shuffle.py

def shuffle_cell_identity(
    spike_counts: NDArray[np.int64],
    encoding_models: NDArray[np.float64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[tuple[NDArray[np.int64], NDArray[np.float64]], None, None]:
    """Shuffle mapping between spike trains and place fields.

    **Primary test for spatial code coherence.** Disrupts the learned
    relationship between a neuron's activity and its encoded spatial location
    by randomly permuting which spike train is associated with which place
    field.

    This shuffle randomly permutes columns (neuron axis) of spike_counts,
    effectively reassigning which place field goes with which spike train.
    The encoding models are returned unchanged.

    Parameters
    ----------
    spike_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts per neuron per time bin.
    encoding_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Firing rate maps (place fields) for each neuron.
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts with neuron identities permuted (columns shuffled).
    encoding_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Original encoding models (unchanged, same object).

    Raises
    ------
    ValueError
        If the number of neurons in ``spike_counts`` (columns) does not
        match the number of rows in ``encoding_models``. Row i of
        ``encoding_models`` is the place field for neuron i (column i of
        ``spike_counts``); a mismatch would silently yield wrong decodes.

    Notes
    -----
    - Randomly permutes columns of spike_counts (neuron axis)
    - Encoding models remain fixed (same object returned each iteration)
    - Equivalent to randomly reassigning which place field goes with which
      spike train
    - Preserves spike counts per neuron per time bin
    - Preserves total spikes per time bin
    - Caution: Can introduce noise if firing rates differ greatly between
      cells

    Alternative implementation (equivalent):
        Instead of shuffling spike_counts columns, shuffle encoding_models
        rows. This yields the same decoded posteriors.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_cell_identity

    >>> spike_counts = np.array([[0, 1, 2], [2, 0, 1]], dtype=np.int64)
    >>> encoding_models = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]])
    >>> for i, (shuffled, models) in enumerate(
    ...     shuffle_cell_identity(spike_counts, encoding_models, n_shuffles=3, rng=42)
    ... ):
    ...     print(
    ...         f"Shuffle {i}: sum={shuffled.sum()}, models unchanged={models is encoding_models}"
    ...     )
    Shuffle 0: sum=6, models unchanged=True
    Shuffle 1: sum=6, models unchanged=True
    Shuffle 2: sum=6, models unchanged=True

    See Also
    --------
    shuffle_place_fields_circular : Circular shift of place fields
    shuffle_time_bins : Temporal order shuffle
    """
    generator = _ensure_rng(rng)
    n_neurons = spike_counts.shape[1]

    if encoding_models.shape[0] != n_neurons:
        raise ValueError(
            f"spike_counts has {n_neurons} neuron(s) (columns) but "
            f"encoding_models has {encoding_models.shape[0]} (rows). These must "
            f"match: row i of encoding_models is the place field for neuron i "
            f"(column i of spike_counts).\n"
            f"Fix: pass spike_counts shaped (n_time_bins, n_neurons) and "
            f"encoding_models shaped (n_neurons, n_bins)."
        )

    for _ in range(n_shuffles):
        # Generate a random permutation of column indices
        perm = generator.permutation(n_neurons)
        # Apply permutation to columns (neuron axis)
        yield spike_counts[:, perm].copy(), encoding_models

shuffle_place_fields_circular ¶

shuffle_place_fields_circular(encoding_models: NDArray[float64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.float64], None, None]

Circularly shift each place field by a random amount.

Conservative null hypothesis. Preserves individual cell spiking properties and local place field structure while disrupting spatial relationships between neurons.

Each neuron's place field is independently shifted by a random amount along the position axis. This preserves the shape of each place field but destroys the spatial relationships between neurons.

Parameters:

Name	Type	Description	Default
`encoding_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Firing rate maps (place fields) for each neuron.	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Place fields with each row circularly shifted by a random amount.

Notes

Each neuron's place field is shifted independently
Preserves the shape of each place field
Destroys spatial relationships between neurons
More conservative than cell identity shuffle
For 2D environments: consider shuffle_place_fields_circular_2d

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_place_fields_circular

>>> encoding_models = np.array(
...     [
...         [1.0, 2.0, 3.0, 4.0],
...         [4.0, 3.0, 2.0, 1.0],
...     ]
... )
>>> for i, shuffled in enumerate(
...     shuffle_place_fields_circular(encoding_models, n_shuffles=3, rng=42)
... ):
...     print(f"Shuffle {i}: shape={shuffled.shape}")
Shuffle 0: shape=(2, 4)
Shuffle 1: shape=(2, 4)
Shuffle 2: shape=(2, 4)

See Also

shuffle_place_fields_circular_2d : 2D circular shift for 2D environments shuffle_cell_identity : Shuffle neuron-to-place-field mapping

Source code in src/neurospatial/stats/shuffle.py

def shuffle_place_fields_circular(
    encoding_models: NDArray[np.float64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Circularly shift each place field by a random amount.

    **Conservative null hypothesis.** Preserves individual cell spiking
    properties and local place field structure while disrupting spatial
    relationships between neurons.

    Each neuron's place field is independently shifted by a random amount
    along the position axis. This preserves the shape of each place field
    but destroys the spatial relationships between neurons.

    Parameters
    ----------
    encoding_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Firing rate maps (place fields) for each neuron.
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Place fields with each row circularly shifted by a random amount.

    Notes
    -----
    - Each neuron's place field is shifted independently
    - Preserves the shape of each place field
    - Destroys spatial relationships between neurons
    - More conservative than cell identity shuffle
    - For 2D environments: consider ``shuffle_place_fields_circular_2d``

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_place_fields_circular

    >>> encoding_models = np.array(
    ...     [
    ...         [1.0, 2.0, 3.0, 4.0],
    ...         [4.0, 3.0, 2.0, 1.0],
    ...     ]
    ... )
    >>> for i, shuffled in enumerate(
    ...     shuffle_place_fields_circular(encoding_models, n_shuffles=3, rng=42)
    ... ):
    ...     print(f"Shuffle {i}: shape={shuffled.shape}")
    Shuffle 0: shape=(2, 4)
    Shuffle 1: shape=(2, 4)
    Shuffle 2: shape=(2, 4)

    See Also
    --------
    shuffle_place_fields_circular_2d : 2D circular shift for 2D environments
    shuffle_cell_identity : Shuffle neuron-to-place-field mapping
    """
    generator = _ensure_rng(rng)
    n_neurons, n_bins = encoding_models.shape

    # Pre-compute column indices for vectorized roll
    col_indices = np.arange(n_bins)

    for _ in range(n_shuffles):
        # Generate random shift amounts for each neuron
        shifts = np.asarray(generator.integers(0, n_bins, size=n_neurons))

        # Vectorized circular shift using advanced indexing
        # For each row i, we want columns (col_indices - shifts[i]) % n_bins
        # Shape: (n_neurons, n_bins)
        shifted_indices = (col_indices - shifts[:, np.newaxis]) % n_bins
        shuffled = encoding_models[np.arange(n_neurons)[:, np.newaxis], shifted_indices]

        yield shuffled

shuffle_place_fields_circular_2d ¶

shuffle_place_fields_circular_2d(encoding_models: NDArray[float64], env: Environment, *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.float64], None, None]

Circularly shift 2D place fields in both dimensions.

For 2D environments, shifts place fields in both x and y dimensions. This is the 2D analog of shuffle_place_fields_circular.

Each neuron's place field is reshaped to a 2D grid, shifted by random amounts in both dimensions, and flattened back to 1D.

Parameters:

Name	Type	Description	Default
`encoding_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Firing rate maps (place fields) for each neuron.	required
`env`	`Environment`	2D environment with grid layout (provides `grid_shape`).	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_models`	`(NDArray[float64], shape(n_neurons, n_bins))`	Place fields with 2D circular shifts applied.

Raises:

Type	Description
`ValueError`	If environment is not 2D or doesn't have grid layout.

Examples:

>>> import numpy as np
>>> from neurospatial import Environment
>>> from neurospatial.stats.shuffle import shuffle_place_fields_circular_2d

>>> positions = np.random.default_rng(42).uniform(0, 10, (100, 2))
>>> env = Environment.from_samples(positions, bin_size=2.0)
>>> encoding_models = np.random.default_rng(42).random(
...     (3, env.n_bins)
... )
>>> for i, shuffled in enumerate(
...     shuffle_place_fields_circular_2d(encoding_models, env, n_shuffles=3, rng=42)
... ):
...     print(f"Shuffle {i}: shape={shuffled.shape}")
Shuffle 0: shape=(3, 36)
Shuffle 1: shape=(3, 36)
Shuffle 2: shape=(3, 36)

See Also

shuffle_place_fields_circular : 1D circular shift shuffle_cell_identity : Shuffle neuron-to-place-field mapping

Source code in src/neurospatial/stats/shuffle.py

def shuffle_place_fields_circular_2d(
    encoding_models: NDArray[np.float64],
    env: Environment,
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Circularly shift 2D place fields in both dimensions.

    For 2D environments, shifts place fields in both x and y dimensions.
    This is the 2D analog of ``shuffle_place_fields_circular``.

    Each neuron's place field is reshaped to a 2D grid, shifted by random
    amounts in both dimensions, and flattened back to 1D.

    Parameters
    ----------
    encoding_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Firing rate maps (place fields) for each neuron.
    env : Environment
        2D environment with grid layout (provides ``grid_shape``).
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_models : NDArray[np.float64], shape (n_neurons, n_bins)
        Place fields with 2D circular shifts applied.

    Raises
    ------
    ValueError
        If environment is not 2D or doesn't have grid layout.

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial import Environment
    >>> from neurospatial.stats.shuffle import shuffle_place_fields_circular_2d

    >>> positions = np.random.default_rng(42).uniform(0, 10, (100, 2))  # doctest: +SKIP
    >>> env = Environment.from_samples(positions, bin_size=2.0)  # doctest: +SKIP
    >>> encoding_models = np.random.default_rng(42).random(
    ...     (3, env.n_bins)
    ... )  # doctest: +SKIP
    >>> for i, shuffled in enumerate(  # doctest: +SKIP
    ...     shuffle_place_fields_circular_2d(encoding_models, env, n_shuffles=3, rng=42)
    ... ):
    ...     print(f"Shuffle {i}: shape={shuffled.shape}")
    Shuffle 0: shape=(3, 36)
    Shuffle 1: shape=(3, 36)
    Shuffle 2: shape=(3, 36)

    See Also
    --------
    shuffle_place_fields_circular : 1D circular shift
    shuffle_cell_identity : Shuffle neuron-to-place-field mapping
    """
    # Validate 2D environment
    if env.n_dims != 2:
        raise ValueError(
            f"shuffle_place_fields_circular_2d requires a 2D environment, "
            f"got {env.n_dims}D"
        )

    # Get grid shape from environment
    if not hasattr(env.layout, "grid_shape") or env.layout.grid_shape is None:
        raise ValueError(
            "shuffle_place_fields_circular_2d requires environment with grid layout "
            "(grid_shape attribute)"
        )

    grid_shape = env.layout.grid_shape

    # Verify encoding models match expected size
    expected_bins = int(np.prod(grid_shape))
    actual_bins = encoding_models.shape[1]
    if actual_bins != expected_bins:
        raise ValueError(
            f"shuffle_place_fields_circular_2d requires encoding models with "
            f"{expected_bins} bins (matching grid_shape {grid_shape}), but got "
            f"{actual_bins} bins. This may occur if the environment has inactive "
            f"bins (masked grid). Consider using shuffle_place_fields_circular "
            f"for masked environments."
        )
    generator = _ensure_rng(rng)
    n_neurons = encoding_models.shape[0]

    # Pre-compute indices for vectorized 2D roll
    # Create meshgrid of indices
    rows, cols = np.arange(grid_shape[0]), np.arange(grid_shape[1])

    for _ in range(n_shuffles):
        # Generate random shift amounts for each neuron in each dimension
        shifts_x = np.asarray(generator.integers(0, grid_shape[0], size=n_neurons))
        shifts_y = np.asarray(generator.integers(0, grid_shape[1], size=n_neurons))

        # Reshape all fields to 2D grids: (n_neurons, grid_shape[0], grid_shape[1])
        fields_2d = encoding_models.reshape(n_neurons, grid_shape[0], grid_shape[1])

        # Vectorized 2D circular shift using advanced indexing
        # For each neuron i, we want:
        #   row_indices = (rows - shifts_x[i]) % grid_shape[0]
        #   col_indices = (cols - shifts_y[i]) % grid_shape[1]
        # Shape: (n_neurons, grid_shape[0])
        shifted_row_indices = (rows - shifts_x[:, np.newaxis]) % grid_shape[0]
        # Shape: (n_neurons, grid_shape[1])
        shifted_col_indices = (cols - shifts_y[:, np.newaxis]) % grid_shape[1]

        # Apply shifts using advanced indexing
        # fields_2d[n, shifted_row_indices[n, :, None], shifted_col_indices[n, None, :]]
        neuron_indices = np.arange(n_neurons)[:, np.newaxis, np.newaxis]
        row_idx = shifted_row_indices[:, :, np.newaxis]
        col_idx = shifted_col_indices[:, np.newaxis, :]

        shifted_2d = fields_2d[neuron_indices, row_idx, col_idx]

        # Flatten back to 1D: (n_neurons, n_bins)
        shuffled = shifted_2d.reshape(n_neurons, -1)

        yield shuffled

shuffle_posterior_circular ¶

shuffle_posterior_circular(posterior: NDArray[float64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.float64], None, None]

Circularly shift posterior at each time bin independently.

Controls for chance linear alignment of position estimates by disrupting trajectory progression while preserving local smoothness. Each time bin's posterior is independently shifted by a random amount along the position axis.

Parameters:

Name	Type	Description	Default
`posterior`	`(NDArray[float64], shape(n_time_bins, n_bins))`	Posterior probability distribution from decoding. Each row should sum to 1.0.	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_posterior`	`(NDArray[float64], shape(n_time_bins, n_bins))`	Posterior with each row circularly shifted by a random amount.

Notes

Each time bin is shifted independently
Preserves the shape of each instantaneous posterior
Destroys temporal continuity of decoded positions
Caution: Can generate position representations that don't exist in original data (edge effects near track boundaries)
Normalization is preserved (each row still sums to 1.0)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_posterior_circular

>>> # Create a normalized posterior (3 time bins, 5 spatial bins)
>>> raw = np.array(
...     [
...         [0.1, 0.2, 0.4, 0.2, 0.1],
...         [0.5, 0.3, 0.1, 0.05, 0.05],
...         [0.05, 0.1, 0.2, 0.4, 0.25],
...     ]
... )
>>> for i, shuffled in enumerate(
...     shuffle_posterior_circular(raw, n_shuffles=3, rng=42)
... ):
...     print(
...         f"Shuffle {i}: shape={shuffled.shape}, sums to 1={np.allclose(shuffled.sum(axis=1), 1.0)}"
...     )
Shuffle 0: shape=(3, 5), sums to 1=True
Shuffle 1: shape=(3, 5), sums to 1=True
Shuffle 2: shape=(3, 5), sums to 1=True

See Also

shuffle_posterior_weighted_circular : Weighted circular shift with edge effect mitigation shuffle_place_fields_circular : Circular shift of place fields

Source code in src/neurospatial/stats/shuffle.py

def shuffle_posterior_circular(
    posterior: NDArray[np.float64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Circularly shift posterior at each time bin independently.

    Controls for chance linear alignment of position estimates by disrupting
    trajectory progression while preserving local smoothness. Each time bin's
    posterior is independently shifted by a random amount along the position
    axis.

    Parameters
    ----------
    posterior : NDArray[np.float64], shape (n_time_bins, n_bins)
        Posterior probability distribution from decoding. Each row should
        sum to 1.0.
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_posterior : NDArray[np.float64], shape (n_time_bins, n_bins)
        Posterior with each row circularly shifted by a random amount.

    Notes
    -----
    - Each time bin is shifted independently
    - Preserves the shape of each instantaneous posterior
    - Destroys temporal continuity of decoded positions
    - Caution: Can generate position representations that don't exist in
      original data (edge effects near track boundaries)
    - Normalization is preserved (each row still sums to 1.0)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_posterior_circular

    >>> # Create a normalized posterior (3 time bins, 5 spatial bins)
    >>> raw = np.array(
    ...     [
    ...         [0.1, 0.2, 0.4, 0.2, 0.1],
    ...         [0.5, 0.3, 0.1, 0.05, 0.05],
    ...         [0.05, 0.1, 0.2, 0.4, 0.25],
    ...     ]
    ... )
    >>> for i, shuffled in enumerate(
    ...     shuffle_posterior_circular(raw, n_shuffles=3, rng=42)
    ... ):
    ...     print(
    ...         f"Shuffle {i}: shape={shuffled.shape}, sums to 1={np.allclose(shuffled.sum(axis=1), 1.0)}"
    ...     )
    Shuffle 0: shape=(3, 5), sums to 1=True
    Shuffle 1: shape=(3, 5), sums to 1=True
    Shuffle 2: shape=(3, 5), sums to 1=True

    See Also
    --------
    shuffle_posterior_weighted_circular : Weighted circular shift with edge
        effect mitigation
    shuffle_place_fields_circular : Circular shift of place fields
    """
    generator = _ensure_rng(rng)
    n_time_bins, n_bins = posterior.shape

    # Pre-compute column indices for vectorized roll
    col_indices = np.arange(n_bins)

    for _ in range(n_shuffles):
        # Generate random shift amounts for each time bin
        shifts = np.asarray(generator.integers(0, n_bins, size=n_time_bins))

        # Vectorized circular shift using advanced indexing
        # For each row i, we want columns (col_indices - shifts[i]) % n_bins
        # Shape: (n_time_bins, n_bins)
        shifted_indices = (col_indices - shifts[:, np.newaxis]) % n_bins
        shuffled = posterior[np.arange(n_time_bins)[:, np.newaxis], shifted_indices]

        yield shuffled

shuffle_posterior_weighted_circular ¶

shuffle_posterior_weighted_circular(posterior: NDArray[float64], *, edge_buffer: int = 5, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.float64], None, None]

Weighted circular shift with edge effect mitigation.

Refined version of posterior shuffle that maintains non-uniformity and reduces edge effects by restricting shifts when the MAP position is near track boundaries. This is more conservative than standard circular shuffle.

Parameters:

Name	Type	Description	Default
`posterior`	`(NDArray[float64], shape(n_time_bins, n_bins))`	Posterior probability distribution from decoding. Each row should sum to 1.0.	required
`edge_buffer`	`int`	Number of bins from track ends where shifts are restricted. When the MAP position is within `edge_buffer` bins of either edge, the shift is restricted to keep the MAP position within bounds (not wrapping to the other end). Default justification (5 bins): For typical 1D track decoding with 2-5 cm spatial bins, 5 bins = 10-25 cm buffer region. This prevents unrealistic position jumps (e.g., rat at start wrapping to end) while allowing sufficient shift variability. Recommended values: Short tracks (<1m): edge_buffer = 3-5 bins Medium tracks (1-2m): edge_buffer = 5-10 bins Long tracks (>2m): edge_buffer = 10-15 bins Adaptive: edge_buffer = max(5, n_bins // 20) For 2D environments or circular tracks, use `shuffle_posterior_circular` instead (no edge restriction needed).	`5`
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_posterior`	`(NDArray[float64], shape(n_time_bins, n_bins))`	Posterior with weighted circular shifts applied.

Notes

More conservative than standard circular shuffle
Shifts are restricted near track ends to mitigate edge effects
When MAP position is within edge_buffer of an edge, shift is constrained to avoid wrapping probability mass to the other end
Preserves row normalization (each row still sums to 1.0)

The edge restriction works as follows: - For each time bin, compute MAP position (argmax) - If MAP is within edge_buffer of left edge (bin < edge_buffer): restrict shift to [-MAP, n_bins - 2*edge_buffer] to keep MAP in center - If MAP is within edge_buffer of right edge (bin >= n_bins - edge_buffer): restrict shift to [-(MAP - center), n_bins - MAP] to keep MAP in center - Otherwise: allow full circular shift [0, n_bins)

Example with n_bins=20 and edge_buffer=5: - MAP at bin 2 (near left): shift limited to [-2, 10) so MAP stays in [0, 11] - MAP at bin 17 (near right): shift limited to [-7, 3) so MAP stays in [10, 19] - MAP at bin 10 (center): full circular shift [0, 20)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_posterior_weighted_circular

>>> # Create a normalized posterior (3 time bins, 10 spatial bins)
>>> raw = np.random.default_rng(42).random((3, 10))
>>> posterior = raw / raw.sum(axis=1, keepdims=True)
>>> for i, shuffled in enumerate(
...     shuffle_posterior_weighted_circular(
...         posterior, edge_buffer=2, n_shuffles=3, rng=42
...     )
... ):
...     print(f"Shuffle {i}: sums to 1={np.allclose(shuffled.sum(axis=1), 1.0)}")
Shuffle 0: sums to 1=True
Shuffle 1: sums to 1=True
Shuffle 2: sums to 1=True

See Also

shuffle_posterior_circular : Standard circular shift without edge restriction shuffle_place_fields_circular : Circular shift of place fields

Source code in src/neurospatial/stats/shuffle.py

def shuffle_posterior_weighted_circular(
    posterior: NDArray[np.float64],
    *,
    edge_buffer: int = 5,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Weighted circular shift with edge effect mitigation.

    Refined version of posterior shuffle that maintains non-uniformity and
    reduces edge effects by restricting shifts when the MAP position is near
    track boundaries. This is more conservative than standard circular shuffle.

    Parameters
    ----------
    posterior : NDArray[np.float64], shape (n_time_bins, n_bins)
        Posterior probability distribution from decoding. Each row should
        sum to 1.0.
    edge_buffer : int, default=5
        Number of bins from track ends where shifts are restricted.
        When the MAP position is within ``edge_buffer`` bins of either edge,
        the shift is restricted to keep the MAP position within bounds
        (not wrapping to the other end).

        **Default justification (5 bins)**: For typical 1D track decoding with
        2-5 cm spatial bins, 5 bins = 10-25 cm buffer region. This prevents
        unrealistic position jumps (e.g., rat at start wrapping to end) while
        allowing sufficient shift variability. Recommended values:

        - Short tracks (<1m): edge_buffer = 3-5 bins
        - Medium tracks (1-2m): edge_buffer = 5-10 bins
        - Long tracks (>2m): edge_buffer = 10-15 bins
        - Adaptive: edge_buffer = max(5, n_bins // 20)

        For 2D environments or circular tracks, use ``shuffle_posterior_circular``
        instead (no edge restriction needed).
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_posterior : NDArray[np.float64], shape (n_time_bins, n_bins)
        Posterior with weighted circular shifts applied.

    Notes
    -----
    - More conservative than standard circular shuffle
    - Shifts are restricted near track ends to mitigate edge effects
    - When MAP position is within ``edge_buffer`` of an edge, shift is
      constrained to avoid wrapping probability mass to the other end
    - Preserves row normalization (each row still sums to 1.0)

    The edge restriction works as follows:
    - For each time bin, compute MAP position (argmax)
    - If MAP is within ``edge_buffer`` of left edge (bin < edge_buffer):
      restrict shift to [-MAP, n_bins - 2*edge_buffer] to keep MAP in center
    - If MAP is within ``edge_buffer`` of right edge (bin >= n_bins - edge_buffer):
      restrict shift to [-(MAP - center), n_bins - MAP] to keep MAP in center
    - Otherwise: allow full circular shift [0, n_bins)

    Example with n_bins=20 and edge_buffer=5:
    - MAP at bin 2 (near left): shift limited to [-2, 10) so MAP stays in [0, 11]
    - MAP at bin 17 (near right): shift limited to [-7, 3) so MAP stays in [10, 19]
    - MAP at bin 10 (center): full circular shift [0, 20)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_posterior_weighted_circular

    >>> # Create a normalized posterior (3 time bins, 10 spatial bins)
    >>> raw = np.random.default_rng(42).random((3, 10))
    >>> posterior = raw / raw.sum(axis=1, keepdims=True)
    >>> for i, shuffled in enumerate(
    ...     shuffle_posterior_weighted_circular(
    ...         posterior, edge_buffer=2, n_shuffles=3, rng=42
    ...     )
    ... ):
    ...     print(f"Shuffle {i}: sums to 1={np.allclose(shuffled.sum(axis=1), 1.0)}")
    Shuffle 0: sums to 1=True
    Shuffle 1: sums to 1=True
    Shuffle 2: sums to 1=True

    See Also
    --------
    shuffle_posterior_circular : Standard circular shift without edge restriction
    shuffle_place_fields_circular : Circular shift of place fields
    """
    generator = _ensure_rng(rng)
    n_time_bins, n_bins = posterior.shape

    # Handle empty posterior
    if n_time_bins == 0:
        for _ in range(n_shuffles):
            yield posterior.copy()
        return

    # Compute MAP positions for each time bin (used to determine shift restrictions)
    map_positions = np.argmax(posterior, axis=1)

    # NOTE: the per-time-bin loop below is intentionally not vectorized. Each
    # iteration draws exactly one integer from ``generator`` in time order,
    # and the *range* of that draw depends on the bin's MAP position (edge
    # bins use a narrower, position-dependent range). Drawing all shifts at
    # once would change the RNG consumption pattern and therefore the output
    # for a fixed seed, breaking the reproducible-shuffle contract that
    # callers rely on for null distributions. The cost is also bounded by
    # n_time_bins (a single candidate event, typically tens of bins), so
    # vectorizing buys little. Left as a loop deliberately.
    for _ in range(n_shuffles):
        shuffled = np.empty_like(posterior)
        for i in range(n_time_bins):
            map_pos = map_positions[i]

            # Determine allowed shift range based on MAP position
            # The key principle: restrict shifts to prevent probability mass from
            # wrapping to the opposite end of the track when near edges.
            if edge_buffer == 0 or n_bins <= 2 * edge_buffer:
                # No edge restriction or buffer spans entire track
                # Allow any shift in [0, n_bins)
                shift_amount = generator.integers(0, n_bins)
            elif map_pos < edge_buffer:
                # Near left edge: restrict shift range to avoid wrapping to far right
                # Allow negative shifts (move left) but limit positive shifts
                min_shift = -map_pos  # Don't shift MAP below 0
                # Limit positive shift: keep shifted MAP within center region
                max_shift = min(n_bins - map_pos, n_bins - 2 * edge_buffer)
                if max_shift <= min_shift:
                    max_shift = min_shift + 1  # Ensure valid range
                shift_amount = generator.integers(min_shift, max_shift)
            elif map_pos >= n_bins - edge_buffer:
                # Near right edge: restrict shift range to avoid wrapping to far left
                # Allow positive shifts (move right) but limit negative shifts
                max_shift = n_bins - map_pos  # Don't shift MAP beyond last bin
                # Limit negative shift: keep shifted MAP within center region
                min_shift = max(-(map_pos - (n_bins - 2 * edge_buffer)), -(map_pos))
                if max_shift <= min_shift:
                    min_shift = max_shift - 1  # Ensure valid range
                shift_amount = generator.integers(min_shift, max_shift)
            else:
                # Not near edge: allow full circular shift
                shift_amount = generator.integers(0, n_bins)

            shuffled[i, :] = np.roll(posterior[i, :], shift_amount)
        yield shuffled

shuffle_spikes_isi ¶

shuffle_spikes_isi(spike_times: NDArray[float64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.float64], None, None]

Shuffle inter-spike intervals to test ISI ordering significance.

Creates surrogate spike trains by randomly permuting the inter-spike intervals (ISIs). The first spike time is preserved, and ISIs are shuffled and cumulatively summed to create new spike times. This preserves the ISI distribution while destroying temporal correlations.

Parameters:

Name	Type	Description	Default
`spike_times`	`(NDArray[float64], shape(n_spikes))`	Sorted array of spike times for a single neuron.	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_times`	`(NDArray[float64], shape(n_spikes))`	Spike times with shuffled ISIs (same first spike time).

Notes

Preserves the first spike time exactly
Preserves the ISI distribution (same ISI values, different order)
Destroys temporal correlations between consecutive ISIs
Output spike times are always sorted (monotonically increasing)
Useful for testing ISI-dependent effects (e.g., burst detection, adaptation)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_spikes_isi

>>> spike_times = np.array([0.1, 0.15, 0.25, 0.4, 0.45, 0.7])
>>> for i, shuffled in enumerate(
...     shuffle_spikes_isi(spike_times, n_shuffles=3, rng=42)
... ):
...     print(
...         f"Shuffle {i}: first spike preserved = {shuffled[0] == spike_times[0]}"
...     )
Shuffle 0: first spike preserved = True
Shuffle 1: first spike preserved = True
Shuffle 2: first spike preserved = True

See Also

shuffle_time_bins : Shuffle temporal order of binned spike counts. generate_poisson_surrogates : Generate Poisson surrogate spike trains.

Source code in src/neurospatial/stats/shuffle.py

def shuffle_spikes_isi(
    spike_times: NDArray[np.float64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Shuffle inter-spike intervals to test ISI ordering significance.

    Creates surrogate spike trains by randomly permuting the inter-spike
    intervals (ISIs). The first spike time is preserved, and ISIs are
    shuffled and cumulatively summed to create new spike times. This
    preserves the ISI distribution while destroying temporal correlations.

    Parameters
    ----------
    spike_times : NDArray[np.float64], shape (n_spikes,)
        Sorted array of spike times for a single neuron.
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_times : NDArray[np.float64], shape (n_spikes,)
        Spike times with shuffled ISIs (same first spike time).

    Notes
    -----
    - Preserves the first spike time exactly
    - Preserves the ISI distribution (same ISI values, different order)
    - Destroys temporal correlations between consecutive ISIs
    - Output spike times are always sorted (monotonically increasing)
    - Useful for testing ISI-dependent effects (e.g., burst detection,
      adaptation)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_spikes_isi

    >>> spike_times = np.array([0.1, 0.15, 0.25, 0.4, 0.45, 0.7])
    >>> for i, shuffled in enumerate(
    ...     shuffle_spikes_isi(spike_times, n_shuffles=3, rng=42)
    ... ):
    ...     print(
    ...         f"Shuffle {i}: first spike preserved = {shuffled[0] == spike_times[0]}"
    ...     )
    Shuffle 0: first spike preserved = True
    Shuffle 1: first spike preserved = True
    Shuffle 2: first spike preserved = True

    See Also
    --------
    shuffle_time_bins : Shuffle temporal order of binned spike counts.
    generate_poisson_surrogates : Generate Poisson surrogate spike trains.
    """
    generator = _ensure_rng(rng)

    # Handle edge cases
    if len(spike_times) <= 1:
        for _ in range(n_shuffles):
            yield spike_times.copy()
        return

    # Compute inter-spike intervals
    isis = np.diff(spike_times)
    first_spike = spike_times[0]

    for _ in range(n_shuffles):
        # Shuffle the ISIs
        shuffled_isis = generator.permutation(isis)
        # Reconstruct spike times from shuffled ISIs
        shuffled_times = np.concatenate(
            [[first_spike], first_spike + np.cumsum(shuffled_isis)]
        )
        yield shuffled_times

shuffle_time_bins ¶

shuffle_time_bins(spike_counts: NDArray[int64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.int64], None, None]

Shuffle temporal order of time bins within an event.

Primary test for sequential structure. Disrupts temporal order while preserving instantaneous firing characteristics.

This shuffle randomly permutes the rows (time bins) of the spike count matrix. Each row is kept intact, so the spike counts per neuron per time bin are preserved, but the temporal sequence is destroyed.

Parameters:

Name	Type	Description	Default
`spike_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts for a single candidate event (PBE/SWR).	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts with time bins in random order.

Notes

Preserves spike counts per neuron per time bin
Preserves total spikes per neuron across event
Destroys temporal sequence information
Underpowered with very few time bins (<5)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_time_bins

>>> spike_counts = np.array([[0, 1], [2, 0], [1, 1]], dtype=np.int64)
>>> for i, shuffled in enumerate(
...     shuffle_time_bins(spike_counts, n_shuffles=3, rng=42)
... ):
...     print(f"Shuffle {i}: shape={shuffled.shape}, sum={shuffled.sum()}")
Shuffle 0: shape=(3, 2), sum=5
Shuffle 1: shape=(3, 2), sum=5
Shuffle 2: shape=(3, 2), sum=5

See Also

shuffle_time_bins_coherent : Coherent shuffle (same permutation for all neurons)

Source code in src/neurospatial/stats/shuffle.py

def shuffle_time_bins(
    spike_counts: NDArray[np.int64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.int64], None, None]:
    """Shuffle temporal order of time bins within an event.

    **Primary test for sequential structure.** Disrupts temporal order while
    preserving instantaneous firing characteristics.

    This shuffle randomly permutes the rows (time bins) of the spike count
    matrix. Each row is kept intact, so the spike counts per neuron per
    time bin are preserved, but the temporal sequence is destroyed.

    Parameters
    ----------
    spike_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts for a single candidate event (PBE/SWR).
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts with time bins in random order.

    Notes
    -----
    - Preserves spike counts per neuron per time bin
    - Preserves total spikes per neuron across event
    - Destroys temporal sequence information
    - Underpowered with very few time bins (<5)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_time_bins

    >>> spike_counts = np.array([[0, 1], [2, 0], [1, 1]], dtype=np.int64)
    >>> for i, shuffled in enumerate(
    ...     shuffle_time_bins(spike_counts, n_shuffles=3, rng=42)
    ... ):
    ...     print(f"Shuffle {i}: shape={shuffled.shape}, sum={shuffled.sum()}")
    Shuffle 0: shape=(3, 2), sum=5
    Shuffle 1: shape=(3, 2), sum=5
    Shuffle 2: shape=(3, 2), sum=5

    See Also
    --------
    shuffle_time_bins_coherent : Coherent shuffle (same permutation for all neurons)
    """
    generator = _ensure_rng(rng)
    n_time_bins = spike_counts.shape[0]

    for _ in range(n_shuffles):
        # Generate a random permutation of row indices
        perm = generator.permutation(n_time_bins)
        # Apply permutation to rows
        yield spike_counts[perm].copy()

shuffle_time_bins_coherent ¶

shuffle_time_bins_coherent(spike_counts: NDArray[int64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.int64], None, None]

Shuffle time bins coherently across all neurons (time-swap shuffle).

Preserves population co-activation structure while disrupting temporal order. This is the same as shuffle_time_bins in implementation, but is named differently to emphasize that the permutation is coherent across all neurons (preserving instantaneous population vectors).

This shuffle is less conservative than independent per-neuron shuffles because pairwise correlations between neurons are maintained within each time bin.

Parameters:

Name	Type	Description	Default
`spike_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts for a single candidate event.	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Spike counts with rows (time bins) permuted coherently.

Notes

All neurons see the same temporal permutation
Preserves instantaneous population vectors
Destroys temporal progression but keeps co-firing structure
Mathematically equivalent to shuffle_time_bins, but the function name clarifies the intent (testing temporal structure while preserving population co-activation)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_time_bins_coherent

>>> spike_counts = np.array([[0, 1, 2], [2, 0, 1], [1, 2, 0]], dtype=np.int64)
>>> for shuffled in shuffle_time_bins_coherent(spike_counts, n_shuffles=2, rng=42):
...     # Each row is preserved exactly (coherent permutation)
...     original_rows = {tuple(row) for row in spike_counts}
...     shuffled_rows = {tuple(row) for row in shuffled}
...     print(f"Row sets match: {original_rows == shuffled_rows}")
Row sets match: True
Row sets match: True

See Also

shuffle_time_bins : Primary temporal shuffle (same implementation) shuffle_cell_identity : Shuffle neuron-to-place-field mapping

Source code in src/neurospatial/stats/shuffle.py

def shuffle_time_bins_coherent(
    spike_counts: NDArray[np.int64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.int64], None, None]:
    """Shuffle time bins coherently across all neurons (time-swap shuffle).

    Preserves population co-activation structure while disrupting temporal
    order. This is the same as ``shuffle_time_bins`` in implementation, but
    is named differently to emphasize that the permutation is coherent
    across all neurons (preserving instantaneous population vectors).

    This shuffle is **less conservative** than independent per-neuron shuffles
    because pairwise correlations between neurons are maintained within each
    time bin.

    Parameters
    ----------
    spike_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts for a single candidate event.
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Spike counts with rows (time bins) permuted coherently.

    Notes
    -----
    - All neurons see the same temporal permutation
    - Preserves instantaneous population vectors
    - Destroys temporal progression but keeps co-firing structure
    - Mathematically equivalent to ``shuffle_time_bins``, but the function
      name clarifies the intent (testing temporal structure while preserving
      population co-activation)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_time_bins_coherent

    >>> spike_counts = np.array([[0, 1, 2], [2, 0, 1], [1, 2, 0]], dtype=np.int64)
    >>> for shuffled in shuffle_time_bins_coherent(spike_counts, n_shuffles=2, rng=42):
    ...     # Each row is preserved exactly (coherent permutation)
    ...     original_rows = {tuple(row) for row in spike_counts}
    ...     shuffled_rows = {tuple(row) for row in shuffled}
    ...     print(f"Row sets match: {original_rows == shuffled_rows}")
    Row sets match: True
    Row sets match: True

    See Also
    --------
    shuffle_time_bins : Primary temporal shuffle (same implementation)
    shuffle_cell_identity : Shuffle neuron-to-place-field mapping
    """
    # This is an exact alias of shuffle_time_bins: permuting whole rows
    # already applies one coherent permutation across all neurons. We
    # delegate rather than duplicate the body so the two cannot silently
    # diverge; the separate public name documents the intent (testing
    # temporal structure while preserving population co-activation).
    yield from shuffle_time_bins(spike_counts, n_shuffles=n_shuffles, rng=rng)

shuffle_trials ¶

shuffle_trials(trial_labels: NDArray[int64], *, n_shuffles: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.int64], None, None]

Shuffle trial labels to test trial identity significance.

Randomly permutes trial labels to create a null distribution for testing whether observed effects depend on trial identity. Useful for testing learning effects, trial-type differences, or temporal dependencies.

Parameters:

Name	Type	Description	Default
`trial_labels`	`(NDArray[int64], shape(n_samples))`	Array of trial labels/indices for each sample (e.g., timepoint). Labels can be any integers (e.g., 0, 1, 2 or trial IDs).	required
`n_shuffles`	`int`	Number of shuffled versions to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`shuffled_labels`	`(NDArray[int64], shape(n_samples))`	Trial labels with permuted order.

Notes

Permutes the entire label array (all samples)
Preserves the distribution of labels (same counts per label)
Destroys the temporal association between samples and trials
Useful for testing trial-type effects or learning-related changes

Examples:

>>> import numpy as np
>>> from neurospatial.stats.shuffle import shuffle_trials

>>> trial_labels = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2])
>>> for i, shuffled in enumerate(
...     shuffle_trials(trial_labels, n_shuffles=3, rng=42)
... ):
...     print(
...         f"Shuffle {i}: unique values preserved = {set(shuffled) == set(trial_labels)}"
...     )
Shuffle 0: unique values preserved = True
Shuffle 1: unique values preserved = True
Shuffle 2: unique values preserved = True

See Also

shuffle_time_bins : Shuffle temporal order of spike counts. shuffle_cell_identity : Shuffle neuron-to-place-field mapping.

Source code in src/neurospatial/stats/shuffle.py

def shuffle_trials(
    trial_labels: NDArray[np.int64],
    *,
    n_shuffles: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.int64], None, None]:
    """Shuffle trial labels to test trial identity significance.

    Randomly permutes trial labels to create a null distribution for testing
    whether observed effects depend on trial identity. Useful for testing
    learning effects, trial-type differences, or temporal dependencies.

    Parameters
    ----------
    trial_labels : NDArray[np.int64], shape (n_samples,)
        Array of trial labels/indices for each sample (e.g., timepoint).
        Labels can be any integers (e.g., 0, 1, 2 or trial IDs).
    n_shuffles : int, default=1000
        Number of shuffled versions to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    shuffled_labels : NDArray[np.int64], shape (n_samples,)
        Trial labels with permuted order.

    Notes
    -----
    - Permutes the entire label array (all samples)
    - Preserves the distribution of labels (same counts per label)
    - Destroys the temporal association between samples and trials
    - Useful for testing trial-type effects or learning-related changes

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.shuffle import shuffle_trials

    >>> trial_labels = np.array([0, 0, 0, 1, 1, 1, 2, 2, 2])
    >>> for i, shuffled in enumerate(
    ...     shuffle_trials(trial_labels, n_shuffles=3, rng=42)
    ... ):
    ...     print(
    ...         f"Shuffle {i}: unique values preserved = {set(shuffled) == set(trial_labels)}"
    ...     )
    Shuffle 0: unique values preserved = True
    Shuffle 1: unique values preserved = True
    Shuffle 2: unique values preserved = True

    See Also
    --------
    shuffle_time_bins : Shuffle temporal order of spike counts.
    shuffle_cell_identity : Shuffle neuron-to-place-field mapping.
    """
    generator = _ensure_rng(rng)
    n_samples = len(trial_labels)

    for _ in range(n_shuffles):
        # Generate a random permutation of all sample indices
        perm = generator.permutation(n_samples)
        yield trial_labels[perm].copy()

generate_inhomogeneous_poisson_surrogates ¶

generate_inhomogeneous_poisson_surrogates(spike_counts: NDArray[int64], *, smoothing_window: int = 3, n_surrogates: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.int64], None, None]

Generate inhomogeneous Poisson surrogate spike trains with smoothed rates.

Creates surrogate spike counts by sampling from Poisson distributions with time-varying rates estimated from smoothed spike counts. This preserves slow rate fluctuations while destroying fine temporal structure.

Tests that structure exceeds time-varying rate expectations. More conservative than homogeneous surrogates - controls for rate modulation while testing for sequence-specific coding.

Parameters:

Name	Type	Description	Default
`spike_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Original spike counts per neuron per time bin.	required
`smoothing_window`	`int`	Size of the uniform smoothing window (in time bins) applied to estimate time-varying rates. Larger values preserve less temporal detail (closer to homogeneous). Must be at least 1.	`3`
`n_surrogates`	`int`	Number of surrogate spike trains to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`surrogate`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Surrogate spike counts sampled from time-varying Poisson distributions.

Notes

Uses uniform 1D filter to smooth spike counts in time
Smoothed counts are used as time-varying Poisson rate parameters
Preserves slow rate fluctuations (e.g., event onsets)
Destroys fine temporal correlations and millisecond-scale patterns
More appropriate than homogeneous surrogates when rate modulation is a confound (e.g., ramping activity during replay events)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.surrogates import (
...     generate_inhomogeneous_poisson_surrogates,
... )

>>> spike_counts = np.array(
...     [[0, 1], [2, 2], [3, 3], [2, 2], [1, 1]], dtype=np.int64
... )
>>> for i, surrogate in enumerate(
...     generate_inhomogeneous_poisson_surrogates(
...         spike_counts, smoothing_window=3, n_surrogates=3, rng=42
...     )
... ):
...     print(f"Surrogate {i}: shape={surrogate.shape}")
Surrogate 0: shape=(5, 2)
Surrogate 1: shape=(5, 2)
Surrogate 2: shape=(5, 2)

See Also

generate_poisson_surrogates : Homogeneous Poisson surrogates generate_jittered_spikes : Temporal jitter surrogates

Source code in src/neurospatial/stats/surrogates.py

def generate_inhomogeneous_poisson_surrogates(
    spike_counts: NDArray[np.int64],
    *,
    smoothing_window: int = 3,
    n_surrogates: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.int64], None, None]:
    """Generate inhomogeneous Poisson surrogate spike trains with smoothed rates.

    Creates surrogate spike counts by sampling from Poisson distributions
    with time-varying rates estimated from smoothed spike counts. This
    preserves slow rate fluctuations while destroying fine temporal structure.

    **Tests that structure exceeds time-varying rate expectations.** More
    conservative than homogeneous surrogates - controls for rate modulation
    while testing for sequence-specific coding.

    Parameters
    ----------
    spike_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Original spike counts per neuron per time bin.
    smoothing_window : int, default=3
        Size of the uniform smoothing window (in time bins) applied to
        estimate time-varying rates. Larger values preserve less temporal
        detail (closer to homogeneous). Must be at least 1.
    n_surrogates : int, default=1000
        Number of surrogate spike trains to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    surrogate : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Surrogate spike counts sampled from time-varying Poisson distributions.

    Notes
    -----
    - Uses uniform 1D filter to smooth spike counts in time
    - Smoothed counts are used as time-varying Poisson rate parameters
    - Preserves slow rate fluctuations (e.g., event onsets)
    - Destroys fine temporal correlations and millisecond-scale patterns
    - More appropriate than homogeneous surrogates when rate modulation
      is a confound (e.g., ramping activity during replay events)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.surrogates import (
    ...     generate_inhomogeneous_poisson_surrogates,
    ... )

    >>> spike_counts = np.array(
    ...     [[0, 1], [2, 2], [3, 3], [2, 2], [1, 1]], dtype=np.int64
    ... )
    >>> for i, surrogate in enumerate(
    ...     generate_inhomogeneous_poisson_surrogates(
    ...         spike_counts, smoothing_window=3, n_surrogates=3, rng=42
    ...     )
    ... ):
    ...     print(f"Surrogate {i}: shape={surrogate.shape}")
    Surrogate 0: shape=(5, 2)
    Surrogate 1: shape=(5, 2)
    Surrogate 2: shape=(5, 2)

    See Also
    --------
    generate_poisson_surrogates : Homogeneous Poisson surrogates
    generate_jittered_spikes : Temporal jitter surrogates
    """
    from scipy.ndimage import uniform_filter1d

    generator = _ensure_rng(rng)
    n_time_bins, n_neurons = spike_counts.shape

    # Handle empty spike counts
    if n_time_bins == 0:
        for _ in range(n_surrogates):
            yield np.zeros((0, n_neurons), dtype=np.int64)
        return

    # Compute smoothed rates (time-varying lambda for Poisson)
    # Use uniform filter along time axis with 'nearest' mode for edge handling
    smoothed_counts = uniform_filter1d(
        spike_counts.astype(np.float64),
        size=smoothing_window,
        axis=0,
        mode="nearest",
    )

    # Ensure non-negative rates (should already be, but be safe)
    smoothed_counts = np.maximum(smoothed_counts, 0.0)

    for _ in range(n_surrogates):
        # Generate surrogate by sampling from Poisson with time-varying rates
        # Each (time_bin, neuron) pair uses its corresponding smoothed rate
        surrogate = generator.poisson(lam=smoothed_counts)
        yield surrogate.astype(np.int64)

generate_jittered_spikes ¶

generate_jittered_spikes(spike_times: NDArray[float64], jitter_std: float, *, n_surrogates: int = 1000, rng: Generator | int | None = None, window: tuple[float, float] | None = None) -> Generator[NDArray[np.float64], None, None]

Generate temporally jittered surrogate spike trains.

Creates surrogate spike trains by adding Gaussian noise to each spike time. This tests whether precise spike timing is significant while preserving the approximate temporal distribution of spikes.

Tests that precise spike timing matters. If observed effects depend on exact spike times, jittered surrogates should produce weaker effects.

Parameters:

Name	Type	Description	Default
`spike_times`	`(NDArray[float64], shape(n_spikes))`	Sorted array of spike times for a single neuron.	required
`jitter_std`	`float`	Standard deviation of the Gaussian jitter added to each spike time. Units should match spike_times (typically seconds).	required
`n_surrogates`	`int`	Number of surrogate spike trains to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`
`window`	`tuple[float, float] \| None`	If provided, clip spike times to stay within (min_time, max_time). Useful to prevent spikes from jittering outside the recording period.	`None`

Yields:

Name	Type	Description
`jittered_times`	`(NDArray[float64], shape(n_spikes))`	Spike times with Gaussian jitter applied (sorted).

Notes

Each spike is independently jittered by adding Gaussian noise
Output spike times are always sorted (monotonically increasing)
Preserves the number of spikes exactly
Preserves approximate temporal distribution at scales > jitter_std
Destroys precise timing relationships at scales < jitter_std
Useful for testing phase locking, synchrony, or sequence timing

Examples:

>>> import numpy as np
>>> from neurospatial.stats.surrogates import generate_jittered_spikes

>>> spike_times = np.array([0.1, 0.15, 0.25, 0.4, 0.45, 0.7])
>>> for i, jittered in enumerate(
...     generate_jittered_spikes(
...         spike_times, jitter_std=0.01, n_surrogates=3, rng=42
...     )
... ):
...     print(f"Surrogate {i}: n_spikes={len(jittered)}")
Surrogate 0: n_spikes=6
Surrogate 1: n_spikes=6
Surrogate 2: n_spikes=6

See Also

generate_poisson_surrogates : Rate-matched surrogates (destroys all timing) neurospatial.stats.shuffle.shuffle_spikes_isi : ISI-preserving shuffle

Source code in src/neurospatial/stats/surrogates.py

def generate_jittered_spikes(
    spike_times: NDArray[np.float64],
    jitter_std: float,
    *,
    n_surrogates: int = 1000,
    rng: np.random.Generator | int | None = None,
    window: tuple[float, float] | None = None,
) -> Generator[NDArray[np.float64], None, None]:
    """Generate temporally jittered surrogate spike trains.

    Creates surrogate spike trains by adding Gaussian noise to each spike
    time. This tests whether precise spike timing is significant while
    preserving the approximate temporal distribution of spikes.

    **Tests that precise spike timing matters.** If observed effects depend
    on exact spike times, jittered surrogates should produce weaker effects.

    Parameters
    ----------
    spike_times : NDArray[np.float64], shape (n_spikes,)
        Sorted array of spike times for a single neuron.
    jitter_std : float
        Standard deviation of the Gaussian jitter added to each spike time.
        Units should match spike_times (typically seconds).
    n_surrogates : int, default=1000
        Number of surrogate spike trains to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)
    window : tuple[float, float] | None, default=None
        If provided, clip spike times to stay within (min_time, max_time).
        Useful to prevent spikes from jittering outside the recording period.

    Yields
    ------
    jittered_times : NDArray[np.float64], shape (n_spikes,)
        Spike times with Gaussian jitter applied (sorted).

    Notes
    -----
    - Each spike is independently jittered by adding Gaussian noise
    - Output spike times are always sorted (monotonically increasing)
    - Preserves the number of spikes exactly
    - Preserves approximate temporal distribution at scales > jitter_std
    - Destroys precise timing relationships at scales < jitter_std
    - Useful for testing phase locking, synchrony, or sequence timing

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.surrogates import generate_jittered_spikes

    >>> spike_times = np.array([0.1, 0.15, 0.25, 0.4, 0.45, 0.7])
    >>> for i, jittered in enumerate(
    ...     generate_jittered_spikes(
    ...         spike_times, jitter_std=0.01, n_surrogates=3, rng=42
    ...     )
    ... ):
    ...     print(f"Surrogate {i}: n_spikes={len(jittered)}")
    Surrogate 0: n_spikes=6
    Surrogate 1: n_spikes=6
    Surrogate 2: n_spikes=6

    See Also
    --------
    generate_poisson_surrogates : Rate-matched surrogates (destroys all timing)
    neurospatial.stats.shuffle.shuffle_spikes_isi : ISI-preserving shuffle
    """
    generator = _ensure_rng(rng)
    n_spikes = len(spike_times)

    # Handle empty spike array
    if n_spikes == 0:
        for _ in range(n_surrogates):
            yield np.array([], dtype=np.float64)
        return

    for _ in range(n_surrogates):
        # Add Gaussian jitter to each spike time
        jitter = generator.normal(0.0, jitter_std, size=n_spikes)
        jittered_times = spike_times + jitter

        # Clip to window if specified
        if window is not None:
            jittered_times = np.clip(jittered_times, window[0], window[1])

        # Sort to maintain monotonicity
        jittered_times = np.sort(jittered_times)

        yield jittered_times

generate_poisson_surrogates ¶

generate_poisson_surrogates(spike_counts: NDArray[int64], *, n_surrogates: int = 1000, rng: Generator | int | None = None) -> Generator[NDArray[np.int64], None, None]

Generate homogeneous Poisson surrogate spike trains.

Creates surrogate spike counts by sampling from Poisson distributions with rates equal to the mean firing rate of each neuron across all time bins. This destroys all temporal structure while preserving average rates.

Tests that observed structure exceeds rate-based expectations. If sequential decoding patterns are significant, they should score higher than patterns from rate-matched surrogates.

Parameters:

Name	Type	Description	Default
`spike_counts`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Original spike counts per neuron per time bin.	required
`n_surrogates`	`int`	Number of surrogate spike trains to generate.	`1000`
`rng`	`Generator \| int \| None`	Random number generator for reproducibility. If Generator: Use directly If int: Seed for `np.random.default_rng()` If None: Use default RNG (not reproducible)	`None`

Yields:

Name	Type	Description
`surrogate`	`(NDArray[int64], shape(n_time_bins, n_neurons))`	Surrogate spike counts sampled from Poisson distributions.

Notes

Computes mean spike count per neuron across all time bins
Each time bin independently samples from Poisson(mean_count_per_neuron)
Destroys all temporal correlations and patterns
Preserves mean firing rates per neuron (statistically)
More conservative than shuffle methods (generates truly independent counts)

Examples:

>>> import numpy as np
>>> from neurospatial.stats.surrogates import generate_poisson_surrogates

>>> spike_counts = np.array(
...     [[0, 1], [2, 0], [1, 1]], dtype=np.int64
... )
>>> for i, surrogate in enumerate(
...     generate_poisson_surrogates(spike_counts, n_surrogates=3, rng=42)
... ):
...     print(f"Surrogate {i}: shape={surrogate.shape}, total={surrogate.sum()}")
Surrogate 0: shape=(3, 2), total=3
Surrogate 1: shape=(3, 2), total=5
Surrogate 2: shape=(3, 2), total=6

See Also

generate_inhomogeneous_poisson_surrogates : Time-varying rate surrogates generate_jittered_spikes : Temporal jitter surrogates

Source code in src/neurospatial/stats/surrogates.py

def generate_poisson_surrogates(
    spike_counts: NDArray[np.int64],
    *,
    n_surrogates: int = 1000,
    rng: np.random.Generator | int | None = None,
) -> Generator[NDArray[np.int64], None, None]:
    """Generate homogeneous Poisson surrogate spike trains.

    Creates surrogate spike counts by sampling from Poisson distributions
    with rates equal to the mean firing rate of each neuron across all time
    bins. This destroys all temporal structure while preserving average rates.

    **Tests that observed structure exceeds rate-based expectations.** If
    sequential decoding patterns are significant, they should score higher
    than patterns from rate-matched surrogates.

    Parameters
    ----------
    spike_counts : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Original spike counts per neuron per time bin.
    n_surrogates : int, default=1000
        Number of surrogate spike trains to generate.
    rng : np.random.Generator | int | None, default=None
        Random number generator for reproducibility.

        - If Generator: Use directly
        - If int: Seed for ``np.random.default_rng()``
        - If None: Use default RNG (not reproducible)

    Yields
    ------
    surrogate : NDArray[np.int64], shape (n_time_bins, n_neurons)
        Surrogate spike counts sampled from Poisson distributions.

    Notes
    -----
    - Computes mean spike count per neuron across all time bins
    - Each time bin independently samples from Poisson(mean_count_per_neuron)
    - Destroys all temporal correlations and patterns
    - Preserves mean firing rates per neuron (statistically)
    - More conservative than shuffle methods (generates truly independent counts)

    Examples
    --------
    >>> import numpy as np
    >>> from neurospatial.stats.surrogates import generate_poisson_surrogates

    >>> spike_counts = np.array(
    ...     [[0, 1], [2, 0], [1, 1]], dtype=np.int64
    ... )  # doctest: +SKIP
    >>> for i, surrogate in enumerate(  # doctest: +SKIP
    ...     generate_poisson_surrogates(spike_counts, n_surrogates=3, rng=42)
    ... ):
    ...     print(f"Surrogate {i}: shape={surrogate.shape}, total={surrogate.sum()}")
    Surrogate 0: shape=(3, 2), total=3
    Surrogate 1: shape=(3, 2), total=5
    Surrogate 2: shape=(3, 2), total=6

    See Also
    --------
    generate_inhomogeneous_poisson_surrogates : Time-varying rate surrogates
    generate_jittered_spikes : Temporal jitter surrogates
    """
    generator = _ensure_rng(rng)
    n_time_bins, n_neurons = spike_counts.shape

    # Handle empty spike counts
    if n_time_bins == 0:
        for _ in range(n_surrogates):
            yield np.zeros((0, n_neurons), dtype=np.int64)
        return

    # Compute mean spike count per neuron (averaged across time bins)
    # This is the lambda parameter for our Poisson distribution
    mean_counts = spike_counts.mean(axis=0)  # Shape: (n_neurons,)

    for _ in range(n_surrogates):
        # Generate surrogate by sampling from Poisson with mean rate
        # Each (time_bin, neuron) pair is sampled independently
        surrogate = generator.poisson(lam=mean_counts, size=(n_time_bins, n_neurons))
        yield surrogate.astype(np.int64)

neurospatial.stats¶

stats ¶

Classes¶

CircularBasisResult dataclass ¶

Attributes¶

design_matrix property ¶

ShuffleTestResult dataclass ¶

Attributes¶

is_significant property ¶

Functions¶

plot ¶

Functions¶

circular_basis ¶

circular_basis_metrics ¶

circular_circular_correlation ¶

circular_linear_correlation ¶

circular_mean ¶

circular_variance ¶

is_modulated ¶

mean_resultant_length ¶

phase_position_correlation ¶

plot_circular_basis_tuning ¶

rayleigh_test ¶

wrap_angle ¶

compute_shuffle_pvalue ¶

compute_shuffle_zscore ¶

shuffle_cell_identity ¶

shuffle_place_fields_circular ¶

shuffle_place_fields_circular_2d ¶

shuffle_posterior_circular ¶

shuffle_posterior_weighted_circular ¶

shuffle_spikes_isi ¶

shuffle_time_bins ¶

shuffle_time_bins_coherent ¶

shuffle_trials ¶

generate_inhomogeneous_poisson_surrogates ¶

generate_jittered_spikes ¶

generate_poisson_surrogates ¶

`neurospatial.stats`¶

CircularBasisResult `dataclass` ¶

design_matrix `property` ¶

ShuffleTestResult `dataclass` ¶

is_significant `property` ¶