inference¶

GPU-accelerated statistical inference for neuroimaging.

This module provides fast permutation testing and bootstrap resampling using optional GPU acceleration via PyTorch. When GPU is unavailable, efficiently uses CPU parallelization.

Inspired by BROCCOLI’s GPU permutation testing (Eklund et al. 2014).

Classes:

Methods:

Modules:

Examples:

>>> import numpy as np
>>> from nltools.algorithms.inference import one_sample_permutation_test

>>> # Simple one-sample test
>>> data = np.random.randn(30)  # 30 subjects
>>> result = one_sample_permutation_test(data, n_permute=5000)
>>> print(f"p-value: {result['p']:.3f}")

>>> # Voxel-wise test with GPU acceleration
>>> data = np.random.randn(30, 50000)  # 30 subjects, 50K voxels
>>> result = one_sample_permutation_test(data, n_permute=10000, device='gpu')
>>> print(f"Significant voxels: {(result['p'] < 0.05).sum()}")

Classes¶

OnlineBootstrapStats¶

OnlineBootstrapStats(shape: tuple[int, ...], save_samples: bool = False, percentiles: tuple[float, float] = (2.5, 97.5))

Memory-efficient online statistics aggregator for bootstrap samples.

Uses Welford’s algorithm for numerically stable online computation of mean and variance. Optionally stores all samples for exact percentile CIs.

Parameters:

Attributes:

Methods:

Examples:

>>> stats = OnlineBootstrapStats(shape=(100,), save_samples=False)
>>> for i in range(1000):
...     sample = np.random.randn(100)
...     stats.update(sample)
>>> results = stats.get_results()
>>> print(results.keys())
dict_keys(['mean', 'std', 'Z', 'p', 'ci_lower', 'ci_upper'])

Methods¶

get_results¶

get_results(tail: int | str = 2) -> dict[str, np.ndarray]

Compute final bootstrap statistics.

Parameters:

Returns:

Examples:

stats = OnlineBootstrapStats(shape=(100,), save_samples=False)
for _ in range(1000):
    stats.update(np.random.randn(100))
results = stats.get_results()
# results.keys() -> mean, std, Z, p, ci_lower, ci_upper

update¶

update(sample: np.ndarray) -> None

Update statistics with a new bootstrap sample.

Uses Welford’s algorithm for numerical stability.

Parameters:

Methods¶

circle_shift¶

circle_shift(data: np.ndarray, shift_amount: int | np.ndarray | None = None, random_state: int | np.random.RandomState | None = None) -> np.ndarray

Circular shift for time-series data.

Performs a circular shift that preserves autocorrelation structure. Useful for permutation tests on autocorrelated time series (e.g., fMRI). For 1D data, shifts by a single amount. For 2D data, shifts each feature (column) independently.

Parameters:

Returns:

Examples:

>>> x = np.array([1, 2, 3, 4, 5])
>>> circle_shift(x, shift_amount=2)
array([4, 5, 1, 2, 3])

>>> X = np.array([[1, 10], [2, 20], [3, 30], [4, 40]])
>>> circle_shift(X, shift_amount=np.array([1, 2]))
array([[ 4, 30],
       [ 1, 40],
       [ 2, 10],
       [ 3, 20]])

correlation_permutation_test¶

correlation_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, metric: str = 'pearson', tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

Correlation permutation test.

Tests whether the correlation between data1 and data2 is significantly different from zero by randomly permuting data1 and computing correlations.

Assumption: Observations are independent (i.i.d.). For autocorrelated time series, use timeseries_correlation_permutation_test with circle_shift or phase_randomize methods instead.

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> x = np.random.randn(100)
>>> y = x + np.random.randn(100) * 0.5  # Correlated
>>> result = correlation_permutation_test(x, y, n_permute=5000)
>>> result['correlation']
0.85
>>> result['p']
0.001

>>> # Multi-feature (2D arrays)
>>> data1 = np.random.randn(100, 10)  # 100 samples, 10 features
>>> data2 = data1 + np.random.randn(100, 10) * 0.3  # Correlated
>>> result = correlation_permutation_test(data1, data2, n_permute=5000)
>>> result['correlation'].shape
(10,)
>>> result['p'].shape
(10,)

>>> # GPU acceleration
>>> result = correlation_permutation_test(data1, data2, n_permute=5000, device='gpu')

distance_correlation¶

distance_correlation(x: np.ndarray, y: np.ndarray, bias_corrected: bool = True, ttest: bool = False) -> dict

Compute the distance correlation between 2 arrays to test for multivariate dependence (linear or non-linear).

Arrays must match on their first dimension. It’s almost always preferable to compute the bias_corrected version which can also optionally perform a ttest. This ttest operates on a statistic thats ~dcorr^2 and will be also returned.

Explanation: Distance correlation involves computing the normalized covariance of two centered euclidean distance matrices. Each distance matrix is the euclidean distance between rows (if x or y are 2d) or scalars (if x or y are 1d). Each matrix is centered prior to computing the covariance either using double-centering or u-centering, which corrects for bias as the number of dimensions increases. U-centering is almost always preferred in all cases. It also permits inference of the normalized covariance between each distance matrix using a one-tailed directional t-test. (Szekely & Rizzo, 2013). While distance correlation is normally bounded between 0 and 1, u-centering can produce negative estimates, which are never significant.

