Repository: dougiesquire/xstatstests (language: python)

Issue:
ks2d2s should work lazily on dask arrays
I thought it did, but it doesn't. Requires some refactoring of the internals


Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 9048144..b0365ac 100755
--- a/README.md
+++ b/README.md
@@ -6,10 +6,31 @@
 [![License: MIT](https://img.shields.io/badge/License-MIT-green.svg)](https://github.com/dougiesquire/xsampletests/blob/master/LICENSE)
 [![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/python/black)
 
-Sample testing for xarray objects. Currently includes the following statistical tests (most are currently simple, often vectorized, wrappers of [`scipy.stats`](https://docs.scipy.org/doc/scipy/reference/stats.html) functions):
-- `xsampletests.ks_2samp_1d`: one-dimensional Kolmogorov-Smirnov test on two samples
-- `xsampletests.ks_2samp_2d`: two-dimensional Kolmogorov-Smirnov test on two samples
-- `xsampletests.anderson_ksamp`: Anderson-Darling test on K samples
+Statistical test on multiple samples stored in xarray objects. Currently includes the following statistical tests (most are currently simple, often vectorized, wrappers of [`scipy.stats`](https://docs.scipy.org/doc/scipy/reference/stats.html) functions):
+
+#### Distribution tests
+- `xsampletests.ks_2samp_1d`: one-dimensional Kolmogorov-Smirnov test on two independent samples
+- `xsampletests.ks_2samp_2d`: two-dimensional Kolmogorov-Smirnov test on two independent samples
+- `xsampletests.anderson_ksamp`: Anderson-Darling test on K independent samples
+- `xsampletests.cramervonmises_2samp` : Cramér-von Mises test on two independent samples
+- `xsampletests.epps_singleton_2samp` : Epps-Singleton test on two independent samples
+- `xsampletests.ansari` : Ansari-Bradley test for equal distribution scale parameters from two independent samples
+- `xsampletests.mood` : Mood test for equal distribution scale parameters from two independent samples
+
+#### Parameter tests
+- `xsampletests.ttest_ind` : t-test for the means of two independent samples
+- `xsampletests.ttest_rel` : t-test for the means of two related samples
+- `xsampletests.bartlett` : Bartlett test for the variances of K independent samples
+- `xsampletests.levene` : Levene test for the variances of K independent samples
+- `xsampletests.fligner` : Fligner-Killeen test for the variances of K independent samples
+- `xsampletests.median_test` : Mood test for the medians of K independent samples
+
+#### Other tests
+- `xsampletests.mannwhitneyu` : Mann-Whitney U rank test on two independent samples
+- `xsampletests.ranksums` : Wilcoxon rank-sum statistic on two independent samples
+- `xsampletests.kruskal` : Kruskal-Wallis H-test on K independent samples
+- `xsampletests.friedmanchisquare` : Friedman chi-squared test on K repeated samples
+- `xsampletests.brunnermunzel` : Brunner-Munzel test on two independent samples
 
 This package was originally called `xks` but was renamed when tests additional to the KS test were added.
 


Test patch:
diff --git a/xsampletests/__init__.py b/xsampletests/__init__.py
index 8d66382..be5ce7d 100644
--- a/xsampletests/__init__.py
+++ b/xsampletests/__init__.py
@@ -1,4 +1,23 @@
-from .core import ks_2samp_1d, ks_2samp_2d, anderson_ksamp
+from .core import (
+    ks_2samp_1d,
+    ks_2samp_2d,
+    anderson_ksamp,
+    ttest_ind,
+    ttest_rel,
+    cramervonmises_2samp,
+    epps_singleton_2samp,
+    mannwhitneyu,
+    ranksums,
+    kruskal,
+    friedmanchisquare,
+    brunnermunzel,
+    ansari,
+    bartlett,
+    levene,
+    fligner,
+    median_test,
+    mood,
+)
 
 from . import _version
 
diff --git a/xsampletests/core.py b/xsampletests/core.py
index b35e641..4ea17c0 100755
--- a/xsampletests/core.py
+++ b/xsampletests/core.py
@@ -6,21 +6,147 @@
 
 SAMPLE_DIM = "xsampletest_sample_dim"
 
-# outputs[0] -> test statistic, outputs[1] -> p-value
+# 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],
+    },
 }
 
 
@@ -64,8 +190,6 @@ def _prep_data(*args, dim, nd):
     return args_prepped, input_core_dims
 
 
-# 1-dimensional tests
-# -------------------
 def _wrap_scipy(func, args, dim, kwargs):
     """Generic xarray wrapper for subset of scipy stats functions"""
 
@@ -80,15 +204,36 @@ def _wrap_scipy_func(*args, scipy_func_info, scipy_kwargs):
             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 xsampletests")
+        if kw == "axis":
+            vectorize = False
+            scipy_kwargs["axis"] = -1
+
     output_core_dims = [[]] * 2
     output_dtypes = ["float32"] * 2
-    kwargs = dict(scipy_func_info=scipy_function_info[func], scipy_kwargs=kwargs)
+    kwargs = dict(scipy_func_info=info, scipy_kwargs=scipy_kwargs)
     statistic, pvalue = xr.apply_ufunc(
         _wrap_scipy_func,
         *args,
@@ -96,7 +241,7 @@ def _wrap_scipy_func(*args, scipy_func_info, scipy_kwargs):
         input_core_dims=input_core_dims,
         output_core_dims=output_core_dims,
         output_dtypes=output_dtypes,
-        vectorize=True,
+        vectorize=vectorize,
         dask="parallelized",
     )
 
@@ -105,15 +250,18 @@ def _wrap_scipy_func(*args, scipy_func_info, scipy_kwargs):
 
 def ks_2samp_1d(ds1, ds2, dim, kwargs={}):
     """
-    One-dimensional Kolmogorov-Smirnov test on two samples. This test compares the
-    underlying continuous distributions F(x) and G(x) of two independent samples.
+    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
+        Sample 1 data. Nans are automatically removed prior to executing the test
     ds2 : xarray Dataset
-        Sample 2 data. Size of two samples can be different
+        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
@@ -124,7 +272,13 @@ def ks_2samp_1d(ds1, ds2, dim, kwargs={}):
     statistics : xarray Dataset
         Dataset with the following variables:
         - "statistic" : The KS statistic
-        - "pvalue" : One-tailed or two-tailed p-value.
+        - "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
     --------
@@ -141,7 +295,7 @@ def ks_2samp_1d(ds1, ds2, dim, kwargs={}):
 
 def anderson_ksamp(*args, dim, kwargs={}):
     """
-    The Anderson-Darling test for k-samples.
+    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
@@ -151,7 +305,8 @@ def anderson_ksamp(*args, dim, kwargs={}):
     Parameters
     ----------
     args : xarray Datasets
-        The k samples of data. Sizes of samples can be different
+        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
@@ -161,21 +316,563 @@ def anderson_ksamp(*args, dim, kwargs={}):
     -------
     statistics : xarray Dataset
         Dataset with the following variables:
-        - "statistic" : Normalized k-sample Anderson-Darling test statistic.
+        - "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%.
+            at 0.1% / 25%
 
-    See also
-    --------
-    scipy.stats.anderson_ksamp
+    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)
 
 
-# 2-dimensional tests
-# -------------------
+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 ks_2samp_2d_np(x1, y1, x2, y2):
     """
     Two-dimensional Kolmogorov-Smirnov test on two samples. For now, returns only
@@ -265,7 +962,7 @@ def _maxdist(x1, y1, x2, y2):
     ]
     assert (x1.shape[-1] == y1.shape[-1]) and (x2.shape[-1] == y2.shape[-1])
     assert (x1.shape[:-1] == x2.shape[:-1]) and (y1.shape[:-1] == y2.shape[:-1])
-    # n1, n2 = x1.shape[-1], x2.shape[-1]
+
     D1 = _maxdist(x1, y1, x2, y2)
     D2 = _maxdist(x2, y2, x1, y1)
     return (D1 + D2) / 2
@@ -276,19 +973,23 @@ def ks_2samp_2d(ds1, ds2, dim):
     Two-dimensional Kolmogorov-Smirnov test on two samples. For now, returns only
     the KS statistic with the expectation that confidence is assigned via resampling.
 
+    This test compares the underlying continuous distributions ds1 and ds2 of two
+    two-dimensional independent samples.
+
     Parameters
     ----------
     ds1 : xarray Dataset
         Sample 1 data containing two variables
     ds2 : xarray Dataset
-        Sample 2 data containing two variables
+        Sample 2 data containing two variables. The sizes of samples 1 and 2 along dim
+        can be different
     dim : str
-        The name of the sample dimension in ds1 and ds2
+        The name of the sample dimension(s) in ds1 and ds2
 
     Returns
     -------
     statistic : xarray Dataset
-        KS statistic estimating the max difference between the join distributions
+        KS statistic estimating the max difference between the joint distributions
 
     References
     ----------
@@ -305,6 +1006,9 @@ def ks_2samp_2d(ds1, ds2, dim):
         ks_2samp_2d_np,
         *args,
         input_core_dims=input_core_dims,
+        output_core_dims=[[]],
+        output_dtypes=["float32"],
+        dask="parallelized",
     )
 
     return res.to_dataset(name="statistic")
diff --git a/xsampletests/tests/fixtures.py b/xsampletests/tests/fixtures.py
index 2ce81aa..40ad72d 100644
--- a/xsampletests/tests/fixtures.py
+++ b/xsampletests/tests/fixtures.py
@@ -19,7 +19,7 @@ def raise_if_computed():
     return a
 
 
-def ds_1var(shape, dask=False):
+def ds_1var(shape, add_nans=False, dask=False):
     """An example DataSet.
     The first dimension is named 'sample' and the data is identical over other dimensions
     """
@@ -32,4 +32,7 @@ def ds_1var(shape, dask=False):
             np.random.random(size=shape[0]), list(range(1, len(shape)))
         )
         data = np.tile(data, (*shape[1:],))
+    if add_nans:
+        n_nan = max(int(shape[0] / 5), 1)
+        data[np.random.choice(shape[0], n_nan, replace=False), ...] = np.nan
     return xr.DataArray(data, coords=coords, dims=dims).to_dataset(name="var")
diff --git a/xsampletests/tests/test_core.py b/xsampletests/tests/test_core.py
index 64e1985..9772357 100755
--- a/xsampletests/tests/test_core.py
+++ b/xsampletests/tests/test_core.py
@@ -12,7 +12,7 @@
 from .fixtures import ds_1var
 
 
-def check_vs_scipy_func(func, args, dask, kwargs={}):
+def check_vs_scipy_func(func, args, kwargs={}):
     """Test values relative to scipy function"""
 
     def _stack_sample_dim(ds):
@@ -55,31 +55,31 @@ def _test_vs_scipy_values(inputs, outputs, func_info, kwargs={}):
     # Test with a single sample dim
 
     outputs = getattr(xst, func)(*args, dim="sample", kwargs=kwargs)
-    if dask is False:
-        _test_vs_scipy_values(args, outputs, function_info, 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"])
-    if dask is False:
-        _test_vs_scipy_values(args, outputs_stack, function_info)
+    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("dask", [True, False])
 @pytest.mark.parametrize("alternative", ["two-sided", "less", "greater"])
 @pytest.mark.parametrize("method", ["auto", "exact", "asymp"])
-def test_ks_2samp_1d(
-    ds1_n_per_sample, ds2_n_per_sample, shape, dask, alternative, method
+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, dask),
-        ds_1var((ds2_n_per_sample,) + shape, dask),
+        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, dask, kwargs)
+    check_vs_scipy_func("ks_2samp_1d", args, kwargs)
 
 
 @pytest.mark.parametrize("k_samples", [2, 3, 5])
@@ -87,22 +87,417 @@ def test_ks_2samp_1d(
     "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("dask", [True, False])
 @pytest.mark.parametrize("midrank", [True, False])
-def test_anderson_ksamp(k_samples, n_per_sample, shape, dask, midrank):
-    args = [ds_1var((n,) + shape, dask) for n in n_per_sample[slice(k_samples)]]
+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, dask, kwargs)
+    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)])
-def test_ks_2samp_2d_identical(samples, shape):
+@pytest.mark.parametrize("dask", [True, False])
+def test_ks_2samp_2d_identical(samples, shape, dask):
     """Check that KS statistical is zero for identical arrays"""
-    ds1_v1 = ds_1var((samples,) + shape).rename({"var": "var_1"})
-    ds1_v2 = ds_1var((samples,) + shape).rename({"var": "var_2"})
+    ds1_v1 = ds_1var((samples,) + shape, dask=dask).rename({"var": "var_1"})
+    ds1_v2 = ds_1var((samples,) + shape, dask=dask).rename({"var": "var_2"})
     ds1 = xr.merge([ds1_v1, ds1_v2])
     ds2 = ds1.copy()
     D = xst.ks_2samp_2d(ds1, ds2, "sample")
 
-    npt.assert_allclose(D["statistic"].values, 0.0)
+    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")