diff --git a/README.md b/README.md index 112a3cddc..ad3e54cd5 100644 --- a/README.md +++ b/README.md @@ -140,6 +140,5 @@ Contributions of any kind are welcome! ## Project History hyppo is a rebranding of mgcpy, which was founded in November 2018. -mgcpy was designed and written by Satish Palaniappan, Sambit -Panda, Junhao Xiong, Sandhya Ramachandran, and Ronak Mehtra. hyppo +mgcpy was designed and written by @tpsatish95, @sampan501, @junhaobearxiong, @sundaysundya, @ananyas713, and @ronakdm. hyppo was designed and written by Sambit Panda and first released in February 2020. diff --git a/benchmarks/condi_indep_power_sampsize.py b/benchmarks/condi_indep_power_sampsize.py new file mode 100644 index 000000000..f37843385 --- /dev/null +++ b/benchmarks/condi_indep_power_sampsize.py @@ -0,0 +1,132 @@ +""" +1D Independence Testing Power vs. Sample Size +=============================================== + +Here, we show finite testing power comparisons between the various tests within hyppo. +For a test to be consistent, we would expect power to converge to 1 as sample size +increases. Tests that converge to 1 quicker have higher finite testing power and +are likely better to use for your use case. The simulation in the bottom right is +used so that we know that we are properly controlling for type I error, which is +important becase otherwise the test would be invalid (power = alpha-level = 0.05). +""" + +import os +import sys + +import matplotlib.pyplot as plt +import numpy as np +import seaborn as sns +from hyppo.conditional import COND_INDEP_TESTS +from hyppo.tools import COND_SIMULATIONS, power +from joblib import Parallel, delayed + +sys.path.append(os.path.realpath("..")) + +# make plots look pretty +sns.set(color_codes=True, style="white", context="talk", font_scale=2) +PALETTE = sns.color_palette("Set1") +sns.set_palette(PALETTE[1:]) + +# constants +MAX_SAMPLE_SIZE = 100 +STEP_SIZE = 5 +SAMP_SIZES = range(5, MAX_SAMPLE_SIZE + STEP_SIZE, STEP_SIZE) +POWER_REPS = 5 + +# simulation titles +SIM_TITLES = [k for k, v in COND_SIMULATIONS.items()] + +# these tests only make sense for > 1 dimension data +remove = ["fcit", "kci"] +COND_INDEP_TESTS = dict([(k, v) for k, v in INDEP_TESTS.items() if k not in remove]) + + +def estimate_power(sim, test, auto=False): + """Compute the mean of the estimated power of 5 replications over sample sizes.""" + if test == "MaxMargin": + test = ["MaxMargin", "Dcorr"] + est_power = np.array( + [ + np.mean( + [ + power(test, pow_type="indep", sim=sim, n=i, p=1, auto=auto) + for _ in range(POWER_REPS) + ] + ) + for i in SAMP_SIZES + ] + ) + np.savetxt( + "../benchmarks/conditional_vs_samplesize/{}_{}.csv".format(sim, test), + est_power, + delimiter=",", + ) + + return est_power + + +# At this point, we would run this bit of code to generate the data for the figure and +# store it under the "vs_sampsize" directory. Since this code takes a very long time, +# we have commented out these lines of codes. If you would like to generate the data, +# uncomment these lines and run the file. + +outputs = Parallel(n_jobs=-1, verbose=100)( + [ + delayed(estimate_featimport)(sim_name, test) + for sim_name in SIMULATIONS.keys() + for test in INDEP_TESTS.keys() + ] +) + + +def plot_power(): + fig, ax = plt.subplots(nrows=4, ncols=5, figsize=(25, 20)) + plt.suptitle( + "Multivariate Independence Testing (Increasing Sample Size)", + y=0.93, + va="baseline", + ) + + for i, row in enumerate(ax): + for j, col in enumerate(row): + count = 5 * i + j + sim = list(SIMULATIONS.keys())[count] + + for test in INDEP_TESTS.keys(): + est_power = np.genfromtxt( + "../benchmarks/vs_samplesize/{}_{}.csv".format(sim, test), + delimiter=",", + ) + + col.plot(SAMP_SIZES, est_power, label=INDEP_TESTS[test].__name__, lw=2) + col.set_xticks([]) + if i == 3: + col.set_xticks([SAMP_SIZES[0], SAMP_SIZES[-1]]) + col.set_ylim(-0.05, 1.05) + col.set_yticks([]) + if j == 0: + col.set_yticks([0, 1]) + col.set_title(SIM_TITLES[count]) + + fig.text(0.5, 0.05, "Sample Size", ha="center") + fig.text( + 0.07, + 0.5, + "Statistical Power", + va="center", + rotation="vertical", + ) + leg = plt.legend( + bbox_to_anchor=(0.5, 0.05), + bbox_transform=plt.gcf().transFigure, + ncol=len(INDEP_TESTS.keys()), + loc="upper center", + ) + leg.get_frame().set_linewidth(0.0) + for legobj in leg.legendHandles: + legobj.set_linewidth(5.0) + plt.subplots_adjust(hspace=0.50) + + +# plot the power +plot_power() diff --git a/docs/changelog/v0.3.3.md b/docs/changelog/v0.3.3.md new file mode 100644 index 000000000..1eec255cb --- /dev/null +++ b/docs/changelog/v0.3.3.md @@ -0,0 +1,17 @@ +# hyppo v0.3.3 + +> **Note:** hyppo v0.3.3 has not been released yet! + +This is a minor release for general bug fixes, documentation improvements, and general package maintenance. + +## Authors + + + + + +## Issues Closed + +## PRs Merged + +* [#232](https://github.com/neurodata/hyppo/pull/310): run pytest in parallel diff --git a/hyppo/conditional/__init__.py b/hyppo/conditional/__init__.py index ff876ea3a..75aad2862 100644 --- a/hyppo/conditional/__init__.py +++ b/hyppo/conditional/__init__.py @@ -1,7 +1,15 @@ +from .cdcorr import ConditionalDcorr from .FCIT import FCIT from .kci import KCI - +from .pcorr import PartialCorr +from .pdcorr import PartialDcorr __all__ = [s for s in dir()] # add imported tests to __all__ -COND_INDEP_TESTS = {"fcit": FCIT, "kci": KCI} +COND_INDEP_TESTS = { + "fcit": FCIT, + "kci": KCI, + "conditionaldcorr": ConditionalDcorr, + "partialcorr": PartialCorr, + "partialdcorr": PartialDcorr, +} diff --git a/hyppo/conditional/_utils.py b/hyppo/conditional/_utils.py new file mode 100644 index 000000000..9d78c1403 --- /dev/null +++ b/hyppo/conditional/_utils.py @@ -0,0 +1,88 @@ +import numpy as np + +from ..tools import check_ndarray_xyz, check_reps, contains_nan, convert_xyz_float64 + + +class _CheckInputs: + """Checks inputs for all independence tests""" + + def __init__(self, x, y, z, reps=None, max_dims=None): + self.x = x + self.y = y + self.z = z + self.reps = reps + self.max_dims = max_dims + + def __call__(self): + check_ndarray_xyz(self.x, self.y, self.z) + contains_nan(self.x) + contains_nan(self.y) + contains_nan(self.z) + self.x, self.y, self.z = self.check_dim_xyz(max_dims=self.max_dims) + self.x, self.y, self.z = convert_xyz_float64(self.x, self.y, self.z) + self._check_min_samples() + self._check_variance() + + if self.reps: + check_reps(self.reps) + + return self.x, self.y, self.z + + def check_dim_xyz(self, max_dims): + """Check and convert x and y to proper dimensions""" + if self.x.ndim == 1: + self.x = self.x[:, np.newaxis] + elif self.x.ndim != 2: + raise ValueError( + "Expected a 2-D array `x`, found shape " "{}".format(self.x.shape) + ) + if self.y.ndim == 1: + self.y = self.y[:, np.newaxis] + elif self.y.ndim != 2: + raise ValueError( + "Expected a 2-D array `y`, found shape " "{}".format(self.y.shape) + ) + if self.z.ndim == 1: + self.z = self.z[:, np.newaxis] + elif self.z.ndim != 2: + raise ValueError( + "Expected a 2-D array `z`, found shape " "{}".format(self.z.shape) + ) + + if max_dims is not None: + _, dx = self.x.shape + _, dy = self.y.shape + _, dz = self.z.shape + + if np.any(np.array([dx, dy, dz]) > max_dims): + raise ValueError( + f"x, y, z must have be univariate and have shape [n,{max_dims}]" + ) + + self._check_nd_indeptest() + + return self.x, self.y, self.z + + def _check_nd_indeptest(self): + """Check if number of samples is the same""" + nx, _ = self.x.shape + ny, _ = self.y.shape + nz, _ = self.z.shape + if not np.all(np.array([nx, ny, nz]) == nx): + raise ValueError( + "Shape mismatch, x, y and z must have shape " + + "[n, p], [n, q] and [n, r]." + ) + + def _check_min_samples(self): + """Check if the number of samples is at least 3""" + nx = self.x.shape[0] + ny = self.y.shape[0] + nz = self.z.shape[0] + + if nx <= 3 or ny <= 3 or nz <= 3: + raise ValueError("Number of samples is too low") + + def _check_variance(self): + if np.var(self.x) == 0 or np.var(self.y) == 0 or np.var(self.z) == 0: + raise ValueError(f"Test cannot be run, one of the inputs has 0 variance {np.var(self.x)}, {np.var(self.y)}, {np.var(self.z)}, {self.z}, {self.z.shape}") diff --git a/hyppo/conditional/base.py b/hyppo/conditional/base.py index 102e24aff..3c3ef0282 100644 --- a/hyppo/conditional/base.py +++ b/hyppo/conditional/base.py @@ -13,7 +13,11 @@ class ConditionalIndependenceTest(ABC): """ - def __init__(self): + def __init__(self, **kwargs): + self.stat = None + self.pvalue = None + self.kwargs = kwargs + super().__init__() @abstractmethod diff --git a/hyppo/conditional/cdcorr.py b/hyppo/conditional/cdcorr.py new file mode 100644 index 000000000..cb6a758ad --- /dev/null +++ b/hyppo/conditional/cdcorr.py @@ -0,0 +1,291 @@ +from functools import partial + +import numpy as np +from scipy.spatial.distance import pdist, squareform + +from ..tools import compute_dist, perm_test +from ._utils import _CheckInputs +from .base import ConditionalIndependenceTest, ConditionalIndependenceTestOutput + + +class ConditionalDcorr(ConditionalIndependenceTest): + r""" + Conditional Distance Covariance/Correlation (CDcov/CDcorr) test statistic and p-value. + + CDcorr is a measure of dependence between two paired random matrices + given a third random matrix of not necessarily equal dimensions :footcite:p:`wang2015conditional`. + The coefficient is 0 if and only if the matrices are independent given + third matrix. + + Parameters + ---------- + compute_distance : str, callable, or None, default: "euclidean" + A function that computes the distance among the samples within each + data matrix. + Valid strings for ``compute_distance`` are, as defined in + :func:`sklearn.metrics.pairwise_distances`, + + - From scikit-learn: [``"euclidean"``, ``"cityblock"``, ``"cosine"``, + ``"l1"``, ``"l2"``, ``"manhattan"``] See the documentation for + :mod:`scipy.spatial.distance` for details + on these metrics. + - From scipy.spatial.distance: [``"braycurtis"``, ``"canberra"``, + ``"chebyshev"``, ``"correlation"``, ``"dice"``, ``"hamming"``, + ``"jaccard"``, ``"kulsinski"``, ``"mahalanobis"``, ``"minkowski"``, + ``"rogerstanimoto"``, ``"russellrao"``, ``"seuclidean"``, + ``"sokalmichener"``, ``"sokalsneath"``, ``"sqeuclidean"``, + ``"yule"``] See the documentation for :mod:`scipy.spatial.distance` for + details on these metrics. + + Set to ``None`` or ``"precomputed"`` if ``x`` and ``y`` are already distance + matrices. To call a custom function, either create the distance matrix + before-hand or create a function of the form ``metric(x, **kwargs)`` + where ``x`` is the data matrix for which pairwise distances are + calculated and ``**kwargs`` are extra arguements to send to your custom + function. + use_cov : bool, + If `True`, then the statistic will compute the covariance rather than the + correlation. + bandwith : str, scalar, 1d-array + The method used to calculate the bandwidth used for kernel density estimate of + the conditional matrix. This can be ‘scott’, ‘silverman’, a scalar constant or a + 1d-array with length ``r`` which is the dimensions of the conditional matrix. + If None (default), ‘scott’ is used. + **kwargs + Arbitrary keyword arguments for ``compute_distance``. + + References + ---------- + .. footbibliography:: + """ + + def __init__( + self, compute_distance="euclidean", use_cov=True, bandwidth=None, **kwargs + ): + self.use_cov = use_cov + self.compute_distance = compute_distance + self.bandwidth = bandwidth + + # Check bandwidth input + if bandwidth is not None: + if not isinstance(bandwidth, (int, float, np.ndarray, str)): + raise ValueError( + "`bandwidth` should be 'scott', 'silverman', a scalar or 1d-array." + ) + if isinstance(bandwidth, str): + if bandwidth.lower() not in ["scott", "silverman"]: + raise ValueError( + f"`bandwidth` must be either 'scott' or 'silverman' not '{bandwidth}'" + ) + self.bandwidth = bandwidth.lower() + + # set is_distance to true if compute_distance is None + self.is_distance = False + if not compute_distance: + self.is_distance = True + + ConditionalIndependenceTest.__init__(self, **kwargs) + + def __repr__(self): + return "ConditionalDcorr" + + def statistic(self, x, y, z): + r""" + Helper function that calculates the CDcov/CDcorr test statistic. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x`` and ``y`` can be + distance matrices and ``z`` can be a similarity matrix where the + shapes must be ``(n, n)``. + + Returns + ------- + stat : float + The computed CDcov/CDcorr statistic. + """ + check_input = _CheckInputs(x, y, z) + x, y, z = check_input() + + if not self.is_distance: + distx, disty = compute_dist( + x, + y, + metric=self.compute_distance, + ) + distz = self._compute_kde(z) + else: + distx = x + disty = y + distz = z + + if self.use_cov: + stat = _cdcov(distx, disty, distz).mean() + else: + stat = _cdcorr(distx, disty, distz).mean() + + self.stat = stat + return stat + + def test( + self, + x, + y, + z, + reps=1000, + workers=1, + random_state=None, + ): + r""" + Calculates the CDcov/CDcorr test statistic and p-value. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x`` and ``y`` can be + distance matrices and ``z`` can be a similarity matrix where the + shapes must be ``(n, n)``. + reps : int, default: 1000 + The number of replications used to estimate the null distribution + when using the permutation test used to calculate the p-value. + workers : int, default: 1 + The number of cores to parallelize the p-value computation over. + Supply ``-1`` to use all cores available to the Process. + random_state : int, default: None + The random_state for permutation testing to be fixed for + reproducibility. + + Returns + ------- + stat : float + The computed CDcov/CDcorr statistic. + pvalue : float + The computed CDcov/CDcorr p-value. + """ + check_input = _CheckInputs(x, y, z, reps=reps) + x, y, z = check_input() + + if not self.is_distance: + x, y = compute_dist(x, y, metric=self.compute_distance, **self.kwargs) + z = self._compute_kde(z) + + self.is_distance = True + + stat, pvalue, null_dist = perm_test( + self.statistic, + x=x, + y=y, + z=z, + reps=reps, + workers=workers, + is_distsim=self.is_distance, + random_state=random_state, + permuter=partial(self._permuter, probs=z), + ) + self.stat = stat + self.pvalue = pvalue + self.null_dist = null_dist + + return ConditionalIndependenceTestOutput(stat, pvalue) + + def _compute_kde(self, data): + n, d = data.shape + + if isinstance(self.bandwidth, (int, float)): + bandwidth_ = np.repeat(self.bandwidth, d) + elif isinstance(self.bandwidth, np.ndarray): + if not d == self.bandwidth.size: + raise ValueError(f"`bandwidth` must be of length {d}.") + bandwidth_ = self.bandwidth + elif self.bandwidth is None or self.bandwidth == "scott": + factor = np.power(n, (-1.0 / (d + 4))) + stds = np.std(data, axis=0, ddof=1) + bandwidth_ = factor * stds + elif self.bandwidth == "silverman": + factor = (n * (d + 2) / 4.0) ** (-1.0 / (d + 4)) + stds = np.std(data, axis=0, ddof=1) + bandwidth_ = factor * stds + + # Compute constants + denom = np.power(2 * np.pi, d / 2) * np.power(bandwidth_.prod(), 0.5) + + # Compute weighted kernel + sim_mat = pdist( + data, "sqeuclidean", w=1 / bandwidth_**2 + ) # Section 4.2 equation + sim_mat = squareform(-0.5 * sim_mat) + np.exp(sim_mat, sim_mat) + sim_mat /= denom + + rng = np.random.default_rng() + sim_mat += rng.standard_normal(size=sim_mat.shape) * 1e-6 + return sim_mat + + def _permuter(self, probs, rng=None, axis=1): + """ + Weighted sampling with replacement + """ + if rng is None: + rng = np.random.default_rng() + + n = probs.shape[1 - axis] + sums = probs.sum(axis=1, keepdims=True) + idx = ((probs / sums).cumsum(axis=1) > rng.random(n)[:, None]).argmax(axis=1) + + return idx + + +def _weighted_center_distmat(distx, weights): + """Centers the distance matrices with weights""" + + n = distx.shape[0] + + scl = np.sum(weights) + row_sum = np.sum(np.multiply(distx, weights), axis=1) / scl + total_sum = weights @ row_sum / scl + + cent_distx = distx - row_sum.reshape(-1, n).T - row_sum.reshape(-1, n) + total_sum + + return cent_distx + + +def _cdcov(distx, disty, distz): + """Calculate the CDcov test statistic""" + n = distx.shape[0] + + cdcov = np.zeros(n) + + for i in range(n): + r = distz[[i]] + cdx = _weighted_center_distmat(distx, distz[i]) + cdy = _weighted_center_distmat(disty, distz[i]) + cdcov[i] = (cdx * cdy * r * r.T).sum() / r.sum() ** 2 + + cdcov *= 12 * np.power(distz.mean(axis=0), 4) + + return cdcov + + +def _cdcorr(distx, disty, distz): + """Calculate the CDcorr test statistic""" + varx = _cdcov(distx, distx, distz) + vary = _cdcov(disty, disty, distz) + covar = _cdcov(distx, disty, distz) + + # stat is 0 with negative variances (would make denominator undefined) + denom = varx * vary + non_positives = denom <= 0 + if np.any(non_positives): + denom[non_positives] = 0 + stat = np.divide( + covar, np.sqrt(denom), out=np.zeros_like(covar), where=np.invert(non_positives) + ) + + return stat diff --git a/hyppo/conditional/pcorr.py b/hyppo/conditional/pcorr.py new file mode 100644 index 000000000..f725473f1 --- /dev/null +++ b/hyppo/conditional/pcorr.py @@ -0,0 +1,179 @@ +import numpy as np +from scipy.stats import t + +from ..tools import perm_test +from ._utils import _CheckInputs +from .base import ConditionalIndependenceTest, ConditionalIndependenceTestOutput + + +class PartialCorr(ConditionalIndependenceTest): + r""" + Conditional Pearson's correlation test. + + Partial correlation is a measure of the association between two univariate + variables given a third univariate variable. + + Parameters + ---------- + compute_distance : str, callable, or None, default: "euclidean" + A function that computes the distance among the samples within each + data matrix. + Valid strings for ``compute_distance`` are, as defined in + :func:`sklearn.metrics.pairwise_distances`, + + - From scikit-learn: [``"euclidean"``, ``"cityblock"``, ``"cosine"``, + ``"l1"``, ``"l2"``, ``"manhattan"``] See the documentation for + :mod:`scipy.spatial.distance` for details + on these metrics. + - From scipy.spatial.distance: [``"braycurtis"``, ``"canberra"``, + ``"chebyshev"``, ``"correlation"``, ``"dice"``, ``"hamming"``, + ``"jaccard"``, ``"kulsinski"``, ``"mahalanobis"``, ``"minkowski"``, + ``"rogerstanimoto"``, ``"russellrao"``, ``"seuclidean"``, + ``"sokalmichener"``, ``"sokalsneath"``, ``"sqeuclidean"``, + ``"yule"``] See the documentation for :mod:`scipy.spatial.distance` for + details on these metrics. + + Set to ``None`` or ``"precomputed"`` if ``x`` and ``y`` are already distance + matrices. To call a custom function, either create the distance matrix + before-hand or create a function of the form ``metric(x, **kwargs)`` + where ``x`` is the data matrix for which pairwise distances are + calculated and ``**kwargs`` are extra arguements to send to your custom + function. + **kwargs + Arbitrary keyword arguments for ``compute_distance``. + + Notes + ----- + The statistic is computed as follows: + + ..math:: + r_{x, y ; z} = \frac{\rho_{xy} - \rho_{xz} \rho_{yz}}{\sqrt{(1 - \rho_{xz}^2)(1 - \rho_{yz}^2)}} + + where :math:`\rho_{xy}` is the Pearson correlation coefficient between :math:`x` and + :math:`y`. The partial correlation test is implemented as a t-test footcite:p:`legendre2000`:. + + References + ---------- + .. footbibliography:: + """ + + def __init__(self, **kwargs): + ConditionalIndependenceTest.__init__(self, **kwargs) + + def statistic(self, x, y, z): + r""" + Helper function that calculates the partial correlation test statistic. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x`` and ``y`` can be + distance matrices and ``z`` can be a similarity matrix where the + shapes must be ``(n, n)``. + + Returns + ------- + stat : float + The computed partial correlation test statistic. + """ + check_input = _CheckInputs(x, y, z, max_dims=1) + x, y, z = check_input() + + corr = np.corrcoef(np.hstack([x, y, z]), rowvar=False) + cov_xy, cov_xz, cov_yz = corr[np.triu_indices_from(corr, k=1)] + + if cov_xz <= 0 or cov_yz <= 0: + stat = 0 + else: + stat = (cov_xy - cov_xz * cov_yz) / np.sqrt( + (1 - cov_xz**2) * (1 - cov_yz**2) + ) + + self.stat = stat + return stat + + def test( + self, + x, + y, + z, + reps=1000, + workers=1, + auto=True, + perm_blocks=None, + random_state=None, + ): + r""" + Calculates the partial correlation test statistic and p-value. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, 1)``, ``(n, 1)`` and + ``(n, 1)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. + reps : int, default: 1000 + The number of replications used to estimate the null distribution + when using the permutation test used to calculate the p-value. + workers : int, default: 1 + The number of cores to parallelize the p-value computation over. + Supply ``-1`` to use all cores available to the Process. + auto : bool, default: True + If True, the p-value is computed by t-distribution approximation. + Parameters ``reps`` and ``workers`` are irrelevant in this case. + Otherwise, :class:`hyppo.tools.perm_test` will be run. + perm_blocks : None or ndarray, default: None + Defines blocks of exchangeable samples during the permutation test. + If None, all samples can be permuted with one another. Requires `n` + rows. At each column, samples with matching column value are + recursively partitioned into blocks of samples. Within each final + block, samples are exchangeable. Blocks of samples from the same + partition are also exchangeable between one another. If a column + value is negative, that block is fixed and cannot be exchanged. + random_state : int, default: None + The random_state for permutation testing to be fixed for + reproducibility. + + Returns + ------- + stat : float + The computed partial correlation test statistic. + pvalue : float + The computed partial correlation test p-value. + """ + check_input = _CheckInputs(x, y, z, reps=reps, max_dims=1) + x, y, z = check_input() + + if auto: # run t-stat + stat = self.statistic(x, y, z) + + n = x.shape[0] + dof = n - 3 + tstat = stat * np.sqrt(dof) / np.sqrt(1 - stat**2) + pvalue = 1 - t.cdf(np.abs(tstat), dof) + + self.stat = stat + self.pvalue = pvalue + self.null_dist = None + else: # run permutation + stat, pvalue, null_dist = perm_test( + self.statistic, + x=x, + y=y, + z=z, + reps=reps, + workers=workers, + is_distsim=False, + perm_blocks=perm_blocks, + random_state=random_state, + ) + + self.stat = stat + self.pvalue = pvalue + self.null_dist = null_dist + + return ConditionalIndependenceTestOutput(stat, pvalue) diff --git a/hyppo/conditional/pdcorr.py b/hyppo/conditional/pdcorr.py new file mode 100644 index 000000000..c2ef29d7d --- /dev/null +++ b/hyppo/conditional/pdcorr.py @@ -0,0 +1,257 @@ +import numpy as np +from numba import jit + +from ..independence.dcorr import _center_distmat +from ..tools import compute_dist, perm_test +from ._utils import _CheckInputs +from .base import ConditionalIndependenceTest, ConditionalIndependenceTestOutput + + +class PartialDcorr(ConditionalIndependenceTest): + r""" + Partial Distance Covariance/Correlation (PDcov/PDcorr) test statistic and p-value. + + PDcorr is a measure of dependence between two paired random matrices + given a third random matrix of not necessarily equal dimensions :footcite:p:`szekelyPartialDistanceCorrelation2014a`. + + Parameters + ---------- + compute_distance : str, callable, or None, default: "euclidean" + A function that computes the distance among the samples within each + data matrix. + Valid strings for ``compute_distance`` are, as defined in + :func:`sklearn.metrics.pairwise_distances`, + + - From scikit-learn: [``"euclidean"``, ``"cityblock"``, ``"cosine"``, + ``"l1"``, ``"l2"``, ``"manhattan"``] See the documentation for + :mod:`scipy.spatial.distance` for details + on these metrics. + - From scipy.spatial.distance: [``"braycurtis"``, ``"canberra"``, + ``"chebyshev"``, ``"correlation"``, ``"dice"``, ``"hamming"``, + ``"jaccard"``, ``"kulsinski"``, ``"mahalanobis"``, ``"minkowski"``, + ``"rogerstanimoto"``, ``"russellrao"``, ``"seuclidean"``, + ``"sokalmichener"``, ``"sokalsneath"``, ``"sqeuclidean"``, + ``"yule"``] See the documentation for :mod:`scipy.spatial.distance` for + details on these metrics. + + Set to ``None`` or ``"precomputed"`` if ``x`` and ``y`` are already distance + matrices. To call a custom function, either create the distance matrix + before-hand or create a function of the form ``metric(x, **kwargs)`` + where ``x`` is the data matrix for which pairwise distances are + calculated and ``**kwargs`` are extra arguements to send to your custom + function. + use_cov : bool, + If `True`, then the statistic will compute the covariance rather than the + correlation. + **kwargs + Arbitrary keyword arguments for ``compute_distance``. + + Notes + ----- + The statistic can be derived as follows: + + Let :math:`x`, :math:`y`, and :math:`z` be :math:`(n, p)` samples of random + variables :math:`X`, :math:`Y` and :math:`Z`. Let :math:`D^x` be the :math:`n \times n` + distance matrix of :math:`x`, :math:`D^y` be the :math:`n \times n` be + the distance matrix of :math:`y`, and :math:`D^z` be the :math:`n \times n` distance + matrix of :math:`z`. Let :math:`C^x`, :math:`C^y`, and :math:`C^z` be the unbiased centered + distance matrices (see :class:`hyppo.independence.Dcorr` for more details). The + partial distance covariance is defined as + + .. math:: + \mathrm{PDcov}_n (x, y; z) &= \frac{1}{n(n-3)} \sum_{i\neq j}^n \left(P_{z^\perp}(x)\right)_{i,j} \left(P_{z^\perp}(y)\right)_{i,j} + where + .. math:: + P_{z^\perp}(x) &= C^x - \frac{(C^x\cdot C^z)}{ C^z \cdot C^z) C^z + is the orthogonal proejction of :math:`C^x` onto the subspace orthogonal to :math:`C^z`. + The partial distance correlation is defined as + + .. math:: + \mathrm{PDcorr}_n (x, y; z) &= \frac{P_{z^\perp}(x)\cdot P_{z^\perp}(y)}{\abs{P_{z^\perp}(x)}\abs{P_{z^\perp}(y)}} + + Equivalently, the partial distance correlation can be also defined as + + .. math:: + \mathrm{CDcorr}_n (x, y; z) &= \frac{R_{xy} - R_{xz} R_{yz}}{\sqrt{(1 - R_{xz}^2)(1 - R_{yz}^2)}} + where :math:`R_{xy}` is the unbiased distance correlation between :math:`x` and :math:`y`. + + References + ---------- + .. footbibliography:: + """ + + def __init__(self, compute_distance="euclidean", use_cov=True, **kwargs): + self.use_cov = use_cov + self.compute_distance = compute_distance + + # set is_distance to true if compute_distance is None + self.is_distance = False + if not compute_distance: + self.is_distance = True + + ConditionalIndependenceTest.__init__(self, **kwargs) + + def __repr__(self): + return "PartialDcorr" + + def statistic(self, x, y, z): + r""" + Helper function that calculates the PDcov/PDcorr test statistic. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x`` and ``y`` can be + distance matrices and ``z`` can be a similarity matrix where the + shapes must be ``(n, n)``. + + Returns + ------- + stat : float + The computed PDcov/PDcorr statistic. + """ + check_input = _CheckInputs(x, y, z) + x, y, z = check_input() + + if not self.is_distance: + distx, disty, distz = compute_dist( + x, + y, + z, + metric=self.compute_distance, + ) + else: + distx = x + disty = y + distz = z + + if self.use_cov: + stat = _pdcov(distx, disty, distz) + else: + stat = _pdcorr(distx, disty, distz) + + self.stat = stat + return stat + + def test( + self, + x, + y, + z, + reps=1000, + workers=1, + random_state=None, + ): + r""" + Calculates the PDcov/PDcorr test statistic and p-value. + + Parameters + ---------- + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x`` and ``y`` can be + distance matrices and ``z`` can be a similarity matrix where the + shapes must be ``(n, n)``. + reps : int, default: 1000 + The number of replications used to estimate the null distribution + when using the permutation test used to calculate the p-value. + workers : int, default: 1 + The number of cores to parallelize the p-value computation over. + Supply ``-1`` to use all cores available to the Process. + random_state : int, default: None + The random_state for permutation testing to be fixed for + reproducibility. + + Returns + ------- + stat : float + The computed PDcov/PDcorr statistic. + pvalue : float + The computed PDcov/PDcorr p-value. + """ + check_input = _CheckInputs(x, y, z, reps=reps) + x, y, z = check_input() + + if not self.is_distance: + x, y, z = compute_dist(x, y, z, metric=self.compute_distance, **self.kwargs) + self.is_distance = True + + stat, pvalue, null_dist = perm_test( + self.statistic, + x=x, + y=y, + z=z, + reps=reps, + workers=workers, + is_distsim=self.is_distance, + random_state=random_state, + ) + self.stat = stat + self.pvalue = pvalue + self.null_dist = null_dist + + return ConditionalIndependenceTestOutput(stat, pvalue) + + +@jit(nopython=True, cache=True) +def _pdcov(distx, disty, distz): + """Calculate the PDcov test statistic""" + N = distx.shape[0] + denom = N * (N - 3) + + distx = _center_distmat(distx, bias=False) + disty = _center_distmat(disty, bias=False) + distz = _center_distmat(distz, bias=False) + + cov_xz = np.sum(distx * distz) + cov_yz = np.sum(disty * distz) + var_z = np.sum(distz * distz) + + if var_z == 0: + stat = 0 + else: + proj_xz = distx - (cov_xz / var_z) * distz + proj_yz = disty - (cov_yz / var_z) * distz + stat = np.sum(proj_xz * proj_yz) / denom + + return stat + + +@jit(nopython=True, cache=True) +def _pdcorr(distx, disty, distz): + """Calculate the PDcorr test statistic""" + + distx = _center_distmat(distx, bias=False) + disty = _center_distmat(disty, bias=False) + distz = _center_distmat(distz, bias=False) + + cov_xy = np.sum(distx * disty) + cov_xz = np.sum(distx * distz) + cov_yz = np.sum(disty * distz) + var_x = np.sum(distx * distx) + var_y = np.sum(disty * disty) + var_z = np.sum(distz * distz) + + if var_x <= 0 or var_y <= 0: + rxy = 0 + else: + rxy = cov_xy / np.real(np.sqrt(var_x * var_y)) + + if var_x <= 0 or var_z <= 0: + rxz = 0 + else: + rxz = cov_xz / np.real(np.sqrt(var_x * var_z)) + + if var_y <= 0 or var_z <= 0: + ryz = 0 + else: + ryz = cov_yz / np.real(np.sqrt(var_y * var_z)) + + stat = (rxy - rxz * ryz) / np.sqrt((1 - rxz**2) * (1 - ryz**2)) + + return stat diff --git a/hyppo/conditional/tests/test_cdcorr.py b/hyppo/conditional/tests/test_cdcorr.py new file mode 100644 index 000000000..23982f775 --- /dev/null +++ b/hyppo/conditional/tests/test_cdcorr.py @@ -0,0 +1,42 @@ +import numpy as np +import pytest +from numpy.testing import assert_almost_equal + +from ...tools import indep_normal, power +from .. import ConditionalDcorr + + +class TestCDcorrStat: + @pytest.mark.parametrize("n", [100, 200]) + @pytest.mark.parametrize("obs_stat", [0.0]) + def test_linear_oned(self, n, obs_stat): + np.random.seed(123456789) + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = ConditionalDcorr().test(x, y, z) + stat2 = ConditionalDcorr().statistic(x, y, z) + + assert_almost_equal(stat1, obs_stat, decimal=2) + assert_almost_equal(stat2, obs_stat, decimal=2) + + @pytest.mark.parametrize("n", [100, 200]) + def test_rep(self, n): + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = ConditionalDcorr().test(x, y, z, random_state=2) + stat2, pvalue2 = ConditionalDcorr().test(x, y, z, random_state=2) + + assert stat1 == stat2 + assert pvalue1 == pvalue2 + + +class TestCDcorrTypeIError: + def test_oned(self): + np.random.seed(123456789) + est_power = power( + "ConditionalDcorr", + sim_type="condi", + sim="independent_normal", + n=100, + p=1, + alpha=0.05, + ) + assert_almost_equal(0.05, est_power, decimal=2) diff --git a/hyppo/conditional/tests/test_pcorr.py b/hyppo/conditional/tests/test_pcorr.py new file mode 100644 index 000000000..abb01ea8b --- /dev/null +++ b/hyppo/conditional/tests/test_pcorr.py @@ -0,0 +1,42 @@ +import numpy as np +import pytest +from numpy.testing import assert_almost_equal + +from ...tools import indep_normal, power +from .. import PartialCorr + + +class TestPcorrStat: + @pytest.mark.parametrize("n", [100, 200]) + @pytest.mark.parametrize("obs_stat", [0.0]) + def test_linear_oned(self, n, obs_stat): + np.random.seed(3) + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = PartialCorr().test(x, y, z) + stat2 = PartialCorr().statistic(x, y, z) + + assert_almost_equal(stat1, obs_stat, decimal=1) + assert_almost_equal(stat2, obs_stat, decimal=1) + + @pytest.mark.parametrize("n", [100, 200]) + def test_rep(self, n): + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = PartialCorr().test(x, y, z, random_state=2) + stat2, pvalue2 = PartialCorr().test(x, y, z, random_state=2) + + assert stat1 == stat2 + assert pvalue1 == pvalue2 + + +class TestPcorrTypeIError: + def test_oned(self): + np.random.seed(1) + est_power = power( + "PartialCorr", + sim_type="condi", + sim="independent_normal", + n=100, + p=1, + alpha=0.05, + ) + assert_almost_equal(0.05, est_power, decimal=2) diff --git a/hyppo/conditional/tests/test_pdcorr.py b/hyppo/conditional/tests/test_pdcorr.py new file mode 100644 index 000000000..79e35dffa --- /dev/null +++ b/hyppo/conditional/tests/test_pdcorr.py @@ -0,0 +1,42 @@ +import numpy as np +import pytest +from numpy.testing import assert_almost_equal + +from ...tools import indep_normal, power +from .. import PartialDcorr + + +class TestPDcorrStat: + @pytest.mark.parametrize("n", [100, 200]) + @pytest.mark.parametrize("obs_stat", [0.0]) + def test_linear_oned(self, n, obs_stat): + np.random.seed(123456789) + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = PartialDcorr().test(x, y, z) + stat2 = PartialDcorr().statistic(x, y, z) + + assert_almost_equal(stat1, obs_stat, decimal=2) + assert_almost_equal(stat2, obs_stat, decimal=2) + + @pytest.mark.parametrize("n", [100, 200]) + def test_rep(self, n): + x, y, z = indep_normal(n, 1) + stat1, pvalue1 = PartialDcorr().test(x, y, z, random_state=2) + stat2, pvalue2 = PartialDcorr().test(x, y, z, random_state=2) + + assert stat1 == stat2 + assert pvalue1 == pvalue2 + + +class TestPDcorrTypeIError: + def test_oned(self): + np.random.seed(123456789) + est_power = power( + "PartialDcorr", + sim_type="condi", + sim="independent_normal", + n=100, + p=1, + alpha=0.05, + ) + assert_almost_equal(est_power, 0.05, decimal=2) diff --git a/hyppo/conditional/tests/test_utils.py b/hyppo/conditional/tests/test_utils.py new file mode 100644 index 000000000..0302585c9 --- /dev/null +++ b/hyppo/conditional/tests/test_utils.py @@ -0,0 +1,59 @@ +from itertools import permutations + +import numpy as np +import pytest +from numpy.testing import assert_raises + +from .._utils import _CheckInputs + + +class TestErrorWarn: + """Tests errors and warnings.""" + + def test_error_notndarray(self): + # raises error if x or y is not a ndarray + x = np.arange(20) + y = np.arange(20) + z = [5] * 20 + + for data in permutations((x, y, z)): + assert_raises(TypeError, _CheckInputs(*data)) + + def test_error_shape(self): + # raises error if number of samples different (n) + x = np.arange(100).reshape(25, 4) + y = x.reshape(10, 10) + z = x.reshape(20, 5) + assert_raises(ValueError, _CheckInputs(x, y, z)) + + def test_error_lowsamples(self): + # raises error if samples are low (< 3) + x = np.arange(3) + y = np.arange(3) + z = np.arange(3) + assert_raises(ValueError, _CheckInputs(x, y, z)) + + def test_error_nans(self): + # raises error if inputs contain NaNs + x = np.arange(20, dtype=float) + x[0] = np.nan + assert_raises(ValueError, _CheckInputs(x, x, x)) + + y = np.arange(20) + assert_raises(ValueError, _CheckInputs(x, y, y)) + + @pytest.mark.parametrize( + "reps", [-1, "1"] # reps is negative # reps is not integer + ) + def test_error_reps(self, reps): + # raises error if reps is negative + x = np.arange(20) + assert_raises(ValueError, _CheckInputs(x, x, x, reps=reps)) + + def test_max_sims(self): + # raises error if data contains more dimentions than allowed + x = np.arange(20).reshape(10, 2) + y = np.arange(10).reshape(10, 1) + z = y + + assert_raises(ValueError, _CheckInputs(x, y, z, max_dims=1)) diff --git a/hyppo/independence/friedman_rafsky.py b/hyppo/independence/friedman_rafsky.py index a3f032d41..aa1fc582b 100644 --- a/hyppo/independence/friedman_rafsky.py +++ b/hyppo/independence/friedman_rafsky.py @@ -122,8 +122,8 @@ def test( self.statistic, x, y, - reps, - workers, + reps=reps, + workers=workers, is_distsim=False, random_state=random_state, ) diff --git a/hyppo/independence/kmerf.py b/hyppo/independence/kmerf.py index 1dd363604..33927eac3 100644 --- a/hyppo/independence/kmerf.py +++ b/hyppo/independence/kmerf.py @@ -186,6 +186,9 @@ def statistic(self, x, y): # get feature importances from gini-based random forest self.importances = self.clf.feature_importances_ + # get feature importances from gini-based random forest + self.importances = self.clf.feature_importances_ + return stat def test(self, x, y, reps=1000, workers=1, auto=True, random_state=None): diff --git a/hyppo/time_series/base.py b/hyppo/time_series/base.py index 7aa2329ce..4bc47ec2a 100644 --- a/hyppo/time_series/base.py +++ b/hyppo/time_series/base.py @@ -6,6 +6,7 @@ from sklearn.utils import check_random_state from ..tools import compute_dist +from sklearn.utils import check_random_state class TimeSeriesTestOutput(NamedTuple): diff --git a/hyppo/tools/__init__.py b/hyppo/tools/__init__.py index f0146c3b1..60bbba7b6 100755 --- a/hyppo/tools/__init__.py +++ b/hyppo/tools/__init__.py @@ -1,4 +1,5 @@ from .common import * +from .conditional_indep_sim import * from .indep_sim import * from .ksample_sim import * from .power import * diff --git a/hyppo/tools/common.py b/hyppo/tools/common.py index 8e47ed591..82f968082 100644 --- a/hyppo/tools/common.py +++ b/hyppo/tools/common.py @@ -47,6 +47,16 @@ def check_ndarray_xy(x, y): raise TypeError("x and y must be ndarrays") +def check_ndarray_xyz(x, y, z): + """Check if x or y is an ndarray of float""" + if ( + not isinstance(x, np.ndarray) + or not isinstance(y, np.ndarray) + or not isinstance(z, np.ndarray) + ): + raise TypeError("x, y, and z must be ndarrays") + + def convert_xy_float64(x, y): """Convert x or y to np.float64 (if not already done)""" # convert x and y to floats @@ -56,6 +66,16 @@ def convert_xy_float64(x, y): return x, y +def convert_xyz_float64(x, y, z): + """Convert x or y or z to np.float64 (if not already done)""" + # convert x and y to floats + x = np.asarray(x).astype(np.float64) + y = np.asarray(y).astype(np.float64) + z = np.asarray(z).astype(np.float64) + + return x, y, z + + def check_reps(reps): """Check if reps is valid""" # check if reps is an integer > than 0 @@ -72,29 +92,35 @@ def check_reps(reps): warnings.warn(msg, RuntimeWarning) -def _check_distmat(x, y): - """Check if x and y are distance matrices.""" - if ( - not np.allclose(x, x.T) - or not np.allclose(y, y.T) - or not np.all((x.diagonal() == 0)) - or not np.all((y.diagonal() == 0)) - ): - raise ValueError( - "x and y must be distance matrices, {is_sym} symmetric and " - "{zero_diag} zeros along the diagonal".format( - is_sym="x is not" - if not np.array_equal(x, x.T) - else "y is not" - if not np.array_equal(y, y.T) - else "both are", - zero_diag="x doesn't have" - if not np.all((x.diagonal() == 0)) - else "y doesn't have" - if not np.all((y.diagonal() == 0)) - else "both have", - ) - ) +def _check_distmat(x, y, z=None): + """Check if x, y, and z are distance matrices.""" + if z is None: + data = (x, y) + labs = np.array(["x", "y"]) + else: + data = (x, y, z) + labs = np.array(["x", "y", "z"]) + + check_sym = np.array([not np.allclose(arr, arr.T) for arr in data]) + check_diag = np.array([not np.all((arr.diagonal() == 0)) for arr in data]) + + if np.any(check_sym) or np.any(check_diag): + inputs = "x and y" if z is None else "x, y and z" + if np.any(check_sym): + names = ", ".join(labs[check_sym]) + verb = "is" if check_sym.sum() == 1 else "are" + sym_msg = f"{names} {verb} not symmetric." + else: + sym_msg = "" + if np.any(check_diag): + names = ", ".join(labs[check_diag]) + verb = "does not" if check_diag.sum() == 1 else "do not" + diag_msg = f"{names} {verb} have zeros along the diagonal." + else: + diag_msg = "" + + msg = f"{inputs} must be distance matrices. {sym_msg} {diag_msg}" + raise ValueError(msg) def _check_kernmat(x, y): @@ -281,18 +307,19 @@ def multi_compute_kern(*args, metric="gaussian", workers=1, **kwargs): return sim_matrices -def compute_dist(x, y, metric="euclidean", workers=1, **kwargs): +def compute_dist(x, y, z=None, metric="euclidean", workers=1, **kwargs): """ Distance matrices for the inputs. Parameters ---------- - x,y : ndarray of float - Input data matrices. ``x`` and ``y`` must have the same number of - samples. That is, the shapes must be ``(n, p)`` and ``(n, q)`` where - `n` is the number of samples and `p` and `q` are the number of - dimensions. Alternatively, ``x`` and ``y`` can be distance matrices, - where the shapes must both be ``(n, n)``. + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number + of samples. That is, the shapes must be ``(n, p)``, ``(n, q)`` and + ``(n, r)`` where `n` is the number of samples and `p`, `q`, and `r` + are the number of dimensions. Alternatively, ``x``, ``y`` and ``z`` + can be distance matrices where the shapes must be ``(n, n)``. ``z`` + is an optional ndarray. metric : str, callable, or None, default: "euclidean" A function that computes the distance among the samples within each data matrix. @@ -311,8 +338,8 @@ def compute_dist(x, y, metric="euclidean", workers=1, **kwargs): ``"yule"``] See the documentation for :mod:`scipy.spatial.distance` for details on these metrics. - Set to ``None`` or ``"precomputed"`` if ``x`` and ``y`` are already distance - matrices. To call a custom function, either create the distance matrix + Set to ``None`` or ``"precomputed"`` if ``x``, ``y``, and ``z`` are already + distance matrices. To call a custom function, either create the distance matrix before-hand or create a function of the form ``metric(x, **kwargs)`` where ``x`` is the data matrix for which pairwise distances are calculated and ``**kwargs`` are extra arguements to send to your custom @@ -327,21 +354,34 @@ def compute_dist(x, y, metric="euclidean", workers=1, **kwargs): Returns ------- - distx, disty : ndarray of float - Distance matrices based on the metric provided by the user. + distx, disty, distz : ndarray of float + Distance matrices based on the metric provided by the user. ``distz`` is + only returned if ``z`` is provided. """ if not metric: metric = "precomputed" if callable(metric): distx = metric(x, **kwargs) disty = metric(y, **kwargs) + distz = None if z is None else metric(z, **kwargs) _check_distmat( - distx, disty + distx, + disty, + distz, ) # verify whether matrix is correct, built into sklearn func else: distx = pairwise_distances(x, metric=metric, n_jobs=workers, **kwargs) disty = pairwise_distances(y, metric=metric, n_jobs=workers, **kwargs) - return distx, disty + distz = ( + None + if z is None + else pairwise_distances(z, metric=metric, n_jobs=workers, **kwargs) + ) + + if z is None: + return distx, disty + else: + return distx, disty, distz def check_perm_blocks(perm_blocks): @@ -470,8 +510,7 @@ def __init__(self, y, perm_blocks=None): self.perm_tree = _PermTree(perm_blocks) def __call__(self, rng=None): - if rng is None: - rng = np.random + rng = check_random_state(rng) if self.perm_tree is None: order = rng.permutation(self.n) else: @@ -481,20 +520,25 @@ def __call__(self, rng=None): # p-value computation -def _perm_stat(calc_stat, x, y, is_distsim=True, permuter=None, random_state=None): +def _perm_stat( + calc_stat, x, y, z=None, is_distsim=True, permuter=None, random_state=None +): """Permute the test statistic""" rng = check_random_state(random_state) if permuter is None: order = rng.permutation(y.shape[0]) else: - order = permuter(rng) + order = permuter(rng=rng) if is_distsim: permy = y[order][:, order] else: permy = y[order] - perm_stat = calc_stat(x, permy) + if z is not None: + perm_stat = calc_stat(x, permy, z) + else: + perm_stat = calc_stat(x, permy) return perm_stat @@ -503,11 +547,13 @@ def perm_test( calc_stat, x, y, + z=None, reps=1000, workers=1, is_distsim=True, perm_blocks=None, random_state=None, + permuter=None, ): """ Permutation test for the p-value of a nonparametric test. @@ -521,13 +567,14 @@ def perm_test( ---------- calc_stat : callable The method used to calculate the test statistic (must use hyppo API). - x,y : ndarray of float - Input data matrices. ``x`` and ``y`` must have the same number of - samples. That is, the shapes must be ``(n, p)`` and ``(n, q)`` where - `n` is the number of samples and `p` and `q` are the number of + x,y,z : ndarray of float + Input data matrices. ``x``, ``y`` and ``z`` must have the same number of + samples. That is, the shapes must be ``(n, p)``, ``(n, q)``, ``(n, r)``where + `n` is the number of samples and `p`, `q` , and `r` are the number of dimensions. Alternatively, ``x`` and ``y`` can be distance or similarity - matrices, - where the shapes must both be ``(n, n)``. + matrices, and ``z`` must be a similarity matrix where the shapes + must be ``(n, n)``. ``z`` is an optional matrix only used for conditional + independence testing. reps : int, default: 1000 The number of replications used to estimate the null distribution when using the permutation test used to calculate the p-value. @@ -545,6 +592,10 @@ def perm_test( test, samples within the same final leaf node are exchangeable and blocks of samples with a common parent node are exchangeable. If a column value is negative, the resulting block is unexchangeable. + permuter : callable, default: None + Defines a custom permutation function. If None, the default permutation + test is used. + Returns ------- stat : float @@ -554,8 +605,12 @@ def perm_test( null_dist : list of float The approximated null distribution of shape ``(reps,)``. """ + # calculate observed test statistic - stat = calc_stat(x, y) + if z is None: + stat = calc_stat(x, y) + else: + stat = calc_stat(x, y, z) # make RandomState seeded array if random_state is not None: @@ -567,12 +622,13 @@ def perm_test( random_state = np.random.randint(np.iinfo(np.int32).max, size=reps) # calculate null distribution - permuter = _PermGroups(y, perm_blocks) + if not callable(permuter): + permuter = _PermGroups(y, perm_blocks) null_dist = np.array( Parallel(n_jobs=workers)( [ - delayed(_perm_stat)(calc_stat, x, y, is_distsim, permuter, rng) + delayed(_perm_stat)(calc_stat, x, y, z, is_distsim, permuter, rng) for rng in random_state ] ) diff --git a/hyppo/tools/conditional_indep_sim.py b/hyppo/tools/conditional_indep_sim.py new file mode 100644 index 000000000..daaa07a76 --- /dev/null +++ b/hyppo/tools/conditional_indep_sim.py @@ -0,0 +1,655 @@ +import numpy as np +from sklearn.utils import check_random_state + +from .indep_sim import _CheckInputs + + +def indep_normal(n, p=1, random_state=None): + r""" + Independent standard normal distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + + .. math:: + X, Y, Z &\sim N(0, 1) + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + Ignored. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + mean = np.array([0, 0, 0]) + cov = np.eye(3) + + x, y, z = rng.multivariate_normal(mean, cov, size=n).T + + return x, y, z + + +def indep_lognormal(n, p=1, random_state=None): + r""" + Independent lognormal and normal distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + .. math:: + X &\sim \text{log} N(0, 1)\\ + Y, Z &\sim N(0, 1) + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + x, y, z = indep_normal(n=n, p=p, random_state=random_state) + x = np.exp(x) + + return x, y, z + + +def indep_binomial(n, p=1, random_state=None): + r""" + Independent binomial distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^p`: + .. math:: + X, Y, Z_i &\sim \text{Binom}(10, 0.5) \\ + Z &= (Z_1, Z_2, \ldots, Z_p) + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + The number of dimensions for conditioning variable ``Z``. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x`` and ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + size = (n, 1) + + z = rng.binomial(10, 0.5, size=(n, p)) + x = rng.binomial(10, 0.5, size=size) + y = rng.binomial(10, 0.5, size=size) + + return x, y, z + + +def cond_indep_normal(n, p=1, random_state=None): + r""" + Conditionally independent normal distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + .. math:: + \mu &= (0, 0, 0)\\ + \Sigma &= \begin{bmatrix} + 1 & 0.36 & 0.6 \\ + 0.36 & 1 & 0.6 \\ + 0.6 & 0.6 & 1 + \end{bmatrix}\\ + (X, Y, Z) &\sim MVN(\mu, \Sigma) + + The conditional covariance matrix is given by: + .. math:: + \Sigma(X, Y | Z) &= \begin{bmatrix} + 0.64 & 0\\ + 0 & 0.64 + \end{bmatrix} + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + Ignored. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + mean = np.array([0, 0, 0]) + cov = np.array([[1, 0.36, 0.6], [0.36, 1, 0.6], [0.6, 0.6, 1]]) + + x, y, z = rng.multivariate_normal(mean, cov, size=n).T + + return x, y, z + + +def cond_indep_lognormal(n, p=1, random_state=None): + r""" + Conditionally independent lognormal and normal distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + Ignored. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + mean = np.array([0, 0, 0]) + cov = np.array([[1, 0.36, 0.6], [0.36, 1, 0.6], [0.6, 0.6, 1]]) + + x, y, z = rng.multivariate_normal(mean, cov, size=n).T + x = np.exp(x) + + return x, y, z + + +def cond_indep_normal_nonlinear(n, p=1, random_state=None): + r""" + Conditionally independent normal distributions. Example 3 from :footcite:p:`wang2015conditional`. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + .. math:: + X_1, Y_1, Z &\sim N(0, 1) \\ + Z_1 &= 0.5 \left( \frac{Z^3}{7} + \frac{Z}{2} \right) \\ + Z_2 &= \frac{Z^3}{2} + \frac{Z}{3} \\ + X_2 &= Z_1 + \tanh(X_1) \\ + X &= X_2 + \frac{X_2^3}{3}\\ + Y_2 &= Z_2 + Y_1 \\ + Y &= Y_2 + \frac{Y_2^3}{3} + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + size = (n, 1) + z = rng.normal(size=size) + z1 = 0.5 * (np.power(z, 3) / 7 + z / 2) + z2 = (np.power(z, 3) / 2 + z) / 3 + + x1 = rng.normal(size=size) + x2 = z1 + np.tanh(x1) + x = x2 + np.power(x2, 3) / 3 + + y1 = rng.normal(size=size) + y2 = z2 + y1 + y = y2 + np.tanh(y2 / 3) + + return x, y, z + + +def cond_indep_binomial(n, p=1, random_state=None): + r""" + Conditionally independent binomial distributions. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}^p`: + + .. math:: + X_1, Y_1, Z_i &\sim \text{Binom}(10, 0.5) \\ + Z &= (Z_1, Z_2, \ldots, Z_p) \\ + X &= X_1 + Z_1 + Z_2 + \cdots + Z_p \\ + Y &= Y_1 + Z_1 + Z_2 + \cdots + Z_p + + Examples 2 and 4 from :footcite:p:`wang2015conditional`. + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + The number of dimensions for conditioning variable ``Z``. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + size = (n, 1) + + z = rng.binomial(10, 0.5, size=(n, p)) + x = rng.binomial(10, 0.5, size=size) + z.sum(axis=1, keepdims=True) + y = rng.binomial(10, 0.5, size=size) + z.sum(axis=1, keepdims=True) + + return x, y, z + + +def correlated_binomial(n, p=1, random_state=None): + r""" + Conditionally dependent binomial distributions. Examples 6 from:footcite:p:`wang2015conditional`. + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + + .. math:: + X_1, Z &\sim \text{Binom}(10, 0.5) \\ + X &= X_1 + Z \\ + Y &= (X_1 - 5)^4 + Z + + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + The number of dimensions for conditioning variable ``Z``. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + x1, z = rng.binomial(10, 0.5, size=(2, n)) + + x = x1 + z + y = (x1 - 5) ** 4 + z + + return x, y, z + + +def correlated_normal(n, p=1, random_state=None): + r""" + Conditionally dependent normal distributions. Example 5 from :footcite:p:`wang2015conditional` + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + .. math:: + \mu &= (0, 0, 0)\\ + \Sigma &= \begin{bmatrix} + 1 & 0.7 & 0.6\\ + 0.7 & 1 & 0.6\\ + 0.6 & 0.6 & 1 + \end{bmatrix}\\ + (X, Y, Z) &\sim MVN(\mu, \Sigma) + + The conditional covariance matrix is given by: + .. math:: + \Sigma(X, Y | Z) &= \begin{bmatrix} + 0.64 & 0.34\\ + 0.34 & 0.64 + \end{bmatrix} + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + mean = np.array([0, 0, 0]) + cov = np.array([[1, 0.7, 0.6], [0.7, 1, 0.6], [0.6, 0.6, 1]]) + + x, y, z = rng.multivariate_normal(mean, cov, size=n).T + + return x, y, z + + +def correlated_normal_nonliear(n, p=1, random_state=None): + r""" + Conditionally dependent normal distributions with nonlinear dependence. + Example 7 from :footcite:p:`wang2015conditional` + + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + .. math:: + X_1, Y_1, Z, \epsilon &\sim N(0, 1) \\ + Z_1 &= 0.5(Z^3/7 + Z/2) \\ + Z_2 &= (Z^3/2 + Z)/3 \\ + X_2 &= Z_1 + \tanh(X_1) \\ + X_3 &= X_2 + X_2^3 / 3 \\ + Y_2 &= Z_2 + Y_1\\ + Y_3 &= Y_2 + \tanh(Y_2 / 3) \\ + X &= X_3 + \cosh\epsilon \\ + Y &= Y_3 + \cosh\epsilon^2 + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + x1, y1, z, eps = rng.normal(size=(n, 4)).T + + z1 = 0.5 * (z**3 / 7 + z / 2) + z2 = (z**3 / 2 + z) / 3 + + x2 = z1 + np.tanh(x1) + x3 = x2 + x2**3 / 3 + + y2 = z2 + y1 + y3 = y2 + np.tanh(y2 / 3) + + x = x3 + np.cosh(eps) + y = y3 + np.cosh(eps**2) + + return x, y, z + + +def correlated_lognormal(n, p=1, random_state=None): + r""" + Example 5 from :footcite:p:`szekelyPartialDistanceCorrelation2014a` + :math:`(X, Y, Z) \in \mathbb{R} \times \mathbb{R} \times \mathbb{R}`: + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + x, y, z = correlated_normal(n, p, random_state) + x = np.exp(x) + + return x, y, z + + +def correlated_t_linear(n, p=4, random_state=None): + r""" + Conditionally dependent t-distributed data with linear dependence. + Example 9 from :footcite:p:`wang2015conditional` + + :math:`(X, Y, Z) \in \mathbb{R}^p \times \mathbb{R} \times \mathbb{R}`: + .. math:: + Z_1, Z_2, \cdots, Z_p, Z_{p+1}, Z_{p+2} &\sim t(1)\\ + X &= (Z_1, Z_2, Z_{p-1}, Z_p + Z_{p+1})\\ + Y &= Z_1 + Z_2 + \cdots + Z_p + Z_{p+1} + Z_{p+2}\\ + Z &= Z_{p+1} + + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + p : int + The number of dimentions for variable ``x``. + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + z_i = rng.standard_t(df=1, size=(n, p + 2)) + + x = z_i[:, :p] + x[:, -1] += z_i[:, p] + y = z_i.sum(axis=1) + z = z_i[:, p] + + return x, y, z + + +def correlated_t_quadratic(n, p=10, random_state=None): + r""" + Conditionally dependent t-distributed data with quadratic dependence. + Example 10 from :footcite:p:`wang2015conditional` + + :math:`(X, Y, Z) \in \mathbb{R}^10 \times \mathbb{R}^2 \times \mathbb{R}`: + .. math:: + Z_1, Z_2, \cdots, Z_{13} &\sim t(1)\\ + X_i &= Z_i, i = 1, 2, \cdots, 9\\ + X_{10} &= Z_{10} + Z_{11}\\ + Y_1 &= Z_1 Z_2 + Z_3 Z_4 + Z_5 Z_{11} + Z_{12}\\ + Y_2 &= Z_6 Z_7 + Z_8 Z_9 + Z_{10} Z_{11} + Z_{13}\\ + Z &= Z_{11} + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + z = rng.standard_t(df=1, size=(13, n)) + + x = z[:11] + x[9] += z[10] + y1 = z[0] * z[1] + z[2] * z[3] + z[4] * z[10] + z[11] + y2 = z[5] * z[6] + z[7] * z[8] + z[9] * z[10] + z[12] + y = np.array([y1, y2]) + z = z[10] + + return x.T, y.T, z + + +def correlated_t_nonlinear(n, p=4, random_state=None): + r""" + Conditionally dependent t-distributed data with nonlinear dependence. + Example 11 from :footcite:p:`wang2015conditional` + + :math:`(X, Y, Z) \in \mathbb{R}^4 \times \mathbb{R}^2 \times \mathbb{R}^2`: + .. math:: + Z_1, \ldots, Z_4 &\sim t(2)\\ + Y_1 &= \sin(Z_1) + \cos(Z_2) + Z_3^2 + Z_4^2\\ + Y_2 &= Z_1^2 + Z_2^2 + Z_3 + Z_4\\ + X &= (Z_1, Z_2, Z_3, Z_4)\\ + Y &= (Y_1, Y_2)\\ + Z &= (Z_1, Z_2) + + Parameters + ---------- + n : int + The number of samples desired by the simulation (>= 5). + + Returns + ------- + x,y,z : ndarray of float + Simulated data matrices. ``x``, ``y``, and ``z``. + + References + ---------- + .. footbibliography:: + """ + check_in = _CheckInputs(n, p=p) + check_in() + + rng = check_random_state(random_state) + + z = rng.standard_t(df=2, size=(4, n)) + y1 = np.sin(z[0]) + np.cos(z[1]) + z[2] ** 2 + z[3] ** 2 + y2 = z[0] ** 2 + z[1] ** 2 + z[2] + z[3] + y = np.array([y1, y2]) + + return z.T, y.T, z[:2].T + + +COND_SIMULATIONS = { + "independent_normal": indep_normal, + "independent_lognormal": indep_lognormal, + "independent_binomial": indep_binomial, + "cond_independent_normal": cond_indep_normal, + "cond_independent_lognormal": cond_indep_lognormal, + "cond_independent_binomial": cond_indep_binomial, + "cond_independent_normal_nonlinear": cond_indep_normal_nonlinear, + "correlated_normal": correlated_normal, + "correlated_lognormal": correlated_lognormal, + "correlated_binomial": correlated_binomial, + "correlated_normal_nonliear": correlated_normal_nonliear, + "correlated_t_linear": correlated_t_linear, + "correlated_t_nonlinear": correlated_t_nonlinear, + "correlated_t_quadratic": correlated_t_quadratic, +} + + +def condi_indep_sim(sim, n, p, random_state=None, **kwargs): + r""" + Conditional independence simulation generator. + + Takes a simulation and the required parameters, and outputs the simulated + data matrices. + + Parameters + ---------- + sim : str + The name of the simulation (from the :mod:`hyppo.tools` module) that is to be + rotated. + n : int + The number of samples desired by the simulation (>= 5). + p : int + The number of dimensions desired by the simulation (>= 1). + **kwargs + Additional keyword arguements for the desired simulation. + + Returns + ------- + x,y, z : ndarray of float + Simulated data matrices. + """ + if sim not in COND_SIMULATIONS.keys(): + raise ValueError( + "sim_name must be one of the following: {}".format( + list(COND_SIMULATIONS.keys()) + ) + ) + else: + sim = COND_SIMULATIONS[sim] + + x, y, z = sim(n, p, random_state, **kwargs) + + if x.ndim == 1: + x = x[:, np.newaxis] + if y.ndim == 1: + y = y[:, np.newaxis] + if z.ndim == 1: + z = z[:, np.newaxis] + + return x, y, z diff --git a/hyppo/tools/power.py b/hyppo/tools/power.py index ea5c688b0..393f65ff1 100644 --- a/hyppo/tools/power.py +++ b/hyppo/tools/power.py @@ -2,17 +2,20 @@ import numpy as np +from ..conditional import COND_INDEP_TESTS +from ..d_variate import MULTI_TESTS from ..independence import INDEP_TESTS from ..ksample import KSAMP_TESTS, KSample, k_sample_transform -from ..d_variate import MULTI_TESTS from .indep_sim import indep_sim from .ksample_sim import gaussian_3samp, rot_ksamp +from .conditional_indep_sim import condi_indep_sim _ALL_SIMS = { "indep": indep_sim, "multi": indep_sim, "ksamp": rot_ksamp, "gauss": gaussian_3samp, + "condi": condi_indep_sim, } _NONPERM_TESTS = { @@ -83,10 +86,32 @@ def _multi_perm_stat(test, sim_type, **kwargs): return obs_stat, perm_stat +def _conditional_perm_stat(test, sim_type, **kwargs): + """ + Generates null and alternate distributions for the conditional independence test. + """ + x, y, z = _sim_gen(sim_type=sim_type, **kwargs) + obs_stat = test.statistic(x, y, z) + + # Separate if block for conditional distance correlation since it requires + # special permutation + if str(test).lower() == "cdcorr": + probs = test._compute_kde(z) + idx = test._permuter(probs) + permy = y[idx] + perm_stat = test.statistic(x, permy, z) + else: + permy = np.random.permutation(y) + perm_stat = test.statistic(x, permy, z) + + return obs_stat, perm_stat + + _PERM_STATS = { "indep": _indep_perm_stat, "ksamp": _ksamp_perm_stat, "multi": _multi_perm_stat, + "condi": _conditional_perm_stat, } @@ -107,13 +132,13 @@ def power(test, sim_type, sim=None, n=100, alpha=0.05, reps=1000, auto=False, ** Parameters ---------- test : str or list - The name of the independence test (from the :mod:`hyppo.independence` module - or the :mod:`hyppo.d_variate` module) that is to be tested. If MaxMargin, - accepts list with first entry "MaxMargin" and second entry the name of - another independence test. + The name of the independence test (from the :mod:`hyppo.independence` module, + the :mod:`hyppo.d_variate` module, or the :mod:`hyppo.conditional` module) + that is to be tested. If MaxMargin, accepts list with first entry "MaxMargin" + and second entry the name of another independence test. For :class:`hyppo.ksample.KSample` put the name of the independence test. For other tests in :mod:`hyppo.ksample` just use the name of the class. - sim_type : "indep", "ksamp", "gauss", "multi" + sim_type : "indep", "ksamp", "gauss", "multi", "condi" Type of power method to calculate. Depends on the type of ``sim``. sim : str, default: None The name of the independence simulation (from the :mod:`hyppo.tools` module). @@ -162,13 +187,16 @@ def power(test, sim_type, sim=None, n=100, alpha=0.05, reps=1000, auto=False, ** test = KSAMP_TESTS[test_name]() elif test_name in MULTI_TESTS.keys(): test = MULTI_TESTS[test_name]() + elif test_name in COND_INDEP_TESTS.keys(): + test = COND_INDEP_TESTS[test_name]() else: raise ValueError( "Test {} not in {}".format( test_name, list(INDEP_TESTS.keys()) + list(KSAMP_TESTS.keys()) - + list(MULTI_TESTS.keys()), + + list(MULTI_TESTS.keys()) + + list(COND_INDEP_TESTS.keys()), ) ) @@ -178,6 +206,8 @@ def power(test, sim_type, sim=None, n=100, alpha=0.05, reps=1000, auto=False, ** perm_type = "ksamp" if sim_type == "multi": perm_type = "multi" + if sim_type == "condi": + perm_type = "condi" if sim_type != "gauss": kwargs["sim"] = sim diff --git a/hyppo/tools/tests/test_conditional_indep_sim.py b/hyppo/tools/tests/test_conditional_indep_sim.py new file mode 100644 index 000000000..cb1a9182b --- /dev/null +++ b/hyppo/tools/tests/test_conditional_indep_sim.py @@ -0,0 +1,66 @@ +import numpy as np +import pytest +from numpy.testing import assert_equal, assert_raises + +from .. import SIMULATIONS, indep_sim + + +class TestIndepSimShape: + @pytest.mark.parametrize("n", [100, 1000]) + @pytest.mark.parametrize("p", [1, 5]) + @pytest.mark.parametrize( + "sim", + SIMULATIONS.keys(), + ) + def test_shapes(self, n, p, sim): + np.random.seed(123456789) + x, y = SIMULATIONS[sim](n, p) + x1, y1 = indep_sim(sim, n, p) + nx, px = x.shape + ny, py = y.shape + nx1, px1 = x1.shape + ny1, py1 = y1.shape + + if sim in [ + "joint_normal", + "logarithmic", + "sin_four_pi", + "sin_sixteen_pi", + "two_parabolas", + "square", + "diamond", + "circle", + "ellipse", + "multiplicative_noise", + "multimodal_independence", + ]: + assert_equal(px, py) + assert_equal(px1, py1) + else: + assert_equal(px, py * p) + assert_equal(px1, py1 * p) + assert_equal(nx, ny) + assert_equal(nx1, ny1) + + +class TestIndepSimErrorWarn: + """Tests errors and warnings.""" + + def test_np_inctype(self): + assert_raises(ValueError, indep_sim, n="a", p=1, sim="linear") + assert_raises(ValueError, indep_sim, n=10, p=7.0, sim="linear") + + def test_low_n(self): + assert_raises(ValueError, indep_sim, n=3, p=1, sim="linear") + + def test_low_p(self): + assert_raises(ValueError, indep_sim, n=5, p=0, sim="linear") + + def test_extra_inctype(self): + assert_raises(ValueError, indep_sim, n=5, p=0, sim="linear", low="a") + + def test_wrong_sim(self): + assert_raises(ValueError, indep_sim, n=100, p=1, sim="abcd") + + def test_joint_lowcov(self): + assert_raises(ValueError, indep_sim, n=100, p=20, sim="joint_normal")