Validated against the dcor and dcor.ttest functions in the ‘energy’ R package and the dcor.distance_correlation, dcor.udistance_correlation_sqr, and dcor.independence.distance_correlation_t_test functions in the dcor Python package.

Parameters:

Returns:

Examples:

>>> import numpy as np
>>> x = np.random.randn(20, 3)
>>> y = x + np.random.randn(20, 3) * 0.1  # Strongly correlated
>>> result = distance_correlation(x, y, bias_corrected=True)
>>> 'dcorr' in result
True
>>> 0 <= result['dcorr'] <= 1
True

double_center¶

double_center(mat: np.ndarray) -> np.ndarray

Double center a 2d array.

Double-centering subtracts row means, column means, and adds the grand mean. This centers both rows and columns around zero.

Parameters:

Returns:

Examples:

>>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
>>> result = double_center(mat)
>>> np.allclose(result.mean(axis=0), 0)
True
>>> np.allclose(result.mean(axis=1), 0)
True

isc_group_permutation_test¶

isc_group_permutation_test(group1: np.ndarray, group2: np.ndarray, *, n_permute: int = 5000, summary: Literal['median', 'mean'] = 'median', method: Literal['permute', 'bootstrap'] = 'permute', summary_statistic: Literal['leave-one-out', 'pairwise'] = 'pairwise', ci_percentile: float = 95, tail: int | str = 2, device: Literal['cpu', 'gpu'] | None = 'cpu', n_jobs: int = -1, random_state: int | None = None, return_null: bool = False, progress_bar: bool = False, exclude_self_corr: bool = True, metric: str = 'correlation') -> dict[str, Any]

Compute ISC difference between groups with permutation testing.

Supports both subject-wise permutation and bootstrap methods with efficient CPU-parallel and optional GPU acceleration. Follows the statistical methods from Chen et al. (2016) for correct group comparison inference.

Parameters:

Returns:

Examples:

>>> # Single-feature ISC group comparison
>>> group1 = np.random.randn(100, 10)  # 10 subjects
>>> group2 = np.random.randn(100, 10)
>>> result = isc_group_permutation_test(group1, group2, n_permute=1000)
>>> print(f"ISC difference: {result['isc_group_difference']:.3f}, p: {result['p']:.3f}")

>>> # Voxel-wise ISC group comparison with GPU acceleration
>>> group1_voxels = np.random.randn(100, 10, 5000)  # 5K voxels
>>> group2_voxels = np.random.randn(100, 10, 5000)
>>> result = isc_group_permutation_test(
...     group1_voxels,
...     group2_voxels,
...     summary_statistic='leave-one-out',
...     device='gpu',  # GPU for LOO computation
...     n_permute=5000
... )
>>> print(f"Significant voxels: {(result['p'] < 0.05).sum()}")

isc_permutation_test¶

isc_permutation_test(data: np.ndarray, *, n_permute: int = 5000, summary: Literal['median', 'mean'] = 'median', summary_statistic: Literal['leave-one-out', 'pairwise'] = 'pairwise', method: Literal['bootstrap', 'circle_shift', 'phase_randomize'] = 'bootstrap', ci_percentile: float = 95, tail: int | str = 2, return_null: bool = False, progress_bar: bool = False, exclude_self_corr: bool = True, metric: str = 'correlation', device: Literal['cpu', 'gpu'] | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None) -> dict[str, Any]

Compute intersubject correlation with permutation testing.

Supports both leave-one-out and pairwise ISC computation modes with GPU acceleration for large voxel-wise problems and CPU-parallel bootstrap resampling.

Parameters:

Returns:

Examples:

>>> # Single-feature ISC
>>> data = np.random.randn(100, 10)  # 100 timepoints, 10 subjects
>>> result = isc_permutation_test(data, n_permute=1000)
>>> print(f"ISC: {result['isc']:.3f}, p: {result['p']:.3f}")

>>> # Voxel-wise ISC with GPU acceleration
>>> data_voxels = np.random.randn(100, 50, 5000)  # 5K voxels
>>> result = isc_permutation_test(
...     data_voxels,
...     summary_statistic='leave-one-out',
...     device='gpu',  # GPU for LOO computation
...     n_permute=5000
... )
>>> print(f"Significant voxels: {(result['p'] < 0.05).sum()}")

>>> # Compare LOO vs pairwise
>>> result_loo = isc_permutation_test(data, summary_statistic='leave-one-out')
>>> result_pair = isc_permutation_test(data, summary_statistic='pairwise')
>>> print(f"LOO: {result_loo['isc']:.3f}, Pairwise: {result_pair['isc']:.3f}")

matrix_permutation_test¶

matrix_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, metric: str = 'pearson', how: str = 'upper', include_diag: bool = False, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, random_state: int | None = None, progress_bar: bool = False) -> dict

Matrix permutation test (Mantel test) for correlating two square matrices.

Tests whether the correlation between elements of two matrices is significant by permuting rows and columns of one matrix symmetrically while keeping the other fixed.

Statistical Method: For each permutation, create random permutation perm, then apply: matrix1[perm][:, perm]. This preserves matrix structure while destroying correlation. Count how often permuted correlation is as extreme as observed.

Assumptions:

• Matrices are square and same size

• Under H₀, row/column ordering is exchangeable

• Symmetric permutation preserves matrix properties (e.g., symmetry)

Parameters:

Returns:

Examples:

