Repository: dougiesquire/xstatstests (language: python)

Issue:
Add normality tests



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index a1a8f44..bc144a7 100755
--- a/README.md
+++ b/README.md
@@ -26,6 +26,12 @@ Statistical tests on xarray objects. Currently includes the following statistica
 - `xstatstests.fligner` : Fligner-Killeen test for the variances of K independent samples
 - `xstatstests.median_test` : Mood test for the medians of K independent samples
 
+#### Normality tests
+- `xstatstests.skewtest` : D’Agostino test for normal sample skewness
+- `xstatstests.kurtosistest` : Anscombe test for normal sample kurtosis
+- `xstatstests.normaltest`: D’Agostino/Pearson test for normal sample skewness and kurtosis
+- `xstatstests.jarque_bera` : Jarque-Bera test for normal sample skewness and kurtosis
+
 #### Other tests
 - `xstatstests.mannwhitneyu` : Mann-Whitney U rank test on two independent samples
 - `xstatstests.ranksums` : Wilcoxon rank-sum statistic on two independent samples


Test patch:
diff --git a/xstatstests/__init__.py b/xstatstests/__init__.py
index be5ce7d..028f61f 100644
--- a/xstatstests/__init__.py
+++ b/xstatstests/__init__.py
@@ -1,6 +1,5 @@
-from .core import (
+from .scipy import (
     ks_2samp_1d,
-    ks_2samp_2d,
     anderson_ksamp,
     ttest_ind,
     ttest_rel,
@@ -17,6 +16,13 @@
     fligner,
     median_test,
     mood,
+    skewtest,
+    kurtosistest,
+    normaltest,
+    jarque_bera,
+)
+from .core import (
+    ks_2samp_2d,
 )
 
 from . import _version
diff --git a/xstatstests/core.py b/xstatstests/core.py
index f8e37f9..badbf08 100755
--- a/xstatstests/core.py
+++ b/xstatstests/core.py
@@ -1,876 +1,7 @@
-import inspect
-
 import numpy as np
 import xarray as xr
-import scipy.stats
-
-SAMPLE_DIM = "xsampletest_sample_dim"
-
-# Summary of how to wrap scipy.stats funcs
-# Enter outputs so that first output is the test statistic and the second
-# output is the p-value
-# n_args = -1 means no limit to the number of args
-scipy_function_info = {
-    "ks_2samp_1d": {
-        "name": "ks_2samp",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "anderson_ksamp": {
-        "name": "anderson_ksamp",
-        "stack_args": True,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": ["statistic", "significance_level"],
-    },
-    "ttest_ind": {
-        "name": "ttest_ind",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis"],
-        "outputs": [0, 1],
-    },
-    "ttest_rel": {
-        "name": "ttest_rel",
-        "stack_args": False,
-        "same_sample_sizes": True,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis"],
-        "outputs": [0, 1],
-    },
-    "cramervonmises_2samp": {
-        "name": "cramervonmises_2samp",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": ["statistic", "pvalue"],
-    },
-    "epps_singleton_2samp": {
-        "name": "epps_singleton_2samp",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "mannwhitneyu": {
-        "name": "mannwhitneyu",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis", "keepdims"],
-        "outputs": ["statistic", "pvalue"],
-    },
-    "ranksums": {
-        "name": "ranksums",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis", "keepdims"],
-        "outputs": [0, 1],
-    },
-    "kruskal": {
-        "name": "kruskal",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis", "keepdims"],
-        "outputs": [0, 1],
-    },
-    "friedmanchisquare": {
-        "name": "friedmanchisquare",
-        "stack_args": False,
-        "same_sample_sizes": True,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "brunnermunzel": {
-        "name": "brunnermunzel",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "ansari": {
-        "name": "ansari",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "bartlett": {
-        "name": "bartlett",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "levene": {
-        "name": "levene",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "fligner": {
-        "name": "fligner",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": True,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "median_test": {
-        "name": "median_test",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": [],
-        "outputs": [0, 1],
-    },
-    "mood": {
-        "name": "mood",
-        "stack_args": False,
-        "same_sample_sizes": False,
-        "remove_nans": False,
-        "disallowed_kwargs": ["axis"],
-        "outputs": [0, 1],
-    },
-}
-
-
-def _prep_data(*args, dim, nd):
-    """Prepare data for 2D tests"""
-    if isinstance(dim, str):
-        dim = [dim]
-
-    if any([(not isinstance(ds, xr.Dataset)) for ds in args]):
-        raise TypeError(
-            f"Input arrays must be xarray Datasets with {nd} variable(s) each"
-        )
-
-    args = xr.broadcast(*[ds.copy() for ds in args], exclude=dim)
-
-    if len(dim) == 1:
-        args = [ds.rename({dim[0]: SAMPLE_DIM}) for ds in args]
-    else:
-        args = [ds.stack({SAMPLE_DIM: dim}) for ds in args]
-
-    args = [ds.drop_vars({SAMPLE_DIM, *dim}, errors="ignore") for ds in args]
-    args = [ds.assign_coords({SAMPLE_DIM: range(ds.sizes[SAMPLE_DIM])}) for ds in args]
-
-    assert all(
-        [len(ds.data_vars) == nd for ds in args]
-    ), f"Input Datasets must have {nd} variables each"
-    assert all(
-        [list(args[0].data_vars) == list(ds.data_vars) for ds in args]
-    ), "Variables of all input Datasets must have the same name(s)"
-
-    # Need to rename sample dim otherwise apply_ufunc tries to align
-    # Expand into list of single variable Datasets
-    args_prepped = []
-    input_core_dims = []
-    for ind, ds in enumerate(args):
-        for var in ds.data_vars:
-            sample_dim = f"{SAMPLE_DIM}{ind+1}"
-            input_core_dims.append([sample_dim])
-            args_prepped.append(ds[var].rename({SAMPLE_DIM: sample_dim}))
-
-    return args_prepped, input_core_dims
-
-
-def _wrap_scipy(func, args, dim, kwargs):
-    """Generic xarray wrapper for subset of scipy stats functions"""
-
-    def _wrap_scipy_func(*args, scipy_func_info, scipy_kwargs):
-        """Parse scipy_function_info and apply scipy function"""
-        func = getattr(scipy.stats, scipy_func_info["name"])
-        getter = scipy_func_info["outputs"]
-
-        if scipy_func_info["remove_nans"]:
-            args = [arg[~np.isnan(arg)] for arg in args]
-        if scipy_func_info["stack_args"]:
-            outputs = func(args, **scipy_kwargs)
-        else:
-            outputs = func(*args, **scipy_kwargs)
-
-        return tuple(
-            [getattr(outputs, g) if isinstance(g, str) else outputs[g] for g in getter]
-        )
-
-    info = scipy_function_info[func]
-
-    args, input_core_dims = _prep_data(*args, dim=dim, nd=1)
-
-    sample_sizes = [ds.sizes[dim[0]] for ds, dim in zip(args, input_core_dims)]
-    same_sample_sizes = all([sample_sizes[0] == samp for samp in sample_sizes[1:]])
-    if info["same_sample_sizes"] & (not same_sample_sizes):
-        raise ValueError(
-            f"`{func}` requires that the sample size is the same for all input arrays"
-        )
-
-    scipy_kwargs = kwargs.copy()
-
-    # Simply vectorize if function does not allow axis argument
-    vectorize = True
-    for kw in info["disallowed_kwargs"]:
-        if kw in scipy_kwargs.keys():
-            raise ValueError(f"`{kw}` kwarg is disallowed by xstatstests")
-        if kw == "axis":
-            vectorize = False
-            scipy_kwargs["axis"] = -1
-
-    output_core_dims = [[]] * 2
-    output_dtypes = ["float32"] * 2
-    kwargs = dict(scipy_func_info=info, scipy_kwargs=scipy_kwargs)
-    statistic, pvalue = xr.apply_ufunc(
-        _wrap_scipy_func,
-        *args,
-        kwargs=kwargs,
-        input_core_dims=input_core_dims,
-        output_core_dims=output_core_dims,
-        output_dtypes=output_dtypes,
-        vectorize=vectorize,
-        dask="parallelized",
-    )
-
-    return xr.merge([statistic.rename("statistic"), pvalue.rename("pvalue")])
-
-
-def ks_2samp_1d(ds1, ds2, dim, kwargs={}):
-    """
-    The one-dimensional Kolmogorov-Smirnov test on two samples.
-
-    This test compares the underlying continuous distributions ds1 and ds2 of
-    two independent samples.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data. Nans are automatically removed prior to executing the test
-    ds2 : xarray Dataset
-        Sample 2 data. Nans are automatically removed prior to executing the test.
-        The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ks_2samp
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The KS statistic
-        - "pvalue" : One-tailed or two-tailed p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.ks_2samp.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-
-    See also
-    --------
-    scipy.stats.ks_2samp
-
-    References
-    ----------
-    Hodges, J.L. Jr., “The Significance Probability of the Smirnov Two-Sample Test,”
-        Arkiv fiur Matematik, 3, No. 43 (1958), 469-86.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def anderson_ksamp(*args, dim, kwargs={}):
-    """
-    The Anderson-Darling test for k independent samples.
-
-    The k-sample Anderson-Darling test is a modification of the one-sample
-    Anderson-Darling test. It tests the null hypothesis that k-samples are drawn
-    from the same population without having to specify the distribution function
-    of that population.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. Nans are automatically removed prior to executing the test.
-        The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ad_ksamp
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : Normalized k-sample Anderson-Darling test statistic
-        - "pvalue" : An approximate significance level at which the null hypothesis
-            for the provided samples can be rejected. The value is floored / capped
-            at 0.1% / 25%
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.anderson_ksamp.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def ttest_ind(ds1, ds2, dim, kwargs={}):
-    """
-    The T-test for the means of two independent samples.
-
-    This is a test for the null hypothesis that 2 independent samples have identical
-    average (expected) values. This test assumes that the populations have identical
-    variances by default.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ttest_ind
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The t-statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.ttest_ind.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def ttest_rel(ds1, ds2, dim, kwargs={}):
-    """
-    The T-test for the means of two related samples.
-
-    This is a test for the null hypothesis that two related or repeated samples
-    have identical average (expected) values.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim must be the same
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ttest_rel
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The t-statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.ttest_rel.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def cramervonmises_2samp(ds1, ds2, dim, kwargs={}):
-    """
-    The two-sample Cramér-von Mises test for goodness of fit.
-
-    This is the two-sample version of the Cramér-von Mises test: for two independent
-    samples,the null hypothesis is that the samples come from the same (unspecified)
-    continuous distribution.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data. Nans are automatically removed prior to executing the test
-    ds2 : xarray Dataset
-        Sample 2 data. Nans are automatically removed prior to executing the test.
-        The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.cramervonmises_2samp
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Cramér-von Mises statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.cramervonmises_2samp.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def epps_singleton_2samp(ds1, ds2, dim, kwargs={}):
-    """
-    The Epps-Singleton (ES) test on two samples.
-
-    This tests the null hypothesis that two samples have the same underlying probability
-    distribution.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data. Nans are automatically removed prior to executing the test
-    ds2 : xarray Dataset
-        Sample 2 data. Nans are automatically removed prior to executing the test.
-        The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.epps_singleton_2samp
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Epps-Singleton statistic
-        - "pvalue" : The p-value based on the asymptotic chi2-distribution
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.epps_singleton_2samp.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def mannwhitneyu(ds1, ds2, dim, kwargs={}):
-    """
-    The Mann-Whitney U rank test on two independent samples.
-
-    The Mann-Whitney U test is a nonparametric test of the null hypothesis that the
-    distribution underlying sample ds1 is the same as the distribution underlying
-    sample ds2. It is often used as a test of difference in location between distributions.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.mannwhitneyu
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Mann-Whitney U statistic corresponding with sample ds1. See the
-            scipy.stats.epps_singleton_2samp documentation for how to calculate the U
-            statistic corresponding with sample ds2
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.mannwhitneyu.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def ranksums(ds1, ds2, dim, kwargs={}):
-    """
-    The Wilcoxon rank-sum statistic for two independent samples.
-
-    The Wilcoxon rank-sum test tests the null hypothesis that two sets of measurements are drawn
-    from the same distribution. The alternative hypothesis is that values in one sample are more
-    likely to be larger than the values in the other sample. This test should be used to compare
-    two samples from continuous distributions. It does not handle ties between ds1 and ds2.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ranksums
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Wilcoxon rank-sum statistic under the large-sample approximation that
-            the rank sum statistic is normally distributed.
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.ranksums.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def kruskal(*args, dim, kwargs={}):
-    """
-    The Kruskal-Wallis H-test for k independent samples.
-
-    The Kruskal-Wallis H-test tests the null hypothesis that the population median of all of the
-    groups are equal. It is a non-parametric version of ANOVA.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.kruskal
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Kruskal-Wallis H statistic, corrected for ties.
-        - "pvalue" : The p-value for the test using the assumption that H has a chi square
-            distribution. The p-value returned is the survival function of the chi square
-            distribution evaluated at H.
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.kruskal.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def friedmanchisquare(*args, dim, kwargs={}):
-    """
-    The Friedman test for k repeated samples.
-
-    The Friedman test tests the null hypothesis that repeated samples of the same
-    individuals have the same distribution. It is often used to test for consistency
-    among samples obtained in different ways. For example, if two sampling techniques
-    are used on the same set of individuals, the Friedman test can be used to determine
-    if the two sampling techniques are consistent.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. Nans are automatically removed prior to executing the test.
-        The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.friedmanchisquare
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The test statistic, correcting for ties.
-        - "pvalue" : The p-value assuming that the test statistic has a chi squared
-            distribution.
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.kruskal.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def brunnermunzel(ds1, ds2, dim, kwargs={}):
-    """
-    The Brunner-Munzel test on two independent samples.
-
-    The Brunner-Munzel test is a nonparametric test of the null hypothesis that when values
-    are taken one by one from each group, the probabilities of getting large values in both
-    groups are equal. Unlike the Wilcoxon-Mann-Whitney’s U test, this does not require the
-    assumption of equivariance of two groups. Note that this does not assume the
-    distributions are same.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.brunnermunzel
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Brunner-Munzer W statistic.
-        - "pvalue" : The one-sided or two-sided p-value assuming an t distribution.
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.brunnermunzel.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def ansari(ds1, ds2, dim, kwargs={}):
-    """
-    The Ansari-Bradley test for equal distribution scale parameters from two independent samples.
-
-    The Ansari-Bradley test is a non-parametric test for the equality of the scale parameter
-    of the distributions from which two samples were drawn. The null hypothesis states that
-    the ratio of the scale of the distribution underlying ds1 to the scale of the distribution
-    underlying ds2 is 1.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data. Nans are automatically removed prior to executing the test
-    ds2 : xarray Dataset
-        Sample 2 data. Nans are automatically removed prior to executing the test.
-        The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.ansari
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Ansari-Bradley test statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.ansari.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
-
-
-def bartlett(*args, dim, kwargs={}):
-    """
-    The Bartlett test for the variances of k independent samples.
-
-    Bartlett’s test tests the null hypothesis that all input samples are from populations
-    with equal variances. For samples from significantly non-normal populations, Levene’s
-    test levene is more robust.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. Nans are automatically removed prior to executing the test.
-        The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.bartlett
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Bartlett test statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.bartlett.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def levene(*args, dim, kwargs={}):
-    """
-    The Levene test for the variances of k independent samples.
-
-    The Levene test tests the null hypothesis that all input samples are from populations
-    with equal variances. Levene’s test is an alternative to Bartlett’s test in the case
-    where there are significant deviations from normality.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. Nans are automatically removed prior to executing the test.
-        The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.levene
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Levene test statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.levene.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def fligner(*args, dim, kwargs={}):
-    """
-    The Fligner-Killeen test for the variances of k independent samples.
-
-    Fligner’s test tests the null hypothesis that all input samples are from populations
-    with equal variances. Fligner-Killeen’s test is distribution free when populations are
-    identical.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. Nans are automatically removed prior to executing the test.
-        The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.fligner
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The Fligner-Killeen test statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.fligner.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def median_test(*args, dim, kwargs={}):
-    """
-    The Mood test for the medians of k independent samples.
-
-    Mood's test tests that two or more samples come from populations with the same median.
-
-    Parameters
-    ----------
-    args : xarray Datasets
-        The k samples of data. The sizes of the samples along dim can be different
-    dim : str
-        The name of the sample dimension(s) in args
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.median_test
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The test statistic
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.median_test.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
-
-    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
-
-
-def mood(ds1, ds2, dim, kwargs={}):
-    """
-    The Mood test for equal distribution scale parameters from two independent samples.
-
-    Mood’s two-sample test for scale parameters is a non-parametric test for the null
-    hypothesis that two samples are drawn from the same distribution with the same scale
-    parameter.
-
-    Parameters
-    ----------
-    ds1 : xarray Dataset
-        Sample 1 data.
-    ds2 : xarray Dataset
-        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
-    dim : str
-        The name of the sample dimension(s) in ds1 and ds2
-    kwargs : dict
-        Any other kwargs to pass to scipy.stats.mood
-
-    Returns
-    -------
-    statistics : xarray Dataset
-        Dataset with the following variables:
-        - "statistic" : The z-score for the hypothesis test
-        - "pvalue" : The p-value
-
-    Notes
-    -----
-    This function is a simple wrapper on the scipy function scipy.stats.mood.
-    Users are recommended to read the scipy documentation prior to using this
-    function.
-    """
 
-    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+from .utils import _prep_data
 
 
 def ks_2samp_2d_np(x1, y1, x2, y2):
diff --git a/xstatstests/scipy.py b/xstatstests/scipy.py
new file mode 100644
index 0000000..535367a
--- /dev/null
+++ b/xstatstests/scipy.py
@@ -0,0 +1,1043 @@
+import inspect
+
+import numpy as np
+import xarray as xr
+import scipy.stats
+
+from .utils import _prep_data
+
+
+# Summary of how to wrap scipy.stats funcs
+# min_args is used only in test_scipy.py
+# Enter outputs so that first output is the test statistic and the second
+# output is the p-value
+# n_args = -1 means no limit to the number of args
+scipy_function_info = {
+    "ks_2samp_1d": {
+        "name": "ks_2samp",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "anderson_ksamp": {
+        "name": "anderson_ksamp",
+        "min_args": 2,
+        "stack_args": True,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": ["statistic", "significance_level"],
+    },
+    "ttest_ind": {
+        "name": "ttest_ind",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "ttest_rel": {
+        "name": "ttest_rel",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": True,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "cramervonmises_2samp": {
+        "name": "cramervonmises_2samp",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": ["statistic", "pvalue"],
+    },
+    "epps_singleton_2samp": {
+        "name": "epps_singleton_2samp",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "mannwhitneyu": {
+        "name": "mannwhitneyu",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis", "keepdims"],
+        "outputs": ["statistic", "pvalue"],
+    },
+    "ranksums": {
+        "name": "ranksums",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis", "keepdims"],
+        "outputs": [0, 1],
+    },
+    "kruskal": {
+        "name": "kruskal",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis", "keepdims"],
+        "outputs": [0, 1],
+    },
+    "friedmanchisquare": {
+        "name": "friedmanchisquare",
+        "min_args": 3,
+        "stack_args": False,
+        "same_sample_sizes": True,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "brunnermunzel": {
+        "name": "brunnermunzel",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "ansari": {
+        "name": "ansari",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "bartlett": {
+        "name": "bartlett",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "levene": {
+        "name": "levene",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "fligner": {
+        "name": "fligner",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "median_test": {
+        "name": "median_test",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+    "mood": {
+        "name": "mood",
+        "min_args": 2,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "skewtest": {
+        "name": "skewtest",
+        "min_args": 1,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "kurtosistest": {
+        "name": "kurtosistest",
+        "min_args": 1,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "normaltest": {
+        "name": "normaltest",
+        "min_args": 1,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": False,
+        "disallowed_kwargs": ["axis"],
+        "outputs": [0, 1],
+    },
+    "jarque_bera": {
+        "name": "jarque_bera",
+        "min_args": 1,
+        "stack_args": False,
+        "same_sample_sizes": False,
+        "remove_nans": True,
+        "disallowed_kwargs": [],
+        "outputs": [0, 1],
+    },
+}
+
+
+def _wrap_scipy(func, args, dim, kwargs):
+    """Generic xarray wrapper for subset of scipy stats functions"""
+
+    def _wrap_scipy_func(*args, scipy_func_info, scipy_kwargs):
+        """Parse scipy_function_info and apply scipy function"""
+        func = getattr(scipy.stats, scipy_func_info["name"])
+        getter = scipy_func_info["outputs"]
+
+        if scipy_func_info["remove_nans"]:
+            args = [arg[~np.isnan(arg)] for arg in args]
+        if scipy_func_info["stack_args"]:
+            outputs = func(args, **scipy_kwargs)
+        else:
+            outputs = func(*args, **scipy_kwargs)
+
+        return tuple(
+            [getattr(outputs, g) if isinstance(g, str) else outputs[g] for g in getter]
+        )
+
+    info = scipy_function_info[func]
+
+    args, input_core_dims = _prep_data(*args, dim=dim, nd=1)
+
+    sample_sizes = [ds.sizes[dim[0]] for ds, dim in zip(args, input_core_dims)]
+    same_sample_sizes = all([sample_sizes[0] == samp for samp in sample_sizes[1:]])
+    if info["same_sample_sizes"] & (not same_sample_sizes):
+        raise ValueError(
+            f"`{func}` requires that the sample size is the same for all input arrays"
+        )
+
+    scipy_kwargs = kwargs.copy()
+
+    # Simply vectorize if function does not allow axis argument
+    vectorize = True
+    for kw in info["disallowed_kwargs"]:
+        if kw in scipy_kwargs.keys():
+            raise ValueError(f"`{kw}` kwarg is disallowed by xstatstests")
+        if kw == "axis":
+            vectorize = False
+            scipy_kwargs["axis"] = -1
+
+    output_core_dims = [[]] * 2
+    output_dtypes = ["float32"] * 2
+    kwargs = dict(scipy_func_info=info, scipy_kwargs=scipy_kwargs)
+    statistic, pvalue = xr.apply_ufunc(
+        _wrap_scipy_func,
+        *args,
+        kwargs=kwargs,
+        input_core_dims=input_core_dims,
+        output_core_dims=output_core_dims,
+        output_dtypes=output_dtypes,
+        vectorize=vectorize,
+        dask="parallelized",
+    )
+
+    return xr.merge([statistic.rename("statistic"), pvalue.rename("pvalue")])
+
+
+def ks_2samp_1d(ds1, ds2, dim, kwargs={}):
+    """
+    The one-dimensional Kolmogorov-Smirnov test on two samples.
+
+    This test compares the underlying continuous distributions ds1 and ds2 of
+    two independent samples.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data. Nans are automatically removed prior to executing the test
+    ds2 : xarray Dataset
+        Sample 2 data. Nans are automatically removed prior to executing the test.
+        The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ks_2samp
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The KS statistic
+        - "pvalue" : One-tailed or two-tailed p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.ks_2samp.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+
+    See also
+    --------
+    scipy.stats.ks_2samp
+
+    References
+    ----------
+    Hodges, J.L. Jr., “The Significance Probability of the Smirnov Two-Sample Test,”
+        Arkiv fiur Matematik, 3, No. 43 (1958), 469-86.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def anderson_ksamp(*args, dim, kwargs={}):
+    """
+    The Anderson-Darling test for k independent samples.
+
+    The k-sample Anderson-Darling test is a modification of the one-sample
+    Anderson-Darling test. It tests the null hypothesis that k-samples are drawn
+    from the same population without having to specify the distribution function
+    of that population.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. Nans are automatically removed prior to executing the test.
+        The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ad_ksamp
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : Normalized k-sample Anderson-Darling test statistic
+        - "pvalue" : An approximate significance level at which the null hypothesis
+            for the provided samples can be rejected. The value is floored / capped
+            at 0.1% / 25%
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.anderson_ksamp.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def ttest_ind(ds1, ds2, dim, kwargs={}):
+    """
+    The T-test for the means of two independent samples.
+
+    This is a test for the null hypothesis that 2 independent samples have identical
+    average (expected) values. This test assumes that the populations have identical
+    variances by default.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ttest_ind
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The t-statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.ttest_ind.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def ttest_rel(ds1, ds2, dim, kwargs={}):
+    """
+    The T-test for the means of two related samples.
+
+    This is a test for the null hypothesis that two related or repeated samples
+    have identical average (expected) values.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim must be the same
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ttest_rel
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The t-statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.ttest_rel.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def cramervonmises_2samp(ds1, ds2, dim, kwargs={}):
+    """
+    The two-sample Cramér-von Mises test for goodness of fit.
+
+    This is the two-sample version of the Cramér-von Mises test: for two independent
+    samples,the null hypothesis is that the samples come from the same (unspecified)
+    continuous distribution.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data. Nans are automatically removed prior to executing the test
+    ds2 : xarray Dataset
+        Sample 2 data. Nans are automatically removed prior to executing the test.
+        The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.cramervonmises_2samp
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Cramér-von Mises statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.cramervonmises_2samp.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def epps_singleton_2samp(ds1, ds2, dim, kwargs={}):
+    """
+    The Epps-Singleton (ES) test on two samples.
+
+    This tests the null hypothesis that two samples have the same underlying probability
+    distribution.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data. Nans are automatically removed prior to executing the test
+    ds2 : xarray Dataset
+        Sample 2 data. Nans are automatically removed prior to executing the test.
+        The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.epps_singleton_2samp
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Epps-Singleton statistic
+        - "pvalue" : The p-value based on the asymptotic chi2-distribution
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.epps_singleton_2samp.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def mannwhitneyu(ds1, ds2, dim, kwargs={}):
+    """
+    The Mann-Whitney U rank test on two independent samples.
+
+    The Mann-Whitney U test is a nonparametric test of the null hypothesis that the
+    distribution underlying sample ds1 is the same as the distribution underlying
+    sample ds2. It is often used as a test of difference in location between distributions.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.mannwhitneyu
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Mann-Whitney U statistic corresponding with sample ds1. See the
+            scipy.stats.epps_singleton_2samp documentation for how to calculate the U
+            statistic corresponding with sample ds2
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.mannwhitneyu.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def ranksums(ds1, ds2, dim, kwargs={}):
+    """
+    The Wilcoxon rank-sum statistic for two independent samples.
+
+    The Wilcoxon rank-sum test tests the null hypothesis that two sets of measurements are drawn
+    from the same distribution. The alternative hypothesis is that values in one sample are more
+    likely to be larger than the values in the other sample. This test should be used to compare
+    two samples from continuous distributions. It does not handle ties between ds1 and ds2.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ranksums
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Wilcoxon rank-sum statistic under the large-sample approximation that
+            the rank sum statistic is normally distributed.
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.ranksums.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def kruskal(*args, dim, kwargs={}):
+    """
+    The Kruskal-Wallis H-test for k independent samples.
+
+    The Kruskal-Wallis H-test tests the null hypothesis that the population median of all of the
+    groups are equal. It is a non-parametric version of ANOVA.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.kruskal
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Kruskal-Wallis H statistic, corrected for ties.
+        - "pvalue" : The p-value for the test using the assumption that H has a chi square
+            distribution. The p-value returned is the survival function of the chi square
+            distribution evaluated at H.
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.kruskal.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def friedmanchisquare(*args, dim, kwargs={}):
+    """
+    The Friedman test for k repeated samples.
+
+    The Friedman test tests the null hypothesis that repeated samples of the same
+    individuals have the same distribution. It is often used to test for consistency
+    among samples obtained in different ways. For example, if two sampling techniques
+    are used on the same set of individuals, the Friedman test can be used to determine
+    if the two sampling techniques are consistent.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. Nans are automatically removed prior to executing the test.
+        The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.friedmanchisquare
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The test statistic, correcting for ties.
+        - "pvalue" : The p-value assuming that the test statistic has a chi squared
+            distribution.
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.kruskal.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def brunnermunzel(ds1, ds2, dim, kwargs={}):
+    """
+    The Brunner-Munzel test on two independent samples.
+
+    The Brunner-Munzel test is a nonparametric test of the null hypothesis that when values
+    are taken one by one from each group, the probabilities of getting large values in both
+    groups are equal. Unlike the Wilcoxon-Mann-Whitney’s U test, this does not require the
+    assumption of equivariance of two groups. Note that this does not assume the
+    distributions are same.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.brunnermunzel
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Brunner-Munzer W statistic.
+        - "pvalue" : The one-sided or two-sided p-value assuming an t distribution.
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.brunnermunzel.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def ansari(ds1, ds2, dim, kwargs={}):
+    """
+    The Ansari-Bradley test for equal distribution scale parameters from two independent samples.
+
+    The Ansari-Bradley test is a non-parametric test for the equality of the scale parameter
+    of the distributions from which two samples were drawn. The null hypothesis states that
+    the ratio of the scale of the distribution underlying ds1 to the scale of the distribution
+    underlying ds2 is 1.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data. Nans are automatically removed prior to executing the test
+    ds2 : xarray Dataset
+        Sample 2 data. Nans are automatically removed prior to executing the test.
+        The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.ansari
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Ansari-Bradley test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.ansari.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def bartlett(*args, dim, kwargs={}):
+    """
+    The Bartlett test for the variances of k independent samples.
+
+    Bartlett’s test tests the null hypothesis that all input samples are from populations
+    with equal variances. For samples from significantly non-normal populations, Levene’s
+    test levene is more robust.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. Nans are automatically removed prior to executing the test.
+        The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.bartlett
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Bartlett test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.bartlett.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def levene(*args, dim, kwargs={}):
+    """
+    The Levene test for the variances of k independent samples.
+
+    The Levene test tests the null hypothesis that all input samples are from populations
+    with equal variances. Levene’s test is an alternative to Bartlett’s test in the case
+    where there are significant deviations from normality.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. Nans are automatically removed prior to executing the test.
+        The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.levene
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Levene test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.levene.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def fligner(*args, dim, kwargs={}):
+    """
+    The Fligner-Killeen test for the variances of k independent samples.
+
+    Fligner’s test tests the null hypothesis that all input samples are from populations
+    with equal variances. Fligner-Killeen’s test is distribution free when populations are
+    identical.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. Nans are automatically removed prior to executing the test.
+        The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.fligner
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Fligner-Killeen test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.fligner.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def median_test(*args, dim, kwargs={}):
+    """
+    The Mood test for the medians of k independent samples.
+
+    Mood's test tests that two or more samples come from populations with the same median.
+
+    Parameters
+    ----------
+    args : xarray Datasets
+        The k samples of data. The sizes of the samples along dim can be different
+    dim : str
+        The name of the sample dimension(s) in args
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.median_test
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.median_test.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], args, dim, kwargs)
+
+
+def mood(ds1, ds2, dim, kwargs={}):
+    """
+    The Mood test for equal distribution scale parameters from two independent samples.
+
+    Mood’s two-sample test for scale parameters is a non-parametric test for the null
+    hypothesis that two samples are drawn from the same distribution with the same scale
+    parameter.
+
+    Parameters
+    ----------
+    ds1 : xarray Dataset
+        Sample 1 data.
+    ds2 : xarray Dataset
+        Sample 2 data. The sizes of samples 1 and 2 along dim can be different
+    dim : str
+        The name of the sample dimension(s) in ds1 and ds2
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.mood
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The z-score for the hypothesis test
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.mood.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds1, ds2], dim, kwargs)
+
+
+def skewtest(ds, dim, kwargs={}):
+    """
+    The D’Agostino test for normal sample skewness.
+
+    This function tests the null hypothesis that the skewness of the population
+    that the sample was drawn from is the same as that of a corresponding normal
+    distribution.
+
+    Parameters
+    ----------
+    ds : xarray Dataset
+        Sample data.
+    dim : str
+        The name of the sample dimension(s) in ds
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.skewtest
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The z-score for the hypothesis test
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.skewtest.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+
+    References
+    ----------
+    - R. B. D’Agostino, A. J. Belanger and R. B. D’Agostino Jr., “A suggestion for
+        using powerful and informative tests of normality”, American Statistician
+        44, pp. 316-321, 1990.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds], dim, kwargs)
+
+
+def kurtosistest(ds, dim, kwargs={}):
+    """
+    The Anscombe test for normal sample kurtosis.
+
+    This function tests the null hypothesis that the kurtosis of the population
+    from which the sample was drawn is that of the normal distribution.
+
+    Parameters
+    ----------
+    ds : xarray Dataset
+        Sample data.
+    dim : str
+        The name of the sample dimension(s) in ds
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.kurtosistest
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The z-score for the hypothesis test
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.kurtosistest.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+
+    References
+    ----------
+    - F. J. Anscombe, W. J. Glynn, “Distribution of the kurtosis statistic b2 for
+        normal samples”, Biometrika, vol. 70, pp. 227-234, 1983.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds], dim, kwargs)
+
+
+def normaltest(ds, dim, kwargs={}):
+    """
+    The D’Agostino/Pearson test for whether a sample differs from a normal distribution.
+
+    This function tests the null hypothesis that a sample comes from a normal
+    distribution. It is based on D’Agostino and Pearson’s test that combines skewness and
+    kurtosis to produce an omnibus test of normality.
+
+    Parameters
+    ----------
+    ds : xarray Dataset
+        Sample data.
+    dim : str
+        The name of the sample dimension(s) in ds
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.normaltest
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : s^2 + k^2, where s is the z-score returned by xstatstests.skewtest
+            and k is the z-score returned by xstatstests.kurtosistest.
+        - "pvalue" : The 2-sided chi squared p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.normaltest.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+
+    References
+    ----------
+    - D’Agostino, R. B. (1971), “An omnibus test of normality for moderate and large sample
+        size”, Biometrika, 58, 341-348
+    - D’Agostino, R. and Pearson, E. S. (1973), “Tests for departure from normality”,
+        Biometrika, 60, 613-622
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds], dim, kwargs)
+
+
+def jarque_bera(ds, dim, kwargs={}):
+    """
+    The Jarque-Bera goodness of fit test on sample data.
+
+    The Jarque-Bera test tests whether the sample data has the skewness and kurtosis matching
+    a normal distribution. Note that this test only works for a large enough number of data
+    samples (>2000) as the test statistic asymptotically has a Chi-squared distribution with 2
+    degrees of freedom.
+
+    Parameters
+    ----------
+    ds : xarray Dataset
+        Sample data. Nans are automatically removed prior to executing the test.
+    dim : str
+        The name of the sample dimension(s) in ds
+    kwargs : dict
+        Any other kwargs to pass to scipy.stats.jarque_bera
+
+    Returns
+    -------
+    statistics : xarray Dataset
+        Dataset with the following variables:
+        - "statistic" : The Jarque-Bera test statistic
+        - "pvalue" : The p-value
+
+    Notes
+    -----
+    This function is a simple wrapper on the scipy function scipy.stats.jarque_bera.
+    Users are recommended to read the scipy documentation prior to using this
+    function.
+    """
+
+    return _wrap_scipy(inspect.stack()[0][3], [ds], dim, kwargs)
diff --git a/xstatstests/tests/test_core.py b/xstatstests/tests/test_core.py
index 0696589..e57e6b2 100755
--- a/xstatstests/tests/test_core.py
+++ b/xstatstests/tests/test_core.py
@@ -1,451 +1,13 @@
 import pytest
 
-import numpy as np
 import numpy.testing as npt
 
 import xarray as xr
 
-import scipy.stats
-
-from xstatstests.core import scipy_function_info
 import xstatstests as xst
 from .fixtures import ds_1var
 
 
-def check_vs_scipy_func(func, args, kwargs={}):
-    """Test values relative to scipy function"""
-
-    def _stack_sample_dim(ds):
-        """Stack sample dim into two dimensions"""
-        new_length = int(ds.sizes["sample"] / 2)
-        return (
-            ds.assign_coords(sample_1=range(2), sample_2=range(new_length))
-            .stack(dim=["sample_1", "sample_2"])
-            .reset_index("sample", drop=True)
-            .rename(sample="dim")
-            .unstack("dim")
-        )
-
-    def _test_vs_scipy_values(inputs, outputs, func_info, kwargs={}):
-        """Test wrapped xstatstests func values relative to scipy"""
-        scipy_func = getattr(scipy.stats, func_info["name"])
-
-        inputs_np_1d = [
-            np.reshape(inp["var"].values, (inp.sizes["sample"], -1))[:, 0]
-            for inp in inputs
-        ]
-        outputs_np = [outputs["statistic"].values, outputs["pvalue"].values]
-
-        if func_info["stack_args"]:
-            scipy_outputs = scipy_func(inputs_np_1d, **kwargs)
-        else:
-            scipy_outputs = scipy_func(*inputs_np_1d, **kwargs)
-
-        getter = func_info["outputs"]
-        outputs_ver = [
-            getattr(scipy_outputs, g) if isinstance(g, str) else scipy_outputs[g]
-            for g in getter
-        ]
-
-        for res, ver in zip(outputs_np, outputs_ver):
-            npt.assert_allclose(res, ver)
-
-    function_info = scipy_function_info[func]
-
-    # Test with a single sample dim
-
-    outputs = getattr(xst, func)(*args, dim="sample", kwargs=kwargs)
-    _test_vs_scipy_values(args, outputs, function_info, kwargs=kwargs)
-
-    # Test with multiple sample dims
-    args_stack = [_stack_sample_dim(ds) for ds in args]
-    outputs_stack = getattr(xst, func)(
-        *args_stack, dim=["sample_1", "sample_2"], kwargs=kwargs
-    )
-    _test_vs_scipy_values(args, outputs_stack, function_info, kwargs=kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-@pytest.mark.parametrize("method", ["auto", "exact", "asymp"])
-def test_ks_2samp_1d_values(
-    ds1_n_per_sample, ds2_n_per_sample, shape, alternative, method
-):
-    """Check ks_2samp_1d relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
-    ]
-    kwargs = dict(alternative=alternative, method=method)
-    check_vs_scipy_func("ks_2samp_1d", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("midrank", [True, False])
-def test_anderson_ksamp_values(k_samples, n_per_sample, shape, midrank):
-    """Check anderson_ksamp relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=False, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    kwargs = dict(midrank=midrank)
-    check_vs_scipy_func("anderson_ksamp", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("equal_var", [True, False])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-@pytest.mark.parametrize("permutations", [0, 1000])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-@pytest.mark.parametrize("trim", [0, 0.4])
-def test_ttest_ind_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    add_nans,
-    equal_var,
-    nan_policy,
-    permutations,
-    alternative,
-    trim,
-):
-    """Check ttest_ind relative to scipy func"""
-    if (nan_policy == "omit") & ((permutations != 0) | (trim != 0)):
-        pytest.skip(
-            "nan_policy='omit' is currently not supported by permutation tests or trimmed tests."
-        )
-    elif (permutations != 0) & (trim != 0):
-        pytest.skip("Permutations are currently not supported with trimming.")
-    else:
-        args = [
-            ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-            ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-        ]
-        kwargs = dict(
-            equal_var=equal_var,
-            nan_policy=nan_policy,
-            permutations=permutations,
-            random_state=0,
-            alternative=alternative,
-            trim=trim,
-        )
-        check_vs_scipy_func("ttest_ind", args, kwargs)
-
-
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-def test_ttest_rel_values(
-    shape,
-    add_nans,
-    nan_policy,
-    alternative,
-):
-    """Check ttest_rel relative to scipy func"""
-    ds1_n_per_sample = ds2_n_per_sample = 10
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-    ]
-    kwargs = dict(
-        nan_policy=nan_policy,
-        alternative=alternative,
-    )
-    check_vs_scipy_func("ttest_rel", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [2, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [2, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("method", ["auto", "exact", "asymptotic"])
-def test_cramervonmises_2samp_values(ds1_n_per_sample, ds2_n_per_sample, shape, method):
-    """Check cramervonmises_2samp relative to scipy func"""
-    if (method == "exact") & ((ds1_n_per_sample > 2) | (ds2_n_per_sample > 2)):
-        pytest.skip("method == 'exact' is very slow for large samples.")
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
-    ]
-    kwargs = dict(method=method)
-    check_vs_scipy_func("cramervonmises_2samp", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("t", [(0.4, 0.8), (0.25, 0.5, 0.75)])
-def test_epps_singleton_2samp_values(ds1_n_per_sample, ds2_n_per_sample, shape, t):
-    """Check epps_singleton_2samp relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
-    ]
-    kwargs = dict(t=t)
-    check_vs_scipy_func("epps_singleton_2samp", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("use_continuity", [True, False])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-@pytest.mark.parametrize("method", ["auto", "exact", "asymptotic"])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-def test_mannwhitneyu_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    add_nans,
-    use_continuity,
-    alternative,
-    method,
-    nan_policy,
-):
-    """Check mannwhitneyu relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-    ]
-    kwargs = dict(
-        use_continuity=use_continuity,
-        alternative=alternative,
-        method=method,
-        nan_policy=nan_policy,
-    )
-    check_vs_scipy_func("mannwhitneyu", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-def test_ranksums_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    add_nans,
-    alternative,
-    nan_policy,
-):
-    """Check ranksums relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-    ]
-    kwargs = dict(
-        alternative=alternative,
-        nan_policy=nan_policy,
-    )
-    check_vs_scipy_func("ranksums", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-def test_kruskal_values(
-    k_samples,
-    n_per_sample,
-    shape,
-    add_nans,
-    nan_policy,
-):
-    """Check kruskal relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=add_nans, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    kwargs = dict(
-        nan_policy=nan_policy,
-    )
-    check_vs_scipy_func("kruskal", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [3, 4, 5])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-def test_friedmanchisquare_values(
-    k_samples,
-    shape,
-):
-    """Check friedmanchisquare relative to scipy func"""
-    n_per_sample = [10, 10, 10, 10, 10]
-    args = [
-        ds_1var((n,) + shape, add_nans=False, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    check_vs_scipy_func("friedmanchisquare", args)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-@pytest.mark.parametrize("distribution", ["t", "normal"])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-def test_brunnermunzel_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    add_nans,
-    alternative,
-    distribution,
-    nan_policy,
-):
-    """Check brunnermunzel relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
-    ]
-    kwargs = dict(
-        alternative=alternative,
-        distribution=distribution,
-        nan_policy=nan_policy,
-    )
-    check_vs_scipy_func("brunnermunzel", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-def test_ansari_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    alternative,
-):
-    """Check ansari relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
-    ]
-    kwargs = dict(
-        alternative=alternative,
-    )
-    check_vs_scipy_func("ansari", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-def test_bartlett_values(
-    k_samples,
-    n_per_sample,
-    shape,
-):
-    """Check bartlett relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=False, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    check_vs_scipy_func("bartlett", args)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("center", ["mean", "median", "trimmed"])
-@pytest.mark.parametrize("proportiontocut", [0.01, 0.05, 0.1])
-def test_levene_values(k_samples, n_per_sample, shape, center, proportiontocut):
-    """Check levene relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=False, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    kwargs = dict(
-        center=center,
-        proportiontocut=proportiontocut,
-    )
-    check_vs_scipy_func("levene", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("center", ["mean", "median", "trimmed"])
-@pytest.mark.parametrize("proportiontocut", [0.01, 0.05, 0.1])
-def test_fligner_values(k_samples, n_per_sample, shape, center, proportiontocut):
-    """Check fligner relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=False, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    kwargs = dict(
-        center=center,
-        proportiontocut=proportiontocut,
-    )
-    check_vs_scipy_func("fligner", args, kwargs)
-
-
-@pytest.mark.parametrize("k_samples", [2, 3, 5])
-@pytest.mark.parametrize(
-    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
-)
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("add_nans", [True, False])
-@pytest.mark.parametrize("ties", ["below", "above", "ignore"])
-@pytest.mark.parametrize("correction", [True, False])
-@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
-def test_median_test_values(
-    k_samples, n_per_sample, shape, add_nans, ties, correction, nan_policy
-):
-    """Check median_test relative to scipy func"""
-    args = [
-        ds_1var((n,) + shape, add_nans=add_nans, dask=False)
-        for n in n_per_sample[slice(k_samples)]
-    ]
-    kwargs = dict(
-        ties=ties,
-        correction=correction,
-        nan_policy=nan_policy,
-    )
-    check_vs_scipy_func("median_test", args, kwargs)
-
-
-@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
-@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
-@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
-@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
-def test_mood_values(
-    ds1_n_per_sample,
-    ds2_n_per_sample,
-    shape,
-    alternative,
-):
-    """Check mood relative to scipy func"""
-    args = [
-        ds_1var((ds1_n_per_sample,) + shape, add_nans=True, dask=False),
-        ds_1var((ds2_n_per_sample,) + shape, add_nans=True, dask=False),
-    ]
-    kwargs = dict(
-        alternative=alternative,
-    )
-    check_vs_scipy_func("mood", args, kwargs)
-
-
 @pytest.mark.parametrize("samples", [10, 50])
 @pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
 @pytest.mark.parametrize("dask", [True, False])
@@ -459,45 +21,3 @@ def test_ks_2samp_2d_identical(samples, shape, dask):
 
     if dask is False:
         npt.assert_allclose(D["statistic"].values, 0.0)
-
-
-@pytest.mark.parametrize(
-    "func", ["ttest_ind", "ttest_rel", "mannwhitneyu", "ranksums", "kruskal", "mood"]
-)
-@pytest.mark.parametrize("dask", [True, False])
-def test_disallowed_error(func, dask):
-    """Check that error is thrown when a disallowed kwarg is provided to scipy funcs"""
-    n_per_sample = [10, 10]
-    shape = (2, 3)
-    args = [ds_1var((n,) + shape, dask) for n in n_per_sample]
-    kws = scipy_function_info[func]["disallowed_kwargs"]
-    for kw in kws:
-        with pytest.raises(ValueError):
-            getattr(xst, func)(*args, dim="sample", kwargs={kw: None})
-
-
-@pytest.mark.parametrize("func", scipy_function_info.keys())
-def test_dask_compute(func):
-    """Check that functions run with dask arrays and don't compute"""
-    if func == "friedmanchisquare":
-        n_samples = 3
-    else:
-        n_samples = 2
-    n_per_sample = n_samples * [10]
-    shape = (2, 3)
-    args = [ds_1var((n,) + shape, True) for n in n_per_sample]
-    getattr(xst, func)(*args, dim="sample")
-
-
-@pytest.mark.parametrize("func", ["ttest_rel", "friedmanchisquare"])
-@pytest.mark.parametrize("dask", [True, False])
-def test_sample_size_error(func, dask):
-    """Check that error is thrown when sample sizes differ for functions that don't allow this"""
-    if func == "friedmanchisquare":
-        n_per_sample = [10, 20, 30]
-    else:
-        n_per_sample = [10, 20]
-    shape = (2, 3)
-    args = [ds_1var((n,) + shape, dask) for n in n_per_sample]
-    with pytest.raises(ValueError):
-        getattr(xst, func)(*args, dim="sample")
diff --git a/xstatstests/tests/test_scipy.py b/xstatstests/tests/test_scipy.py
new file mode 100644
index 0000000..3984708
--- /dev/null
+++ b/xstatstests/tests/test_scipy.py
@@ -0,0 +1,560 @@
+import pytest
+
+import numpy as np
+import numpy.testing as npt
+
+import scipy.stats
+
+import xstatstests as xst
+from xstatstests.scipy import scipy_function_info
+from .fixtures import ds_1var
+
+
+def check_vs_scipy_func(func, args, kwargs={}):
+    """Test values relative to scipy function"""
+
+    def _stack_sample_dim(ds):
+        """Stack sample dim into two dimensions"""
+        new_length = int(ds.sizes["sample"] / 2)
+        return (
+            ds.assign_coords(sample_1=range(2), sample_2=range(new_length))
+            .stack(dim=["sample_1", "sample_2"])
+            .reset_index("sample", drop=True)
+            .rename(sample="dim")
+            .unstack("dim")
+        )
+
+    def _test_vs_scipy_values(inputs, outputs, func_info, kwargs={}):
+        """Test wrapped xstatstests func values relative to scipy"""
+        scipy_func = getattr(scipy.stats, func_info["name"])
+
+        inputs_np_1d = [
+            np.reshape(inp["var"].values, (inp.sizes["sample"], -1))[:, 0]
+            for inp in inputs
+        ]
+        outputs_np = [outputs["statistic"].values, outputs["pvalue"].values]
+
+        if func_info["stack_args"]:
+            scipy_outputs = scipy_func(inputs_np_1d, **kwargs)
+        else:
+            scipy_outputs = scipy_func(*inputs_np_1d, **kwargs)
+
+        getter = func_info["outputs"]
+        outputs_ver = [
+            getattr(scipy_outputs, g) if isinstance(g, str) else scipy_outputs[g]
+            for g in getter
+        ]
+
+        for res, ver in zip(outputs_np, outputs_ver):
+            npt.assert_allclose(res, ver)
+
+    function_info = scipy_function_info[func]
+
+    # Test with a single sample dim
+
+    outputs = getattr(xst, func)(*args, dim="sample", kwargs=kwargs)
+    _test_vs_scipy_values(args, outputs, function_info, kwargs=kwargs)
+
+    # Test with multiple sample dims
+    args_stack = [_stack_sample_dim(ds) for ds in args]
+    outputs_stack = getattr(xst, func)(
+        *args_stack, dim=["sample_1", "sample_2"], kwargs=kwargs
+    )
+    _test_vs_scipy_values(args, outputs_stack, function_info, kwargs=kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+@pytest.mark.parametrize("method", ["auto", "exact", "asymp"])
+def test_ks_2samp_1d_values(
+    ds1_n_per_sample, ds2_n_per_sample, shape, alternative, method
+):
+    """Check ks_2samp_1d relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    kwargs = dict(alternative=alternative, method=method)
+    check_vs_scipy_func("ks_2samp_1d", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("midrank", [True, False])
+def test_anderson_ksamp_values(k_samples, n_per_sample, shape, midrank):
+    """Check anderson_ksamp relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=False, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    kwargs = dict(midrank=midrank)
+    check_vs_scipy_func("anderson_ksamp", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("equal_var", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+@pytest.mark.parametrize("permutations", [0, 1000])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+@pytest.mark.parametrize("trim", [0, 0.4])
+def test_ttest_ind_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    add_nans,
+    equal_var,
+    nan_policy,
+    permutations,
+    alternative,
+    trim,
+):
+    """Check ttest_ind relative to scipy func"""
+    if (nan_policy == "omit") & ((permutations != 0) | (trim != 0)):
+        pytest.skip(
+            "nan_policy='omit' is currently not supported by permutation tests or trimmed tests."
+        )
+    elif (permutations != 0) & (trim != 0):
+        pytest.skip("Permutations are currently not supported with trimming.")
+    else:
+        args = [
+            ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+            ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+        ]
+        kwargs = dict(
+            equal_var=equal_var,
+            nan_policy=nan_policy,
+            permutations=permutations,
+            random_state=0,
+            alternative=alternative,
+            trim=trim,
+        )
+        check_vs_scipy_func("ttest_ind", args, kwargs)
+
+
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+def test_ttest_rel_values(
+    shape,
+    add_nans,
+    nan_policy,
+    alternative,
+):
+    """Check ttest_rel relative to scipy func"""
+    ds1_n_per_sample = ds2_n_per_sample = 10
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        nan_policy=nan_policy,
+        alternative=alternative,
+    )
+    check_vs_scipy_func("ttest_rel", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [2, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [2, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("method", ["auto", "exact", "asymptotic"])
+def test_cramervonmises_2samp_values(ds1_n_per_sample, ds2_n_per_sample, shape, method):
+    """Check cramervonmises_2samp relative to scipy func"""
+    if (method == "exact") & ((ds1_n_per_sample > 2) | (ds2_n_per_sample > 2)):
+        pytest.skip("method == 'exact' is very slow for large samples.")
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    kwargs = dict(method=method)
+    check_vs_scipy_func("cramervonmises_2samp", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("t", [(0.4, 0.8), (0.25, 0.5, 0.75)])
+def test_epps_singleton_2samp_values(ds1_n_per_sample, ds2_n_per_sample, shape, t):
+    """Check epps_singleton_2samp relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    kwargs = dict(t=t)
+    check_vs_scipy_func("epps_singleton_2samp", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("use_continuity", [True, False])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+@pytest.mark.parametrize("method", ["auto", "exact", "asymptotic"])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_mannwhitneyu_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    add_nans,
+    use_continuity,
+    alternative,
+    method,
+    nan_policy,
+):
+    """Check mannwhitneyu relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        use_continuity=use_continuity,
+        alternative=alternative,
+        method=method,
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("mannwhitneyu", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_ranksums_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    add_nans,
+    alternative,
+    nan_policy,
+):
+    """Check ranksums relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        alternative=alternative,
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("ranksums", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_kruskal_values(
+    k_samples,
+    n_per_sample,
+    shape,
+    add_nans,
+    nan_policy,
+):
+    """Check kruskal relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=add_nans, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    kwargs = dict(
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("kruskal", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [3, 4, 5])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+def test_friedmanchisquare_values(
+    k_samples,
+    shape,
+):
+    """Check friedmanchisquare relative to scipy func"""
+    n_per_sample = [10, 10, 10, 10, 10]
+    args = [
+        ds_1var((n,) + shape, add_nans=False, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    check_vs_scipy_func("friedmanchisquare", args)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+@pytest.mark.parametrize("distribution", ["t", "normal"])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_brunnermunzel_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    add_nans,
+    alternative,
+    distribution,
+    nan_policy,
+):
+    """Check brunnermunzel relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        alternative=alternative,
+        distribution=distribution,
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("brunnermunzel", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+def test_ansari_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    alternative,
+):
+    """Check ansari relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    kwargs = dict(
+        alternative=alternative,
+    )
+    check_vs_scipy_func("ansari", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+def test_bartlett_values(
+    k_samples,
+    n_per_sample,
+    shape,
+):
+    """Check bartlett relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=False, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    check_vs_scipy_func("bartlett", args)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("center", ["mean", "median", "trimmed"])
+@pytest.mark.parametrize("proportiontocut", [0.01, 0.05, 0.1])
+def test_levene_values(k_samples, n_per_sample, shape, center, proportiontocut):
+    """Check levene relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=False, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    kwargs = dict(
+        center=center,
+        proportiontocut=proportiontocut,
+    )
+    check_vs_scipy_func("levene", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("center", ["mean", "median", "trimmed"])
+@pytest.mark.parametrize("proportiontocut", [0.01, 0.05, 0.1])
+def test_fligner_values(k_samples, n_per_sample, shape, center, proportiontocut):
+    """Check fligner relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=False, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    kwargs = dict(
+        center=center,
+        proportiontocut=proportiontocut,
+    )
+    check_vs_scipy_func("fligner", args, kwargs)
+
+
+@pytest.mark.parametrize("k_samples", [2, 3, 5])
+@pytest.mark.parametrize(
+    "n_per_sample", [[10, 10, 10, 10, 10], [10, 20, 30, 40, 50], [50, 40, 30, 20, 10]]
+)
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("ties", ["below", "above", "ignore"])
+@pytest.mark.parametrize("correction", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_median_test_values(
+    k_samples, n_per_sample, shape, add_nans, ties, correction, nan_policy
+):
+    """Check median_test relative to scipy func"""
+    args = [
+        ds_1var((n,) + shape, add_nans=add_nans, dask=False)
+        for n in n_per_sample[slice(k_samples)]
+    ]
+    kwargs = dict(
+        ties=ties,
+        correction=correction,
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("median_test", args, kwargs)
+
+
+@pytest.mark.parametrize("ds1_n_per_sample", [10, 30])
+@pytest.mark.parametrize("ds2_n_per_sample", [10, 20])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+def test_mood_values(
+    ds1_n_per_sample,
+    ds2_n_per_sample,
+    shape,
+    alternative,
+):
+    """Check mood relative to scipy func"""
+    args = [
+        ds_1var((ds1_n_per_sample,) + shape, add_nans=False, dask=False),
+        ds_1var((ds2_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    kwargs = dict(
+        alternative=alternative,
+    )
+    check_vs_scipy_func("mood", args, kwargs)
+
+
+@pytest.mark.parametrize("ds_n_per_sample", [10, 30])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+def test_skewtest_values(
+    ds_n_per_sample,
+    shape,
+    add_nans,
+    nan_policy,
+    alternative,
+):
+    """Check skewtest relative to scipy func"""
+    args = [
+        ds_1var((ds_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        nan_policy=nan_policy,
+        alternative=alternative,
+    )
+    check_vs_scipy_func("skewtest", args, kwargs)
+
+
+@pytest.mark.parametrize("ds_n_per_sample", [20, 30])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+@pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
+def test_kurtosistest_values(
+    ds_n_per_sample,
+    shape,
+    add_nans,
+    nan_policy,
+    alternative,
+):
+    """Check kurtosistest relative to scipy func"""
+    args = [
+        ds_1var((ds_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        nan_policy=nan_policy,
+        alternative=alternative,
+    )
+    check_vs_scipy_func("kurtosistest", args, kwargs)
+
+
+@pytest.mark.parametrize("ds_n_per_sample", [20, 30])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+@pytest.mark.parametrize("add_nans", [True, False])
+@pytest.mark.parametrize("nan_policy", ["propagate", "omit"])
+def test_normaltest_values(
+    ds_n_per_sample,
+    shape,
+    add_nans,
+    nan_policy,
+):
+    """Check normaltest relative to scipy func"""
+    args = [
+        ds_1var((ds_n_per_sample,) + shape, add_nans=add_nans, dask=False),
+    ]
+    kwargs = dict(
+        nan_policy=nan_policy,
+    )
+    check_vs_scipy_func("normaltest", args, kwargs)
+
+
+@pytest.mark.parametrize("ds_n_per_sample", [100, 2000])
+@pytest.mark.parametrize("shape", [(), (2,), (2, 3)])
+def test_jarque_bera_values(
+    ds_n_per_sample,
+    shape,
+):
+    """Check jarque_bera relative to scipy func"""
+    args = [
+        ds_1var((ds_n_per_sample,) + shape, add_nans=False, dask=False),
+    ]
+    check_vs_scipy_func("jarque_bera", args)
+
+
+@pytest.mark.parametrize(
+    "func", ["ttest_ind", "ttest_rel", "mannwhitneyu", "ranksums", "kruskal", "mood"]
+)
+@pytest.mark.parametrize("dask", [True, False])
+def test_disallowed_error(func, dask):
+    """Check that error is thrown when a disallowed kwarg is provided to scipy funcs"""
+    n_per_sample = [10, 10]
+    shape = (2, 3)
+    args = [ds_1var((n,) + shape, dask) for n in n_per_sample]
+    kws = scipy_function_info[func]["disallowed_kwargs"]
+    for kw in kws:
+        with pytest.raises(ValueError):
+            getattr(xst, func)(*args, dim="sample", kwargs={kw: None})
+
+
+@pytest.mark.parametrize("func", scipy_function_info.keys())
+def test_dask_compute(func):
+    """Check that functions run with dask arrays and don't compute"""
+    n_per_sample = [10, 10, 10]
+    n_per_sample = n_per_sample[: scipy_function_info[func]["min_args"]]
+    shape = (2, 3)
+    args = [ds_1var((n,) + shape, True) for n in n_per_sample]
+    getattr(xst, func)(*args, dim="sample")
+
+
+@pytest.mark.parametrize("func", ["ttest_rel", "friedmanchisquare"])
+@pytest.mark.parametrize("dask", [True, False])
+def test_sample_size_error(func, dask):
+    """Check that error is thrown when sample sizes differ for functions that don't allow this"""
+    n_per_sample = [10, 20, 30]
+    n_per_sample = n_per_sample[: scipy_function_info[func]["min_args"]]
+    shape = (2, 3)
+    args = [ds_1var((n,) + shape, dask) for n in n_per_sample]
+    with pytest.raises(ValueError):
+        getattr(xst, func)(*args, dim="sample")
diff --git a/xstatstests/utils.py b/xstatstests/utils.py
new file mode 100644
index 0000000..918a20c
--- /dev/null
+++ b/xstatstests/utils.py
@@ -0,0 +1,43 @@
+import xarray as xr
+
+SAMPLE_DIM = "xsampletest_sample_dim"
+
+
+def _prep_data(*args, dim, nd):
+    """Prepare data for 2D tests"""
+    if isinstance(dim, str):
+        dim = [dim]
+
+    if any([(not isinstance(ds, xr.Dataset)) for ds in args]):
+        raise TypeError(
+            f"Input arrays must be xarray Datasets with {nd} variable(s) each"
+        )
+
+    args = xr.broadcast(*[ds.copy() for ds in args], exclude=dim)
+
+    if len(dim) == 1:
+        args = [ds.rename({dim[0]: SAMPLE_DIM}) for ds in args]
+    else:
+        args = [ds.stack({SAMPLE_DIM: dim}) for ds in args]
+
+    args = [ds.drop_vars({SAMPLE_DIM, *dim}, errors="ignore") for ds in args]
+    args = [ds.assign_coords({SAMPLE_DIM: range(ds.sizes[SAMPLE_DIM])}) for ds in args]
+
+    assert all(
+        [len(ds.data_vars) == nd for ds in args]
+    ), f"Input Datasets must have {nd} variables each"
+    assert all(
+        [list(args[0].data_vars) == list(ds.data_vars) for ds in args]
+    ), "Variables of all input Datasets must have the same name(s)"
+
+    # Need to rename sample dim otherwise apply_ufunc tries to align
+    # Expand into list of single variable Datasets
+    args_prepped = []
+    input_core_dims = []
+    for ind, ds in enumerate(args):
+        for var in ds.data_vars:
+            sample_dim = f"{SAMPLE_DIM}{ind+1}"
+            input_core_dims.append([sample_dim])
+            args_prepped.append(ds[var].rename({SAMPLE_DIM: sample_dim}))
+
+    return args_prepped, input_core_dims