diff --git a/examples/README.txt b/examples/README.txt new file mode 100644 index 0000000..c4e7115 --- /dev/null +++ b/examples/README.txt @@ -0,0 +1,6 @@ +Examples Gallery +================ + +.. contents:: Contents + :local: + :depth: 1 diff --git a/examples/compute_auroc_pvalue.py b/examples/compute_auroc_pvalue.py new file mode 100644 index 0000000..9222706 --- /dev/null +++ b/examples/compute_auroc_pvalue.py @@ -0,0 +1,84 @@ +""" +=============================================================================== +Compute the p-value associated to AUROC +=============================================================================== + +Compute the p-value associated to the area under the curve (AUC) of the +receiver operating characteristic (ROC), for a binary (two-class) +classification problem. + +Warning: this example requires scikit-learn dependency (0.22 at least). +""" +# Authors: Quentin Barthélemy + +import numpy as np +from sklearn.datasets import make_classification +from sklearn.model_selection import train_test_split +from sklearn.linear_model import LogisticRegression +from sklearn.metrics import roc_curve, roc_auc_score +from pypermut.ml import permutation_metric, standard_error_auroc +from pypermut.misc import pvals_to_stars +from matplotlib import pyplot as plt + + +############################################################################### +# Artificial data and classifier +# ------------------------------ + +X, y = make_classification(n_classes=2, n_samples=100, n_features=100, + n_informative=10, n_redundant=10, + n_clusters_per_class=5, random_state=1) +Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.5, + random_state=2) +# Logistic regression: a very good old-fashioned classifier +clf = LogisticRegression().fit(Xtrain, ytrain) +ytest_probas = clf.predict_proba(Xtest)[:, 1] + + +############################################################################### +# Compute ROC curve and AUC +# ------------------------- + +# In this example, we focus on the AUROC metric. However, this analysis can be +# applied to any metric measuring the model performance. It can be any function +# from sklearn.metrics, like roc_auc_score, log_loss, etc. + +# ROC curve +fpr, tpr, _ = roc_curve(ytest, ytest_probas) + +fig, ax = plt.subplots() +ax.set(title="ROC curve", xlabel='FPR', ylabel='TPR') +ax.plot(fpr, tpr, 'C1', label='Logistic regression') +ax.plot([0, 1], [0, 1], 'b--', label='Chance level') +plt.xlim([0.0, 1.0]) +plt.ylim([0.0, 1.05]) +plt.legend(loc="lower right") +plt.show() + +# AUC +auroc = roc_auc_score(ytest, ytest_probas) +print('AUROC = {:.3f}\n'.format(auroc)) + + +############################################################################### +# Compute standard error and p-value of AUC +# ----------------------------------------- + +n_perm = 5000 +np.random.seed(17) + +# standard error of AUC +auroc_se = standard_error_auroc(ytest, ytest_probas, auroc) + +# p-value of AUC, by permutations +auroc_, pval = permutation_metric(ytest, ytest_probas, roc_auc_score, + side='right', n=n_perm) +print('AUROC = {:.3f} +/- {:.3f}, p={:.2e} ({})' + .format(auroc_, auroc_se, pval, pvals_to_stars(pval))) + + +############################################################################### +# Reference +# --------- +# [1] Hanley & McNeil, "The meaning and use of the area under a receiver +# operating characteristic (ROC) curve", Radiology, 1982. diff --git a/examples/compute_exact_ttest_ind.py b/examples/compute_exact_ttest_ind.py new file mode 100644 index 0000000..2c8301b --- /dev/null +++ b/examples/compute_exact_ttest_ind.py @@ -0,0 +1,142 @@ +""" +=============================================================================== +Compute exact Student's t-test for independent samples +=============================================================================== + +Compute classical, permutated and exact right-sided Student's t-tests, for +independent / unpaired samples. +""" +# Authors: Quentin Barthélemy + +import numpy as np +from scipy.stats import t, ttest_ind, mannwhitneyu, shapiro, levene +from pypermut.stats import permutation_ttest_ind, permutation_mannwhitneyu +from pypermut.misc import pvals_to_stars +from matplotlib import pyplot as plt + + +############################################################################### +# Univariate data +# --------------- + +# Artificial small samples +s1 = np.array([7.5, 7.9, 8.3, 8.5, 8.7, 9.2, 9.6, 10.2]) +s2 = np.array([5.9, 6.5, 6.9, 7.4, 7.8, 7.9, 8.1, 8.5, 8.7, 9.1, 9.9]) + +# Plot samples +fig, ax = plt.subplots() +ax.set(title='Boxplot of samples', xlabel='Samples', ylabel='Values') +box = ax.boxplot(np.array([s1, s2], dtype=object), medianprops={'color':'k'}, + showmeans=True, meanline=True, meanprops={'color':'r'}) +ax.plot([1.1, 1.9], [m.get_ydata()[0] for m in box.get('means')], c='r') +ax.scatter(np.ones(len(s1)), s1) +ax.scatter(2*np.ones(len(s2)), s2) +plt.show() + + +############################################################################### +# Student's t-tests for independent samples +# ----------------------------------------- + +# The null hypothesis is that these two samples s1 and s2 are drawn from the +# same Gaussian distribution. + +# Classical test +t_class, p_class = ttest_ind(s1, s2, alternative='greater') +print('Classical t-test:\n t={:.2f}, p={:.3e} ({})' + .format(t_class, p_class, pvals_to_stars(p_class))) + +# Permutated test +np.random.seed(17) +n_perm = 10000 +t_perm, p_perm = permutation_ttest_ind(s1, s2, n=n_perm, side='one') +print('Permutation t-test:\n t={:.2f}, p={:.3e} ({})' + .format(t_perm[0], p_perm[0], pvals_to_stars(p_perm[0]))) + +# Exact test +t_exact, p_exact, t_dist = permutation_ttest_ind(s1, s2, n='all', side='one', + return_dist=True) +print('Exact t-test:\n t={:.2f}, p={:.3e} ({})\n' + .format(t_exact[0], p_exact[0], pvals_to_stars(p_exact[0]))) + +# Plot the difference between approximate and exact right-sided p-values +fig = plt.figure(figsize=(14, 5)) +ax1 = fig.add_subplot(121, title='Approximate p-value', xlabel='t') +ax1.axvline(x=t_exact, c='r', label='Real statistic t') +df = len(s1) + len(s2) - 2 +x = np.linspace(t.ppf(0.001, df), t.ppf(0.999, df), 100) +ax1.plot(x, t.pdf(x, df), label='Approximate t-distribution') +ax1.fill_between(x, t.pdf(x, df), where=(x>=t_exact), facecolor='r', + label='Approximate p-value') +plt.legend(loc='upper left') +ax2 = fig.add_subplot(122, title='Exact p-value', xlabel='t', sharey=ax1) +ax2.axvline(x=t_exact, c='r', label='Real statistic t') +y, x, _ = ax2.hist(t_dist, 50, density=True, alpha=0.4, + label='Exact t-distribution') +ax2.fill_between(x[1:], y, where=(x[1:]>=t_exact), facecolor='r', step='pre', + label='Exact p-value') +ax2.set_xlim(ax1.get_xlim()) +plt.legend(loc='upper left') +plt.show() + +# Classical t-test detects a trend in the difference between these samples, +# but is not able to reject the null hypothesis (with alpha = 0.05). +# Permutation and exact tests reject the null hypothesis, the later giving the +# exact p-value. Amazing! + + +############################################################################### +# Assumptions +# ----------- + +# But, wait! Is Student's t-test valid for such data? +# We have not checked the assumptions related to this parametric test. + +alpha = 0.05 + +# Check homoscedasticity (homogeneity of variances) assumption, with a Levene's +# test +stat, p = levene(s1, s2) +print('Levene test:\n Equal var={}, p={:.3f} ({})' + .format(p < alpha, p, pvals_to_stars(p))) + +# Ok, but we could use a Welch's t-test, which does not assume equal variance. + +# Check normality / Gaussianity assumption, with Shapiro-Wilk tests +stat, p = shapiro(s1) +print('Shapiro-Wilk test, for first sample:\n Normality={}, p={:.3f} ({})' + .format(p < alpha, p, pvals_to_stars(p))) +stat, p = shapiro(s2) +print('Shapiro-Wilk test, for second sample:\n Normality={}, p={:.3f} ({})\n' + .format(p < alpha, p, pvals_to_stars(p))) + +# Consequently, neither Student's t-test nor Welch's t-test can be applied on +# these non-Gaussian data. We must use the non-parametric version. + + +############################################################################### +# Mann-Whitney U tests for unpaired samples +# ----------------------------------------- + +# The null hypothesis is that these two samples s1 and s2 are drawn from the +# same distribution, without assumption on its shape. + +# Classical test +U_class, p_class = mannwhitneyu(s1, s2, alternative='less') +print('Classical Mann-Whitney test:\n U={:.1f}, p={:.3e} ({})' + .format(U_class, p_class, pvals_to_stars(p_class))) + +# Permutated test +U_perm, p_perm = permutation_mannwhitneyu(s1, s2, n=n_perm) +print('Permutation Mann-Whitney test:\n U={:.1f}, p={:.3e} ({})' + .format(U_perm[0], p_perm[0], pvals_to_stars(p_perm[0]))) + +# Exact test +U_exact, p_exact = permutation_mannwhitneyu(s1, s2, n='all') +print('Exact Mann-Whitney test:\n U={:.1f}, p={:.3e} ({})' + .format(U_exact[0], p_exact[0], pvals_to_stars(p_exact[0]))) + +# Finally, we must conclude that we can't reject the null hypothesis. + + +############################################################################### diff --git a/examples/correct_multiple_correlations.py b/examples/correct_multiple_correlations.py new file mode 100644 index 0000000..28de028 --- /dev/null +++ b/examples/correct_multiple_correlations.py @@ -0,0 +1,103 @@ +""" +=============================================================================== +Correct Pearson and Spearman correlations on Anscombe's data +=============================================================================== + +Compare classical and permutated tests on Anscombe's data sets [1], +for Pearson and Spearman correlations after correction for multiple tests. +""" +# Authors: Quentin Barthélemy + +import numpy as np +from scipy.stats import pearsonr, spearmanr +from pypermut.stats import permutation_corr +from pypermut.misc import print_results, print_pvals +from matplotlib import pyplot as plt + + +############################################################################### +# Multivariate data +# ----------------- + +# Anscombe's trio: only the first three data sets of the Anscombe's quartet, +# because the last set does not share the same values x. +x = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5]) +Y = np.array([ + [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84,4.82, 5.68], + [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74], + [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73] +]).T + +# Anscombe's anti-trio: to highlight the difference between classical and +# permutated tests, the number of variables is artificially increased. +Y = np.c_[Y, -Y] +vlabels = ['y{}'.format(v + 1) for v in range(Y.shape[1])] + +# Plot Anscombe's trio and anti-trio +fig, ax = plt.subplots() +ax.set(title="Anscombe's trio and anti-trio", xlabel='x', ylabel='Y') +[plt.scatter(x, y, label='y{}'.format(i + 1)) for i, y in enumerate(Y.T)] +ax.legend() +plt.show() + + +############################################################################### +# Classical statistics +# -------------------- + +# Pearson test +r_class = [pearsonr(y, x) for y in Y.T] +print('Pearson tests:') +print_results(r_class, vlabels, 'r') + +# Spearman test +rho_class = [spearmanr(y, x) for y in Y.T] +print('Spearman tests:') +print_results(rho_class, vlabels, 'rho') + +# Correction for multiple tests, using Bonferroni's method: multiply p-values +# by the number of tests. +n_tests = Y.shape[1] + +p_corrected = np.array(r_class)[:, 1] * n_tests +print('\nPearson tests, after Bonferroni correction:') +print_pvals(p_corrected, vlabels) + +p_corrected = np.array(rho_class)[:, 1] * n_tests +print('Spearman tests, after Bonferroni correction:') +print_pvals(p_corrected, vlabels) + +# Bonferroni correction is a crucial step to adjust p-values of each variable +# for multiple tests, avoiding type-I errors (ie, spurious correlations). +# However, it makes the hypothesis that tests are independent, that is not the +# case here. + + +############################################################################### +# Permutation statistics +# ---------------------- + +n_perm = 10000 +np.random.seed(17) + +# Permutated Pearson test, corrected by Rmax +r_perm, p_perm = permutation_corr(Y, x=x, n=n_perm, side='two') +print('\nPermutated Pearson tests, with {} permutations:'.format(n_perm)) +print_results(np.c_[r_perm, p_perm], vlabels, 'r') + +# Permutated Spearman test, corrected by Rmax +rho_perm, p_perm = permutation_corr(Y, x=x, n=n_perm, corr='spearman', + side='two') +print('Permutated Spearman tests, with {} permutations:'.format(n_perm)) +print_results(np.c_[rho_perm, p_perm], vlabels, 'rho') + +# The Rmax method of the permutation correlation tests is used for adjusting +# the p-values in a way that controls the family-wise error rate. +# It is more powerful than Bonferroni correction when different variables are +# correlated: it avoids type-I and type-II errors. + + +############################################################################### +# Reference +# --------- +# [1] Anscombe, "Graphs in Statistical Analysis", Am Stat, 1973. diff --git a/examples/correct_multiple_ttest_ind.py b/examples/correct_multiple_ttest_ind.py new file mode 100644 index 0000000..93fd2c1 --- /dev/null +++ b/examples/correct_multiple_ttest_ind.py @@ -0,0 +1,111 @@ +""" +=============================================================================== +Correct multiple Student's t-tests for independent samples +=============================================================================== + +Correct multiple Student's t-tests for trivariate independent samples, +comparing three methods: Bonferroni, permutations, and Hotelling. +""" +# Authors: Quentin Barthélemy + +import numpy as np +from scipy.stats import ttest_ind, f +from pypermut.stats import permutation_ttest_ind +from pypermut.misc import print_results, print_pvals, pvals_to_stars +from matplotlib import pyplot as plt +from mpl_toolkits.mplot3d import Axes3D + + +############################################################################### +# Multivariate data +# ----------------- + +# Artificial trivariate Gaussian samples, with a shift in the third variable +n_meas, n_vars = 50, 3 +np.random.seed(42) +X = np.random.randn(n_meas, n_vars) +Y = np.random.randn(n_meas, n_vars) +Y[:, 2] -= 0.8 +vlabels = ['var{}'.format(v) for v in range(n_vars)] + +# Plot trivariate samples +fig = plt.figure(figsize=(12, 5)) +ax1 = fig.add_subplot(121, projection='3d', title='3D visualization', + xlabel=vlabels[0], ylabel=vlabels[1], zlabel=vlabels[2]) +ax1.scatter(X[:, 0], X[:, 1], X[:, 2], label='1st sample') +ax1.scatter(Y[:, 0], Y[:, 1], Y[:, 2], marker='^', label='2nd sample') +plt.legend(loc='center right') +ax2 = fig.add_subplot(122, projection='3d', title='2D projection', + xlabel=vlabels[0], ylabel=vlabels[1], zlabel=vlabels[2]) +for D, c in zip([X, Y], ['C0o', 'C1^']): + ax2.plot(D[:, 0], D[:, 1], c, zdir='z', zs=ax1.get_zlim()[0], ms=4, alpha=0.5) + ax2.plot(D[:, 0], D[:, 2], c, zdir='y', zs=ax1.get_ylim()[-1], ms=4, alpha=0.5) + ax2.plot(D[:, 1], D[:, 2], c, zdir='x', zs=ax1.get_xlim()[0], ms=4, alpha=0.5) +plt.show() + + +############################################################################### +# Student's t-tests for independent samples +# ----------------------------------------- + +# Classical tests +t_class = [ttest_ind(x, y, alternative='greater') for x, y in zip(X.T, Y.T)] +print('Classical t-tests:') +print_results(t_class, vlabels, 't') + +# Bonferroni correction, to control the family-wise error rate for multiple +# tests +p_corrected = np.array(t_class)[:, 1] * n_vars +print('Classical t-tests, after Bonferroni correction:') +print_pvals(p_corrected, vlabels) + +# Permutated tests +n_perm = 10000 +t_perm, p_perm = permutation_ttest_ind(X, Y, n=n_perm, side='one') +print('\nPermutation t-tests:') +print_results(np.c_[t_perm, p_perm], vlabels, 't') + +# Classical t-tests corrected by Bonferroni detect a trend in the difference +# between the third variables of samples, while permutation tests find a +# significant difference. + + +############################################################################### +# Hotelling's two-sample t-squared test, for multivariate samples +# --------------------------------------------------------------- + +# Hotelling's two-sample t-squared test is the multivariate generalization of +# the Student's t-test [1]. + +def hotelling(X, Y): + n_meas_X, n_meas_Y = X.shape[0], Y.shape[0] + n_vars = X.shape[1] + assert Y.shape[1] == n_vars + + n_meas = n_meas_X + n_meas_Y + diff = np.mean(X, axis=0) - np.mean(Y, axis=0) + Cov_X = np.cov(X, rowvar=False) + Cov_Y = np.cov(Y, rowvar=False) + Cov_pooled = ((n_meas_X-1) * Cov_X + (n_meas_Y-1) * Cov_Y) / (n_meas-2) + Cov_pooled_inv = np.linalg.pinv(Cov_pooled) + tsq = (n_meas_X*n_meas_Y) / n_meas * diff.T @ Cov_pooled_inv @ diff + + Fstat = tsq * (n_meas - n_vars - 1) / ((n_meas-2) * n_vars) + Fdist = f(n_vars, n_meas - 1 - n_vars) + pval = 1 - Fdist.cdf(Fstat) + + return tsq, pval + +tsq, p_multiv = hotelling(X, Y) +print('\nHotelling t-squared test:\n t^2={:.2f}, p={:.3e} ({})' + .format(tsq, p_multiv, pvals_to_stars(p_multiv))) + +# Hotelling's test detects a trend in the difference between samples too, and +# is not able to say on which dimension. + + +############################################################################### +# Reference +# --------- +# [1] Hotelling, "The Generalization of Student’s Ratio", Ann Math Statist, +# 1931. diff --git a/tests/test_core.py b/tests/test_core.py new file mode 100644 index 0000000..9c6e82e --- /dev/null +++ b/tests/test_core.py @@ -0,0 +1,102 @@ +""" Tests for module core. + +To execute tests: +>>> py.test -k test_core +""" + +import pytest +from itertools import combinations +import pypermut.core as core + + +@pytest.mark.parametrize("n", range(2, 9)) +def test_permutations_measurements(n): + """ Test permutations along measurements: + uniqueness and null permutation. + + Parameter + --------- + n : int + Length of samples that will be tested. Do not use values over 8, since + it would take too long to permute completely. + """ + assert n > 1, 'Need at least 2 values.' + + max_perms = core.count_permutation_measurements(n) + # Get all permutations + perms = [] + for p in range(max_perms): + perm = core.get_permutation_measurements(n, p) + perms.append(perm) + + # Uniqueness: there should be exactly max_perms permutations, + # and using a set will help determine if there are repeated results. + permutations = set(tuple(perm) for perm in perms) + assert len(permutations) == max_perms, 'Generated permutations are not unique.' + + # Null permutation + assert perms[0] == list(range(n)), 'First permutation is not the null permutation.' + + +@pytest.mark.parametrize("n", range(2, 15)) +def test_permutations_paired_samples(n): + """ Test permutations for paired samples: + uniqueness and null/full permutations. + + Parameter + --------- + n : int + Length of samples that will be tested. + """ + assert n > 1, 'Need at least 2 values.' + + max_perms = core.count_permutations_paired_samples(2, n) + # Get all permutations + perms = [] + for p in range(max_perms): + perm = core.get_permutation_2_paired_samples(n, p) + perms.append(perm) + + # Uniqueness: there should be exactly max_perms permutations, + # and using a set will help determine if there are repeated results. + permutations = set(tuple(map(tuple, perm)) for perm in perms) + assert len(permutations) == max_perms, 'Generated permutations are not unique.' + + # Null and full permutations + assert all(v == 1 for v in perms[0]), 'First permutation is not the null permutation.' + assert all(v == -1 for v in perms[-1]), 'Last permutation is not the full permutation.' + + +@pytest.mark.parametrize("n1", range(2, 14)) +@pytest.mark.parametrize("n2", range(2, 10)) +def test_permutations_unpaired_samples(n1, n2): + """ Test permutations for unpaired samples: + uniqueness and null permutation. + + Parameters + ---------- + n1 : int + Length of first sample that will be tested. + + n2 : int + Length of second sample that will be tested. + """ + assert n1 > 1, 'Need at least 2 values.' + assert n2 > 1, 'Need at least 2 values.' + + max_perms = core.count_permutations_unpaired_samples([n1, n2]) + # Get all combinations + combs = list(combinations(range(n1 + n2), n1)) + # Get all permutations + perms = [] + for p in range(max_perms): + perm = core.get_permutation_unpaired_samples([n1, n2], list(combs[p])) + perms.append(perm) + + # Uniqueness: there should be exactly max_perms permutations, + # and using a set will help determine if there are repeated results. + permutations = set(tuple(perm) for perm in perms) + assert len(permutations) == max_perms, 'Generated permutations are not unique.' + + # Null permutation + assert perms[0] == list(range(n1 + n2)), 'First permutation is not the null permutation.' diff --git a/tests/test_stats.py b/tests/test_stats.py new file mode 100644 index 0000000..d560b2c --- /dev/null +++ b/tests/test_stats.py @@ -0,0 +1,914 @@ +""" Tests for module stats. + +To execute tests: +>>> py.test -k test_stats +""" + +import warnings +import pytest +import numpy as np +import scipy.stats as stats +import pypermut.core as core +import pypermut.stats as ppstats + +np.random.seed(17) + + +@pytest.mark.parametrize("n", range(2, 9)) +@pytest.mark.parametrize("corr_func", ('pearson', 'spearman', 'spearman_fast', 'spearman_fast_nonunique')) +def test_permutation_measurements_best(n, corr_func): + """ This function tests the statistic r and the exact p-value p of + correlation tests, when no permutation can improve the statistic. + + Parameters + ---------- + n : int + Length of samples that will be tested. Do not use values over 8, since + it would take too long to permute completely. + + corr_func : string + Define the correlation type: + 'pearson' for the Pearson product-moment correlation coefficient r, + 'spearman' for the Spearman rank-order correlation coefficient rho. + """ + assert n > 1, 'Need at least 2 values.' + + # Silence warning in spearman_fast_non_unique + if corr_func == 'spearman_fast_nonunique': + warnings.simplefilter('ignore', UserWarning) + + # Create a vector that increases monotonically: 0,1,2,...,n-3,n-2,n-1 + # Consequently, no permutation can improve the correlation. + y = np.arange(n) + + # There are n! different permutations. + max_perms = core.count_permutation_measurements(n) + + # Use a correlation permutation test with all permutations: + # it should give a p-value of exactly 1/(n!) + r, p, dist = ppstats.permutation_corr(y, + n='all', + side='one', + corr=corr_func, + return_dist=True) + + assert dist.shape[0] == max_perms, 'Null distribution should contain n! statistics.' + assert np.isclose(r[0], 1., atol=1e-9), 'Best permutation case should give r of exactly 1.' + assert np.isclose(p[0], 1./max_perms, atol=1e-9), 'Best permutation case should give p-value of exactly 1/(n!).' + + +@pytest.mark.parametrize("n", range(2, 15)) +def test_permutation_paired_samples_best(n): + """ This function tests the exact p-value of tests for two related samples, + when no permutation can improve the statistic. + + Parameter + --------- + n : int + Length of samples that will be tested. + """ + assert n > 1, 'Need at least 2 values.' + + # Create a first random vector x, with positive values + x = np.random.randn(n) + 5 + x[x < 0] = 0 + + # Create a second vector y, a quasi-shift of x, but without overlap with x. + # Consequently, no permutation can improve the statistic. + y = x + 0.001 * np.random.randn(n) - 10 + + # There are 2^n different permutations. + max_perms = core.count_permutations_paired_samples(2, n) + + # Use a Student's permutation t-test for related samples with all possible + # permutations: it should give a p-value of exactly 1/(2^n) + t, p, dist = ppstats.permutation_ttest_rel(x, y, + n='all', + side='one', + return_dist=True) + assert dist.shape[0] == max_perms, 'Student t-test rel: null distribution should contain 2^n statistics.' + assert np.isclose(p[0], 1./max_perms, atol=1e-9), 'Student t-test rel: best permutation case should give p-value of exactly 1/(2^n).' + + # Use a Wilcoxon permutation T test for paired samples with all possible + # permutations: it should give a p-value of exactly 2/(2^n) + # (because T stat is computed on the absolute differences, so full + # permutation gives the same T as null permutation) + T, p, dist = ppstats.permutation_wilcoxon(x, y, + n='all', + zero_method='wilcox', + return_dist=True) + assert dist.shape[0] == max_perms, 'Wilcoxon T test: null distribution should contain 2^n statistics.' + assert np.isclose(p[0], 2./max_perms, atol=1e-9), 'Wilcoxon T test: best permutation case should give p-value of exactly 1/(2^n).' + + +@pytest.mark.parametrize("n1", range(2, 10)) +@pytest.mark.parametrize("n2", range(2, 10)) +def test_permutation_unpaired_samples_best(n1, n2): + """ This function tests the exact p-value of tests for two independent + samples, when no permutation can improve the statistic. + + Parameters + ---------- + n1 : int + Length of first sample that will be tested. + + n2 : int + Length of second sample that will be tested. + """ + assert n1 > 1, 'Need at least 2 values.' + assert n2 > 1, 'Need at least 2 values.' + + # Create a first random vector x, with positive values + x = np.random.randn(n1) + 5 + x[x < 0] = 0 + + # Create a second vector y, with negative values, without overlap with x. + # Consequently, no permutation can improve the statistic. + y = np.random.randn(n2) - 5 + y[y > 0] = 0 + + # There are (n1+n2)!/(n1!n2!) different permutations. + max_perms = core.count_permutations_unpaired_samples([n1, n2]) + + # Use a Student's permutation t-test for independent samples with all + # possible permutations: it should give a p-value of exactly + # 1/((n1+n2)!/(n1!n2!)) + t, p, dist = ppstats.permutation_ttest_ind(x, y, + n='all', + side='one', + return_dist=True) + assert dist.shape[0] == max_perms, 'Student t-test ind: null distribution should contain (n1+n2)!/(n1!n2!) statistics.' + # TODO: assert p-value + # NOT DONE, because percentileofscore randomly detects the presence of stat + # of null permutation in vector dist... + #assert np.isclose(p[0], 1./max_perms, atol=1e-9), 'Student t-test ind: best permutation case should give p-value of exactly 1/((n1+n2)!/(n1!n2!)).' + + # Use a Mann-Whitney permutation U test for unpaired samples with all + # possible permutations: it should give a p-value of exactly + # 2/((n1+n2)!/(n1!n2!)) + # (because computing U = min(U_x, U_y) gives two identical U) + U, p, dist = ppstats.permutation_mannwhitneyu(x, y, + n='all', + return_dist=True) + assert dist.shape[0] == max_perms, 'Mann-Whitney U test: null distribution should contain (n1+n2)!/(n1!n2!) statistics.' + assert np.isclose(p[0], 2./max_perms, atol=1e-9), 'Mann-Whitney U test: best permutation case should give p-value of exactly 1/((n1+n2)!/(n1!n2!)).' + + +@pytest.mark.parametrize("corr_func", ('pearson', 'spearman', 'spearman_fast', 'spearman_fast_nonunique')) +def test_permutation_corr(corr_func, n_reps=20, n_meas=25, n_vars=10, + trend_amplitude=3, side='one', alpha=0.05): + """ This function validates the permutation_corr function by generating + random samples on which statistics and p-values are computed. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + n_meas : int, optional + Number of measurements in samples X and Y. + + n_vars : int, optional + Number of variables in samples X and Y. + + trend_amplitude : float, optional + Amplitude of the generated temporal trend over unitary flat noise. + + alpha : float, optional + Type-I error used for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + # Silence warning in spearman_fast_non_unique + if corr_func == 'spearman_fast_nonunique': + warnings.simplefilter('ignore', UserWarning) + + ### check dimensions + ppstats.permutation_corr(np.random.randn(n_meas), + x=np.random.randn(n_meas), + n=1) + ppstats.permutation_corr(np.random.randn(n_meas, n_vars), + x=np.random.randn(n_meas), + n=1) + with pytest.raises(ValueError): + ppstats.permutation_corr(np.random.randn(n_meas, n_vars), + x=np.random.randn(n_meas+1), + n=1) # not same n_meas + ppstats.permutation_corr(np.random.randn(n_meas, n_vars), + x=np.random.randn(n_meas, n_vars), + n=1) # x not univariate + ppstats.permutation_corr(np.random.randn(n_meas, n_vars, 3), + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + # generate a random data matrix + Y = np.random.random((n_meas, n_vars)) + # add tendency only to the last variable: + # this trend is negative, so should be captured by one-sided test + Y[:, -1] -= trend_amplitude * np.arange(n_meas, dtype=float) / n_meas + Y = -Y + + R_pp, pvals = ppstats.permutation_corr(Y, + n=100, + with_replacement=False, + side=side, + corr=corr_func) + + # assert r values with respect scipy + x = np.arange(n_meas) + R_sp = np.zeros(n_vars) + for v, y in enumerate(Y.T): + if 'pearson' in corr_func: + r, _ = stats.pearsonr(x, y) + R_sp[v] = r + else: + R_sp[v] = stats.spearmanr(x, y).correlation + assert np.allclose(R_pp, R_sp, atol=1e-9), 'Statistics R should be equivalent to scipy.' + + # the trend was only added to the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +def test_permutation_ttest_rel(n_reps=50, n_meas=100, n_vars=10, + diff_mean=-2, diff_std=0.8, alpha=0.05): + """ This function validates the permutation_ttest_rel function by + generating random samples on which statistics and p-values are computed. + + H0: difference between pairs follows a Gaussian distribution. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + n_meas : int, optional + Number of measurements in samples X and Y. + + n_vars : int, optional + Number of variables in samples X and Y. + + diff_mean : float, optional + Value of the mean of the distribution in Y different from X. + + diff_std : float, optional + Value of the standard deviation of the distribution in Y different from + X. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + ### check dimensions + ppstats.permutation_ttest_rel(np.random.randn(n_meas), + np.random.randn(n_meas), + n=1) + ppstats.permutation_ttest_rel(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas, n_vars), + n=1) + with pytest.raises(ValueError): + ppstats.permutation_ttest_rel(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas+1, n_vars), + n=1) # not same n_meas + ppstats.permutation_ttest_rel(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas, n_vars+1), + n=1) # not same n_vars + ppstats.permutation_ttest_rel(np.random.randn(n_meas, n_vars, 3), + np.random.randn(n_meas, n_vars, 3), + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + # X = random Gaussian data, mean=0, std=1 + X = np.random.randn(n_meas, n_vars) + # Y = X + Gaussian noise, to be compliant with H0 + # ie, pairs difference follows a Gaussian distribution around zero + # because should give a low value t + Y = X + 0.1 * np.random.randn(n_meas, n_vars) + # add an offset to last variable, to reject H0 + # because should give a high value t + Y[:, -1] = diff_std * np.random.randn(n_meas) + diff_mean + + t_pp, pvals = ppstats.permutation_ttest_rel(X, Y, + n=100, + with_replacement=False, + side='one') + + # assert t values with respect scipy + t_sp = np.zeros(n_vars) + for v, (x, y) in enumerate(zip(X.T, Y.T)): + t_sp[v] = stats.ttest_rel(x, y).statistic + assert np.allclose(t_pp, t_sp, atol=1e-9), 'Statistics t should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +@pytest.mark.parametrize("equal_var", (True, False)) +def test_permutation_ttest_ind(equal_var, n_reps=50, n_meas_X=100, + n_meas_Y=110, n_vars=10, + diff_mean=-2, alpha=0.05): + """ This function validates the permutation_ttest_ind function by + generating random samples on which statistics and p-values are computed. + + H0: samples are drawn from the same Gaussian distribution, ie with equal + means and variances (excepted for Welch's version). + + Parameters + ---------- + equal_var : bool + If True, it performs the standard independent two samples test that + assumes equal variances. + If False, it performs the Welch's t-test, which does not assume equal + variance. + + n_reps : int, optional + Number of repetitions. + + n_meas_X : int, optional + Number of measurements in X. + + n_meas_Y : int, optional + Number of measurements in Y. + + n_vars : int, optional + Number of variables in X and Y. + + diff_mean : float, optional + Value of the mean of the distribution in Y different from X. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate + """ + ### check dimensions + ppstats.permutation_ttest_ind(np.random.randn(n_meas_X), + np.random.randn(n_meas_Y), + n=1) + ppstats.permutation_ttest_ind(np.random.randn(n_meas_X, n_vars), + np.random.randn(n_meas_Y, n_vars), + n=1) + with pytest.raises(ValueError): + ppstats.permutation_ttest_ind(np.random.randn(n_meas_X, n_vars), + np.random.randn(n_meas_Y, n_vars+1), + n=1) # not same n_vars + ppstats.permutation_ttest_ind(np.random.randn(n_meas_X, n_vars, 3), + np.random.randn(n_meas_Y, n_vars, 3), + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + if equal_var: + std = 1 + else: + # std = random uniform value between 1 and 10 + std = 10 * np.random.rand(n_vars) + 1 + + for i in range(n_reps): + # X = random Gaussian data, mean=0, std=1 + X = np.random.randn(n_meas_X, n_vars) + # Y = random Gaussian data, mean=0, std=1 if equal_var=True + # or std=random number if equal_var=False, + # to be compliant with H0, because should give a low value t + Y = std * np.random.randn(n_meas_Y, n_vars) + # add an offset to last variable, to reject H0 + # because should give a high value t + Y[:, -1] = np.random.randn(n_meas_Y) + diff_mean + + t_pp, pvals = ppstats.permutation_ttest_ind(X, Y, + n=100, + with_replacement=False, + side='one', + equal_var=equal_var) + + # assert t values with respect scipy + t_sp = np.zeros(n_vars) + for v, (x, y) in enumerate(zip(X.T, Y.T)): + t_sp[v] = stats.ttest_ind(x, y, equal_var=equal_var).statistic + assert np.allclose(t_pp, t_sp, atol=1e-9), 'Statistics t should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +@pytest.mark.parametrize("zero_method", ('pratt', 'wilcox', 'zsplit')) +def test_permutation_wilcoxon(zero_method, n_reps=50, n_meas=100, + n_vars=10, diff_mean=2, diff_std=1, + alpha=0.05): + """ This function validates the permutation_wilcoxon function by generating + random samples on which statistics and p-values are computed. + + H0: difference between pairs follows a symmetric distribution around zero. + H1: difference between pairs does not follow a symmetric distribution + around zero. + + Parameters + ---------- + zero_method : {'pratt', 'wilcox', 'zsplit'} + Method for zero-differences processing. + + n_reps : int, optional + Number of repetitions. + + n_meas : int, optional + Number of measurements in samples X and Y. + + n_vars : int, optional + Number of variables in samples X and Y. + + diff_mean : float, optional + Value of the mean of the distribution in Y different from X. + + diff_std : float, optional + Value of the standard deviation of the distribution in Y different from + X. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate + """ + ### check dimensions + ppstats.permutation_wilcoxon(np.random.randn(n_meas), + np.random.randn(n_meas), + n=1) + ppstats.permutation_wilcoxon(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas, n_vars), + n=1) + with pytest.raises(ValueError): + ppstats.permutation_wilcoxon(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas+1, n_vars), + n=1) # not same n_meas + ppstats.permutation_wilcoxon(np.random.randn(n_meas, n_vars), + np.random.randn(n_meas, n_vars+1), + n=1) # not same n_vars + ppstats.permutation_wilcoxon(np.random.randn(n_meas, n_vars, 3), + np.random.randn(n_meas, n_vars, 3), + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + # X = random uniform data between 0 and 1 + X = np.random.rand(n_meas, n_vars) + # Y = X + Gaussian noise, to be compliant with H0 + # ie, pairs difference follows a symmetric distribution around zero + # because should give a high value T + Y = X + 0.01 * np.random.randn(n_meas, n_vars) + # add an offset to last variable, to reject H0 + # because should give a low value T + Y[:, -1] = diff_std * np.random.randn(n_meas) + diff_mean + + T_pp, pvals = ppstats.permutation_wilcoxon(X, Y, + n=100, + with_replacement=False, + zero_method=zero_method) + + # assert T values with respect scipy + T_sp = np.zeros(n_vars) + for v, (x, y) in enumerate(zip(X.T, Y.T)): + T_sp[v] = stats.wilcoxon(x, y, zero_method=zero_method).statistic + assert np.allclose(T_pp, T_sp, atol=1e-9), 'Statistics T should be equivalent to scipy.' + + # the offset was only added to the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +def test_permutation_mannwhitneyu(n_reps=50, n_meas_X=100, n_meas_Y=110, + n_vars=10, diff_mean=1, diff_std=1, + alpha=0.05): + """ This function validates the permutation_mannwhitneyu function by + generating random samples on which statistics and p-values are computed. + + H0: samples are drawn from the same distribution, that may not be normal. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + n_meas_X : int, optional + Number of measurements in X. + + n_meas_Y : int, optional + Number of measurements in Y. + + n_vars : int, optional + Number of variables in X and Y. + + diff_mean : float, optional + Value of the mean of the distribution in Y different from X. + + diff_std : float, optional + Value of the standard deviation of the distribution in Y different from + X. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + ### check dimensions + ppstats.permutation_mannwhitneyu(np.random.randn(n_meas_X), + np.random.randn(n_meas_Y), + n=1) + ppstats.permutation_mannwhitneyu(np.random.randn(n_meas_X, n_vars), + np.random.randn(n_meas_Y, n_vars), + n=1) + with pytest.raises(ValueError): + ppstats.permutation_mannwhitneyu(np.random.randn(n_meas_X, n_vars), + np.random.randn(n_meas_Y, n_vars+1), + n=1) # not same n_vars + ppstats.permutation_mannwhitneyu(np.random.randn(n_meas_X, n_vars, 3), + np.random.randn(n_meas_Y, n_vars, 3), + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + # X = random Gaussian data, mean=0, std=1 + X = np.random.randn(n_meas_X, n_vars) + # Y = random Gaussian data, mean=0, std=1, to be compliant with H0 + # because should give a high value U + Y = np.random.randn(n_meas_Y, n_vars) + # add an offset to last variable, to reject H0 + # because should give a low value U + Y[:, -1] = diff_std * np.random.randn(n_meas_Y) + diff_mean + + U_pp, pvals = ppstats.permutation_mannwhitneyu(X, Y, + n=100, + with_replacement=False) + + # TODO: assert U values with respect scipy + # NOT DONE, because scipy implementation is really "different" ... +# U_sp = np.zeros(n_vars) +# for v, (x, y) in enumerate(zip(X.T, Y.T)): +# U_sp[v] = stats.mannwhitneyu(x, y, alternative='less').statistic +# assert np.allclose(U_pp, U_sp, atol=1e-9), 'Statistics U should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +def test_permutation_f_oneway(n_reps=50, list_meas=[80, 100, 50, 70], + n_vars=3, diff_mean=1, diff_std=1, + alpha=0.05): + """ This function validates the permutation_f_oneway function by + generating random samples on which statistics and p-values are computed. + + H0: means of Gaussian samples are all equal. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + list_meas : list of int, optional + List of number of measurements of each sample. + Its length defines the number of samples / groups. + + n_vars : int, optional + Number of variables in samples. + + diff_mean : float, optional + Value of the mean of the distribution in the last variable different + from others. + + diff_std : float, optional + Value of the standard deviation of the distribution in the last + variable different from others. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + ### check dimensions + ppstats.permutation_f_oneway( + *[np.random.randn(n_meas) for n_meas in list_meas], + n=1) + ppstats.permutation_f_oneway( + *[np.random.randn(n_meas, n_vars) for n_meas in list_meas], + n=1) + with pytest.raises(ValueError): + ppstats.permutation_f_oneway(*[np.random.randn(list_meas[0], n_vars)], + n=1) # only one group + list_vars = np.random.randint(1, 10, size=len(list_meas)) + ppstats.permutation_f_oneway( + *[np.random.randn(n_meas, n_vars_) + for n_meas, n_vars_ in zip(list_meas, list_vars)], + n=1) # not same n_vars + ppstats.permutation_f_oneway( + *[np.random.randn(n_meas, n_vars, 3) for n_meas in list_meas], + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + args = [] + for j, n_meas in enumerate(list_meas): + # all samples = random Gaussian data, mean=0, std=1 + # to be compliant with H0, because should give a low value F + X = np.random.randn(n_meas, n_vars) + if j == len(list_meas)-1: + # add an offset to last sample of the last variable, + # to reject H0, because should give a high value F + X[:, -1] = diff_std * np.random.randn(n_meas) + diff_mean + args.append(X) + + F_pp, pvals = ppstats.permutation_f_oneway(*args, n=500) + + # assert F values with respect scipy + F_sp = np.zeros(n_vars) + for v in range(n_vars): + args_ = [arg[:, v] for arg in args] + F_sp[v] = stats.f_oneway(*args_).statistic + assert np.allclose(F_pp, F_sp, atol=1e-9), 'Statistics F should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +def test_permutation_kruskal(n_reps=50, list_meas=[80, 100, 50, 70], + n_vars=3, diff_mean=1, diff_std=1, + alpha=0.05): + """ This function validates the permutation_kruskal function by + generating random samples on which statistics and p-values are computed. + + H0: medians of samples are all equal. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + list_meas : list of int, optional + List of number of measurements of each sample. + Its length defines the number of samples / groups. + + n_vars : int, optional + Number of variables in samples. + + diff_mean : float, optional + Value of the mean of the distribution in the last variable different + from others. + + diff_std : float, optional + Value of the standard deviation of the distribution in the last + variable different from others. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + ### check dimensions + ppstats.permutation_kruskal( + *[np.random.randn(n_meas) for n_meas in list_meas], + n=1) + ppstats.permutation_kruskal( + *[np.random.randn(n_meas, n_vars) for n_meas in list_meas], + n=1) + with pytest.raises(ValueError): + ppstats.permutation_kruskal(*[np.random.randn(list_meas[0], n_vars)], + n=1) # only one group + list_vars = np.random.randint(1, 10, size=len(list_meas)) + ppstats.permutation_kruskal( + *[np.random.randn(n_meas, n_vars_) + for n_meas, n_vars_ in zip(list_meas, list_vars)], + n=1) # not same n_vars + ppstats.permutation_kruskal( + *[np.random.randn(n_meas, n_vars, 3) for n_meas in list_meas], + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + args = [] + for j, n_meas in enumerate(list_meas): + # all samples = random Gaussian data, mean=0, std=1 + # to be compliant with H0, because should give a low value H + X = np.random.randn(n_meas, n_vars) + if j == len(list_meas)-1: + # add an offset to last sample of the last variable, + # to reject H0, because should give a high value H + X[:, -1] = diff_std * np.random.randn(n_meas) + diff_mean + args.append(X) + + H_pp, pvals = ppstats.permutation_kruskal(*args, n=500) + + # assert H values with respect scipy + H_sp = np.zeros(n_vars) + for v in range(n_vars): + args_ = [arg[:, v] for arg in args] + H_sp[v] = stats.kruskal(*args_).statistic + assert np.allclose(H_pp, H_sp, atol=1e-9), 'Statistics H should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha + + +def test_permutation_friedmanchisquare(n_reps=50, n_meas=10, n_vars=3, + n_groups=5, diff_mean=3, diff_std=0.5, + alpha=0.05): + """ This function validates the permutation_friedmanchisquare function by + generating random samples on which statistics and p-values are computed. + + H0: medians of samples are all equal. + + Parameters + ---------- + n_reps : int, optional + Number of repetitions. + + n_meas : int, optional + Number of measurements in samples. + + n_vars : int, optional + Number of variables in samples. + + n_groups : int, optional + Number of samples / groups. + + diff_mean : float, optional + Value of the mean of the distribution in the last variable different + from others. + + diff_std : float, optional + Value of the standard deviation of the distribution in the last + variable different from others. + + alpha : float, optional + Type-I error use for the generation of false and true positive rates. + + Returns + ------- + fpr : float + False positive rate. + + tpr : float + True positive rate. + """ + ### check dimensions + ppstats.permutation_friedmanchisquare( + *[np.random.randn(n_meas, n_vars) for g in range(n_groups)], + n=1) + with pytest.raises(ValueError): + ppstats.permutation_friedmanchisquare( + *[np.random.randn(n_meas) for g in range(n_groups)], + n=1) # not 2D inputs # TODO + ppstats.permutation_friedmanchisquare( + *[np.random.randn(n_meas) for g in range(2)], + n=1) # only two groups + list_vars = np.random.randint(1, 10, size=len(n_groups)) + ppstats.permutation_friedmanchisquare( + *[np.random.randn(n_meas, n_vars_) for n_vars_ in list_vars], + n=1) # not same n_vars + ppstats.permutation_friedmanchisquare( + *[np.random.randn(n_meas, n_vars, 3) for g in range(n_groups)], + n=1) # more than 2 dims + + ### check false positive and true positive rates + false_positives = [] + true_positives = [] + + for i in range(n_reps): + args = [] + for j in range(n_groups): + # all samples = random Gaussian data, mean=0, std=1 + # to be compliant with H0, because should give a low value chi2 + X = np.random.randn(n_meas, n_vars) + if j == n_groups-1: + # add an offset to last sample of the last variable, + # to reject H0, because should give a high value chi2 + X[:, -1] = diff_std * np.random.randn(n_meas) + diff_mean + args.append(X) + + chi2_pp, pvals = ppstats.permutation_friedmanchisquare(*args, n=500) + + # assert chi2 values with respect scipy + chi2_sp = np.zeros(n_vars) + for v in range(n_vars): + args_ = [arg[:, v] for arg in args] + chi2_sp[v] = stats.friedmanchisquare(*args_).statistic + assert np.allclose(chi2_pp, chi2_sp, atol=1e-9), 'Statistics chi2 should be equivalent to scipy.' + + # the distribution is different only for the last index + false_positives.append(sum([int(p < alpha) for p in pvals[:-1]])) + true_positives.append(int(pvals[-1] < alpha)) + + # assert false positive and true positive rates + fpr = sum(false_positives) / ((n_vars - 1) * n_reps) + tpr = sum(true_positives) / n_reps + assert fpr < alpha + assert tpr > 1. - alpha +