>>> import numpy as np
>>> from nltools.algorithms.inference import matrix_permutation_test
>>>
>>> # Create two correlated similarity matrices
>>> np.random.seed(42)
>>> n = 50
>>> true_pattern = np.random.randn(n)
>>> data1 = np.corrcoef(true_pattern + np.random.randn(n) * 0.1)
>>> data2 = np.corrcoef(true_pattern + np.random.randn(n) * 0.1)
>>>
>>> # Test if matrices are correlated
>>> result = matrix_permutation_test(data1, data2, n_permute=1000)
>>> print(f"Correlation: {result['correlation']:.3f}, p = {result['p']:.4f}")

one_sample_permutation_test¶

one_sample_permutation_test(data: np.ndarray, *, n_permute: int = 5000, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

One-sample permutation test using sign-flipping.

Tests whether the mean of data is significantly different from zero by randomly flipping the sign of each observation. This is the permutation test equivalent of a one-sample t-test.

Assumption: Symmetric error distribution around zero. For highly skewed distributions, consider alternative methods (e.g., bootstrap resampling).

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> data = np.random.randn(30)
>>> result = one_sample_permutation_test(data, n_permute=5000)
>>> result['p']
0.23

>>> # Voxel-wise test with GPU
>>> data = np.random.randn(30, 10000)  # 30 subjects, 10K voxels
>>> result = one_sample_permutation_test(data, n_permute=5000, device='gpu')
>>> result['mean'].shape
(10000,)
>>> result['p'].shape
(10000,)

>>> # Single-threaded (for debugging)
>>> result = one_sample_permutation_test(data, n_permute=5000, device=None)

phase_randomize¶

phase_randomize(data: np.ndarray, *, device: str | None = 'cpu', random_state: int | np.random.RandomState | None = None) -> np.ndarray

FFT-based phase randomization for time-series data.

Preserves the power spectrum (autocorrelation) but destroys nonlinear temporal structure by randomizing Fourier phases. Used to test whether data was generated by a linear Gaussian process or contains nonlinear dynamics.

Parameters:

Returns:

Examples:

>>> x = np.sin(np.linspace(0, 10*np.pi, 100))  # Sine wave
>>> x_rand = phase_randomize(x, random_state=42)
>>> # Power spectrum preserved:
>>> np.allclose(np.abs(np.fft.rfft(x))**2, np.abs(np.fft.rfft(x_rand))**2)
True

>>> # GPU acceleration for large datasets:
>>> x_large = np.random.randn(10000)
>>> x_rand_gpu = phase_randomize(x_large, device='gpu', random_state=42)

timeseries_correlation_permutation_test¶

timeseries_correlation_permutation_test(data1: np.ndarray, data2: np.ndarray, *, method: Literal['circle_shift', 'phase_randomize'] = 'circle_shift', n_permute: int = 5000, metric: Literal['pearson', 'spearman', 'kendall'] = 'pearson', tail: int | str = 2, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, return_null: bool = False, random_state: int | np.random.RandomState | None = None, progress_bar: bool = False) -> dict

Time-series correlation permutation test.

Unlike standard permutation tests that shuffle data independently, this test uses time-series-aware permutation methods that preserve temporal structure (circle_shift) or power spectrum (phase_randomize).

Use this test when data contains temporal autocorrelation. Standard permutation tests inflate Type I error for autocorrelated data.

Parameters:

Returns:

Examples:

>>> x = np.sin(np.linspace(0, 10*np.pi, 100))
>>> y = np.cos(np.linspace(0, 10*np.pi, 100))
>>> result = timeseries_correlation_permutation_test(
...     x, y, method='circle_shift', n_permute=1000, random_state=42
... )
>>> result['correlation']  # Strong negative correlation
-0.999...
>>> result['p'] < 0.05  # Significant
True

>>> # GPU acceleration
>>> result = timeseries_correlation_permutation_test(
...     x, y, method='phase_randomize', device='gpu', n_permute=5000
... )

two_sample_permutation_test¶

two_sample_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

Two-sample permutation test using group label shuffling.

Tests whether two independent groups have different means by randomly permuting group labels. This is the permutation test equivalent of an independent samples t-test.

Assumption: Exchangeability under the null hypothesis (group assignments are arbitrary). Valid for independent samples from similar distributions.

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> data1 = np.random.randn(20)  # Group 1: 20 subjects
>>> data2 = np.random.randn(25)  # Group 2: 25 subjects
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000)
>>> result['p']
0.45

>>> # Voxel-wise test with GPU
>>> data1 = np.random.randn(20, 10000)  # 20 subjects, 10K voxels
>>> data2 = np.random.randn(25, 10000)  # 25 subjects, 10K voxels
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000, device='gpu')
>>> result['mean_diff'].shape
(10000,)
>>> result['p'].shape
(10000,)

>>> # Single-threaded (for debugging)
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000, device=None)

u_center¶

u_center(mat: np.ndarray) -> np.ndarray

U-center a 2d array. U-centering is a bias-corrected form of double-centering.

U-centering corrects for bias that occurs with double-centering as the number of dimensions increases. The diagonal is explicitly set to zero.

Parameters:

Returns:

Examples:

>>> mat = np.random.randn(5, 5)
>>> result = u_center(mat)
>>> np.allclose(np.diag(result), 0)
True

Modules¶

bootstrap¶

Bootstrap inference utilities with CPU/GPU support.

Attributes:

Classes:

Classes¶

OnlineBootstrapStats¶

OnlineBootstrapStats(shape: tuple[int, ...], save_samples: bool = False, percentiles: tuple[float, float] = (2.5, 97.5))

Memory-efficient online statistics aggregator for bootstrap samples.

Uses Welford’s algorithm for numerically stable online computation of mean and variance. Optionally stores all samples for exact percentile CIs.

Parameters:

Attributes:

####### Attributes##

Methods:

Examples:

>>> stats = OnlineBootstrapStats(shape=(100,), save_samples=False)
>>> for i in range(1000):
...     sample = np.random.randn(100)
...     stats.update(sample)
>>> results = stats.get_results()
>>> print(results.keys())
dict_keys(['mean', 'std', 'Z', 'p', 'ci_lower', 'ci_upper'])

M2¶

M2 = np.zeros(shape, dtype=(np.float64))

######## mean

mean = np.zeros(shape, dtype=(np.float64))

######## n

n = 0

######## percentiles

percentiles = percentiles

######## samples

samples = [] if save_samples else None

######## save_samples

save_samples = save_samples

######## shape

shape = shape

####### Functions##

get_results¶

get_results(tail: int | str = 2) -> dict[str, np.ndarray]

Compute final bootstrap statistics.

Parameters:

Parameters:

Returns:

Examples:

stats = OnlineBootstrapStats(shape=(100,), save_samples=False)
for _ in range(1000):
    stats.update(np.random.randn(100))
results = stats.get_results()
# results.keys() -> mean, std, Z, p, ci_lower, ci_upper

######## update

update(sample: np.ndarray) -> None

Update statistics with a new bootstrap sample.

Uses Welford’s algorithm for numerical stability.

Methods¶

correlation¶

Correlation permutation test implementations.

This module provides CPU-parallel and GPU-batched implementations of correlation permutation tests for assessing statistical significance of correlations.

Methods:

Classes¶

Methods¶

correlation_permutation_test¶

correlation_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, metric: str = 'pearson', tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

Correlation permutation test.

Tests whether the correlation between data1 and data2 is significantly different from zero by randomly permuting data1 and computing correlations.

Assumption: Observations are independent (i.i.d.). For autocorrelated time series, use timeseries_correlation_permutation_test with circle_shift or phase_randomize methods instead.

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> x = np.random.randn(100)
>>> y = x + np.random.randn(100) * 0.5  # Correlated
>>> result = correlation_permutation_test(x, y, n_permute=5000)
>>> result['correlation']
0.85
>>> result['p']
0.001

>>> # Multi-feature (2D arrays)
>>> data1 = np.random.randn(100, 10)  # 100 samples, 10 features
>>> data2 = data1 + np.random.randn(100, 10) * 0.3  # Correlated
>>> result = correlation_permutation_test(data1, data2, n_permute=5000)
>>> result['correlation'].shape
(10,)
>>> result['p'].shape
(10,)

>>> # GPU acceleration
>>> result = correlation_permutation_test(data1, data2, n_permute=5000, device='gpu')

intersubject¶

Intersubject correlation, functional connectivity, and phase synchrony.

Methods:

Methods¶

isc¶

isc(data, *, n_samples = 5000, summary = 'median', method = 'bootstrap', ci_percentile = 95, exclude_self_corr = True, tail = 2, metric = 'correlation', return_null = False, n_jobs = -1, random_state = None, progress_bar = False)

Compute pairwise intersubject correlation from observations by subjects array.

This function computes pairwise intersubject correlations (ISC) using the median as recommended by Chen et al., 2016). However, if the mean is preferred, we compute the mean correlation after performing the fisher r-to-z transformation and then convert back to correlations to minimize artificially inflating the correlation values.

There are currently three different methods to compute p-values. These include the classic methods for computing permuted time-series by either circle-shifting the data or phase-randomizing the data (see Lancaster et al., 2018). These methods create random surrogate data while preserving the temporal autocorrelation inherent to the signal. By default, we use the subject-wise bootstrap method from Chen et al., 2016. Instead of recomputing the pairwise ISC using circle_shift or phase_randomization methods, this approach uses the computationally more efficient method of bootstrapping the subjects and computing a new pairwise similarity matrix with randomly selected subjects with replacement. If the same subject is selected multiple times, we set the perfect correlation to a nan with (exclude_self_corr=True). We compute the p-values using the percentile method using the same method in Brainiak.

Chen, G., Shin, Y. W., Taylor, P. A., Glen, D. R., Reynolds, R. C., Israel, R. B., & Cox, R. W. (2016). Untangling the relatedness among correlations, part I: nonparametric approaches to inter-subject correlation analysis at the group level. NeuroImage, 142, 248-259.

Hall, P., & Wilson, S. R. (1991). Two guidelines for bootstrap hypothesis testing. Biometrics, 757-762.

Lancaster, G., Iatsenko, D., Pidde, A., Ticcinelli, V., & Stefanovska, A. (2018). Surrogate data for hypothesis testing of physical systems. Physics Reports, 748, 1-60.

This function is a wrapper around isc_permutation_test from the inference module, which provides optimized implementations with CPU-parallel and GPU acceleration support.

Parameters:

Returns:

isc_group¶

isc_group(group1, group2, *, n_samples = 5000, summary = 'median', method = 'permute', ci_percentile = 95, exclude_self_corr = True, return_null = False, tail = 2, metric = 'correlation', n_jobs = -1, random_state = None, progress_bar = False)

Compute difference in intersubject correlation between groups.

This function computes pairwise intersubject correlations (ISC) using the median as recommended by Chen et al., 2016). However, if the mean is preferred, we compute the mean correlation after performing the fisher r-to-z transformation and then convert back to correlations to minimize artificially inflating the correlation values.

There are currently two different methods to compute p-values. By default, we use the subject-wise permutation method recommended Chen et al., 2016. This method combines the two groups and computes pairwise similarity both within and between the groups. Then the group labels are permuted and the mean difference between the two groups are recomputed to generate a null distribution. The second method uses subject-wise bootstrapping, where a new pairwise similarity matrix with randomly selected subjects with replacement is created separately for each group and the ISC difference between these groups is used to generate a null distribution. If the same subject is selected multiple times, we set the perfect correlation to a nan with (exclude_self_corr=True). We compute the p-values using the percentile method (Hall & Wilson, 1991).

Chen, G., Shin, Y. W., Taylor, P. A., Glen, D. R., Reynolds, R. C., Israel, R. B., & Cox, R. W. (2016). Untangling the relatedness among correlations, part I: nonparametric approaches to inter-subject correlation analysis at the group level. NeuroImage, 142, 248-259.

Hall, P., & Wilson, S. R. (1991). Two guidelines for bootstrap hypothesis testing. Biometrics, 757-762.

This function is a thin wrapper around isc_group_permutation_test from the inference module (which provides optimized CPU parallelization and optional GPU acceleration), pinning the classic pairwise behavior and the n_samples vocabulary.

Parameters:

Returns:

isfc¶

isfc(data, method = 'average', n_jobs = -1)

Compute intersubject functional connectivity (ISFC) from a list of observation x feature matrices.

This function uses the leave one out approach to compute ISFC (Simony et al., 2016). For each subject, compute the cross-correlation between each voxel/roi with the average of the rest of the subjects data. In other words, compute the mean voxel/ROI response for all participants except the target subject. Then compute the correlation between each ROI within the target subject with the mean ROI response in the group average.

Simony, E., Honey, C. J., Chen, J., Lositsky, O., Yeshurun, Y., Wiesel, A., & Hasson, U. (2016). Dynamic reconfiguration of the default mode network during narrative comprehension. Nature communications, 7, 12141.

This function now uses the optimized implementation from the inference module, which provides efficient cross-correlation computation between matrix columns. CPU parallelization is available via joblib when n_jobs > 1 or n_jobs=-1. Each subject’s ISFC computation is independent and can be parallelized efficiently.

Parameters:

Returns:

isps¶

isps(data, *, sampling_freq = 0.5, low_cut = 0.04, high_cut = 0.07, order = 5, pairwise = False)

Compute dynamic intersubject phase synchrony (ISPS) from an observations-by-subjects array.

This function computes the instantaneous intersubject phase synchrony for a single voxel/roi timeseries. Requires multiple subjects. This method is largely based on that described by Glerean et al., 2012 and performs a hilbert transform on narrow bandpass filtered timeseries (butterworth) data to get the instantaneous phase angle. The function returns a dictionary containing the average phase angle, the average vector length, and parametric p-values computed using the rayleigh test using circular statistics (Fisher, 1993). If pairwise=True, then it will compute these on the pairwise phase angle differences, if pairwise=False, it will compute these on the actual phase angles. This is called inter-site phase coupling or inter-trial phase coupling respectively in the EEG literatures.

This function requires narrow band filtering your data. As a default we use the recommendations by (Glerean et al., 2012) of .04-.07Hz. This is similar to the “slow-4” band (0.025–0.067 Hz) described by (Zuo et al., 2010; Penttonen & Buzsáki, 2003), but excludes the .03 band, which has been demonstrated to contain aliased respiration signals (Birn, 2006).

Birn RM, Smith MA, Bandettini PA, Diamond JB. 2006. Separating respiratory-variation-related fluctuations from neuronal-activity- related fluctuations in fMRI. Neuroimage 31:1536–1548.

Buzsáki, G., & Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science, 304(5679), 1926-1929.

Fisher, N. I. (1995). Statistical analysis of circular data. cambridge university press.

Glerean, E., Salmi, J., Lahnakoski, J. M., Jääskeläinen, I. P., & Sams, M. (2012). Functional magnetic resonance imaging phase synchronization as a measure of dynamic functional connectivity. Brain connectivity, 2(2), 91-101.

Parameters:

Returns:

isc¶

Intersubject Correlation (ISC) with GPU-Accelerated Permutation Testing.

This module provides both leave-one-out (LOO) and pairwise ISC computation with efficient CPU-parallel and GPU-batched implementations. Follows the statistical methods from Chen et al. (2016) for correct bootstrap resampling of correlation matrices.

Methods:

Methods¶

isc_group_permutation_test¶

isc_group_permutation_test(group1: np.ndarray, group2: np.ndarray, *, n_permute: int = 5000, summary: Literal['median', 'mean'] = 'median', method: Literal['permute', 'bootstrap'] = 'permute', summary_statistic: Literal['leave-one-out', 'pairwise'] = 'pairwise', ci_percentile: float = 95, tail: int | str = 2, device: Literal['cpu', 'gpu'] | None = 'cpu', n_jobs: int = -1, random_state: int | None = None, return_null: bool = False, progress_bar: bool = False, exclude_self_corr: bool = True, metric: str = 'correlation') -> dict[str, Any]

Compute ISC difference between groups with permutation testing.

Supports both subject-wise permutation and bootstrap methods with efficient CPU-parallel and optional GPU acceleration. Follows the statistical methods from Chen et al. (2016) for correct group comparison inference.

Parameters:

Returns:

Examples:

>>> # Single-feature ISC group comparison
>>> group1 = np.random.randn(100, 10)  # 10 subjects
>>> group2 = np.random.randn(100, 10)
>>> result = isc_group_permutation_test(group1, group2, n_permute=1000)
>>> print(f"ISC difference: {result['isc_group_difference']:.3f}, p: {result['p']:.3f}")

>>> # Voxel-wise ISC group comparison with GPU acceleration
>>> group1_voxels = np.random.randn(100, 10, 5000)  # 5K voxels
>>> group2_voxels = np.random.randn(100, 10, 5000)
>>> result = isc_group_permutation_test(
...     group1_voxels,
...     group2_voxels,
...     summary_statistic='leave-one-out',
...     device='gpu',  # GPU for LOO computation
...     n_permute=5000
... )
>>> print(f"Significant voxels: {(result['p'] < 0.05).sum()}")

isc_permutation_test¶

isc_permutation_test(data: np.ndarray, *, n_permute: int = 5000, summary: Literal['median', 'mean'] = 'median', summary_statistic: Literal['leave-one-out', 'pairwise'] = 'pairwise', method: Literal['bootstrap', 'circle_shift', 'phase_randomize'] = 'bootstrap', ci_percentile: float = 95, tail: int | str = 2, return_null: bool = False, progress_bar: bool = False, exclude_self_corr: bool = True, metric: str = 'correlation', device: Literal['cpu', 'gpu'] | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None) -> dict[str, Any]

Compute intersubject correlation with permutation testing.

Supports both leave-one-out and pairwise ISC computation modes with GPU acceleration for large voxel-wise problems and CPU-parallel bootstrap resampling.

Parameters:

Returns:

Examples:

>>> # Single-feature ISC
>>> data = np.random.randn(100, 10)  # 100 timepoints, 10 subjects
>>> result = isc_permutation_test(data, n_permute=1000)
>>> print(f"ISC: {result['isc']:.3f}, p: {result['p']:.3f}")

>>> # Voxel-wise ISC with GPU acceleration
>>> data_voxels = np.random.randn(100, 50, 5000)  # 5K voxels
>>> result = isc_permutation_test(
...     data_voxels,
...     summary_statistic='leave-one-out',
...     device='gpu',  # GPU for LOO computation
...     n_permute=5000
... )
>>> print(f"Significant voxels: {(result['p'] < 0.05).sum()}")

>>> # Compare LOO vs pairwise
>>> result_loo = isc_permutation_test(data, summary_statistic='leave-one-out')
>>> result_pair = isc_permutation_test(data, summary_statistic='pairwise')
>>> print(f"LOO: {result_loo['isc']:.3f}, Pairwise: {result_pair['isc']:.3f}")

matrix¶

Matrix permutation test implementations (Mantel test).

This module provides CPU-parallel implementations of matrix permutation tests for testing correlation between two square matrices, as well as matrix utility functions for distance correlation and matrix centering operations.

Attributes:

Methods:

Methods¶

distance_correlation¶

distance_correlation(x: np.ndarray, y: np.ndarray, bias_corrected: bool = True, ttest: bool = False) -> dict

Compute the distance correlation between 2 arrays to test for multivariate dependence (linear or non-linear).

Arrays must match on their first dimension. It’s almost always preferable to compute the bias_corrected version which can also optionally perform a ttest. This ttest operates on a statistic thats ~dcorr^2 and will be also returned.

Explanation: Distance correlation involves computing the normalized covariance of two centered euclidean distance matrices. Each distance matrix is the euclidean distance between rows (if x or y are 2d) or scalars (if x or y are 1d). Each matrix is centered prior to computing the covariance either using double-centering or u-centering, which corrects for bias as the number of dimensions increases. U-centering is almost always preferred in all cases. It also permits inference of the normalized covariance between each distance matrix using a one-tailed directional t-test. (Szekely & Rizzo, 2013). While distance correlation is normally bounded between 0 and 1, u-centering can produce negative estimates, which are never significant.

Validated against the dcor and dcor.ttest functions in the ‘energy’ R package and the dcor.distance_correlation, dcor.udistance_correlation_sqr, and dcor.independence.distance_correlation_t_test functions in the dcor Python package.

Parameters:

Returns:

Examples:

>>> import numpy as np
>>> x = np.random.randn(20, 3)
>>> y = x + np.random.randn(20, 3) * 0.1  # Strongly correlated
>>> result = distance_correlation(x, y, bias_corrected=True)
>>> 'dcorr' in result
True
>>> 0 <= result['dcorr'] <= 1
True

double_center¶

double_center(mat: np.ndarray) -> np.ndarray

Double center a 2d array.

Double-centering subtracts row means, column means, and adds the grand mean. This centers both rows and columns around zero.

Parameters:

Returns:

Examples:

>>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
>>> result = double_center(mat)
>>> np.allclose(result.mean(axis=0), 0)
True
>>> np.allclose(result.mean(axis=1), 0)
True

matrix_permutation_test¶

matrix_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, metric: str = 'pearson', how: str = 'upper', include_diag: bool = False, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, random_state: int | None = None, progress_bar: bool = False) -> dict

Matrix permutation test (Mantel test) for correlating two square matrices.

Tests whether the correlation between elements of two matrices is significant by permuting rows and columns of one matrix symmetrically while keeping the other fixed.

Statistical Method: For each permutation, create random permutation perm, then apply: matrix1[perm][:, perm]. This preserves matrix structure while destroying correlation. Count how often permuted correlation is as extreme as observed.

Assumptions:

• Matrices are square and same size

• Under H₀, row/column ordering is exchangeable

• Symmetric permutation preserves matrix properties (e.g., symmetry)

Parameters:

Returns:

Examples:

>>> import numpy as np
>>> from nltools.algorithms.inference import matrix_permutation_test
>>>
>>> # Create two correlated similarity matrices
>>> np.random.seed(42)
>>> n = 50
>>> true_pattern = np.random.randn(n)
>>> data1 = np.corrcoef(true_pattern + np.random.randn(n) * 0.1)
>>> data2 = np.corrcoef(true_pattern + np.random.randn(n) * 0.1)
>>>
>>> # Test if matrices are correlated
>>> result = matrix_permutation_test(data1, data2, n_permute=1000)
>>> print(f"Correlation: {result['correlation']:.3f}, p = {result['p']:.4f}")

u_center¶

u_center(mat: np.ndarray) -> np.ndarray

U-center a 2d array. U-centering is a bias-corrected form of double-centering.

U-centering corrects for bias that occurs with double-centering as the number of dimensions increases. The diagonal is explicitly set to zero.

Parameters:

Returns:

Examples:

>>> mat = np.random.randn(5, 5)
>>> result = u_center(mat)
>>> np.allclose(np.diag(result), 0)
True

one_sample¶

One-sample permutation test implementations.

This module provides CPU-parallel and GPU-batched implementations of the one-sample permutation test (sign-flipping test).

Methods:

Classes¶

Methods¶

one_sample_permutation_test¶

one_sample_permutation_test(data: np.ndarray, *, n_permute: int = 5000, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

One-sample permutation test using sign-flipping.

Tests whether the mean of data is significantly different from zero by randomly flipping the sign of each observation. This is the permutation test equivalent of a one-sample t-test.

Assumption: Symmetric error distribution around zero. For highly skewed distributions, consider alternative methods (e.g., bootstrap resampling).

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> data = np.random.randn(30)
>>> result = one_sample_permutation_test(data, n_permute=5000)
>>> result['p']
0.23

>>> # Voxel-wise test with GPU
>>> data = np.random.randn(30, 10000)  # 30 subjects, 10K voxels
>>> result = one_sample_permutation_test(data, n_permute=5000, device='gpu')
>>> result['mean'].shape
(10000,)
>>> result['p'].shape
(10000,)

>>> # Single-threaded (for debugging)
>>> result = one_sample_permutation_test(data, n_permute=5000, device=None)

timeseries¶

Time-series permutation test implementations.

This module provides GPU-accelerated implementations of time-series permutation tests that preserve temporal structure:

• circle_shift: Circular shift permutation (preserves autocorrelation)

• phase_randomize: FFT-based phase randomization (preserves power spectrum)

• timeseries_correlation_permutation_test: Correlation test with timeseries methods

Methods:

Classes¶

Methods¶

circle_shift¶

circle_shift(data: np.ndarray, shift_amount: int | np.ndarray | None = None, random_state: int | np.random.RandomState | None = None) -> np.ndarray

Circular shift for time-series data.

Performs a circular shift that preserves autocorrelation structure. Useful for permutation tests on autocorrelated time series (e.g., fMRI). For 1D data, shifts by a single amount. For 2D data, shifts each feature (column) independently.

Parameters:

Returns:

Examples:

>>> x = np.array([1, 2, 3, 4, 5])
>>> circle_shift(x, shift_amount=2)
array([4, 5, 1, 2, 3])

>>> X = np.array([[1, 10], [2, 20], [3, 30], [4, 40]])
>>> circle_shift(X, shift_amount=np.array([1, 2]))
array([[ 4, 30],
       [ 1, 40],
       [ 2, 10],
       [ 3, 20]])

phase_randomize¶

phase_randomize(data: np.ndarray, *, device: str | None = 'cpu', random_state: int | np.random.RandomState | None = None) -> np.ndarray

FFT-based phase randomization for time-series data.

Preserves the power spectrum (autocorrelation) but destroys nonlinear temporal structure by randomizing Fourier phases. Used to test whether data was generated by a linear Gaussian process or contains nonlinear dynamics.

Parameters:

Returns:

Examples:

>>> x = np.sin(np.linspace(0, 10*np.pi, 100))  # Sine wave
>>> x_rand = phase_randomize(x, random_state=42)
>>> # Power spectrum preserved:
>>> np.allclose(np.abs(np.fft.rfft(x))**2, np.abs(np.fft.rfft(x_rand))**2)
True

>>> # GPU acceleration for large datasets:
>>> x_large = np.random.randn(10000)
>>> x_rand_gpu = phase_randomize(x_large, device='gpu', random_state=42)

timeseries_correlation_permutation_test¶

timeseries_correlation_permutation_test(data1: np.ndarray, data2: np.ndarray, *, method: Literal['circle_shift', 'phase_randomize'] = 'circle_shift', n_permute: int = 5000, metric: Literal['pearson', 'spearman', 'kendall'] = 'pearson', tail: int | str = 2, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, return_null: bool = False, random_state: int | np.random.RandomState | None = None, progress_bar: bool = False) -> dict

Time-series correlation permutation test.

Unlike standard permutation tests that shuffle data independently, this test uses time-series-aware permutation methods that preserve temporal structure (circle_shift) or power spectrum (phase_randomize).

Use this test when data contains temporal autocorrelation. Standard permutation tests inflate Type I error for autocorrelated data.

Parameters:

Returns:

Examples:

>>> x = np.sin(np.linspace(0, 10*np.pi, 100))
>>> y = np.cos(np.linspace(0, 10*np.pi, 100))
>>> result = timeseries_correlation_permutation_test(
...     x, y, method='circle_shift', n_permute=1000, random_state=42
... )
>>> result['correlation']  # Strong negative correlation
-0.999...
>>> result['p'] < 0.05  # Significant
True

>>> # GPU acceleration
>>> result = timeseries_correlation_permutation_test(
...     x, y, method='phase_randomize', device='gpu', n_permute=5000
... )

two_sample¶

Two-sample permutation test implementations.

This module provides CPU-parallel and GPU-batched implementations of the two-sample permutation test (group permutation test).

Methods:

Classes¶

Methods¶

two_sample_permutation_test¶

two_sample_permutation_test(data1: np.ndarray, data2: np.ndarray, *, n_permute: int = 5000, tail: int | str = 2, return_null: bool = False, device: str | None = 'cpu', n_jobs: int = -1, max_gpu_memory_gb: float | None = None, random_state: int | None = None, progress_bar: bool = False) -> dict

Two-sample permutation test using group label shuffling.

Tests whether two independent groups have different means by randomly permuting group labels. This is the permutation test equivalent of an independent samples t-test.

Assumption: Exchangeability under the null hypothesis (group assignments are arbitrary). Valid for independent samples from similar distributions.

Parameters:

Returns:

Examples:

>>> # Single feature (default CPU parallelization)
>>> data1 = np.random.randn(20)  # Group 1: 20 subjects
>>> data2 = np.random.randn(25)  # Group 2: 25 subjects
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000)
>>> result['p']
0.45

>>> # Voxel-wise test with GPU
>>> data1 = np.random.randn(20, 10000)  # 20 subjects, 10K voxels
>>> data2 = np.random.randn(25, 10000)  # 25 subjects, 10K voxels
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000, device='gpu')
>>> result['mean_diff'].shape
(10000,)
>>> result['p'].shape
(10000,)

>>> # Single-threaded (for debugging)
>>> result = two_sample_permutation_test(data1, data2, n_permute=5000, device=None)

utils¶

Utility functions for permutation testing.

This module contains shared helper functions used across different permutation test implementations.

Attributes:

Methods¶

validation¶

Shared validation utilities for algorithms module.

This module provides common validation functions to reduce code duplication and ensure consistent error handling across the algorithms module.

Methods:

Methods¶

validate_array_shape¶

validate_array_shape(array: np.ndarray, expected_ndim: int, name: str = 'array') -> None

Validate array dimensionality.

Parameters:

validate_array_shape_range¶

validate_array_shape_range(array: np.ndarray, min_ndim: int, max_ndim: int, name: str = 'array') -> None

Validate array dimensionality is within a range.

Parameters:

validate_bootstrap_data¶

validate_bootstrap_data(data: np.ndarray, method: str) -> None

Validate input data for bootstrapping.

Parameters:

validate_bootstrap_method¶

validate_bootstrap_method(method: str, simple_methods: list[str], fitted_methods: list[str]) -> None

Validate bootstrap method name.

Parameters:

validate_device_parameter¶

validate_device_parameter(device: str | None, *, allow_auto: bool = False) -> None

Validate device parameter.

Parameters:

validate_device_parameter_matrix¶

validate_device_parameter_matrix(device: str | None) -> None

Validate device parameter for matrix operations.

Parameters:

validate_how_parameter¶

validate_how_parameter(how: str) -> None

Validate ‘how’ parameter for matrix operations.

Parameters:

validate_metric_parameter¶

validate_metric_parameter(metric: str, allowed: list[str], name: str = 'metric') -> None

Validate metric parameter.

Parameters:

validate_percentiles¶

validate_percentiles(percentiles: tuple[float, float]) -> None

Validate percentile values for confidence intervals.

Parameters:

validate_same_shape¶

validate_same_shape(array1: np.ndarray, array2: np.ndarray, name1: str = 'array1', name2: str = 'array2') -> None

Validate two arrays have same shape.

Parameters:

validate_shape_compatibility¶

validate_shape_compatibility(X: np.ndarray, y: np.ndarray, X_name: str = 'X', y_name: str = 'y') -> None

Validate that X and y have compatible shapes for regression.

Parameters:

validate_square_matrix¶

validate_square_matrix(matrix: np.ndarray, name: str = 'matrix') -> None

Validate matrix is square.

Parameters:

validate_tail_parameter¶

validate_tail_parameter(tail: int | str) -> str

Validate the public tail vocabulary and normalize to the internal form.

The public vocabulary (v0.6.0) is deliberately two-valued — the direction of a one-tailed test is fixed by the test’s convention, never chosen from the data (a data-driven direction would silently halve every p-value):

Parameters:

Returns: