diff --git a/.circleci/config.yml b/.circleci/config.yml
index 0ef7b1a4..d9411042 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -2,7 +2,6 @@ version: 2.1
 
 orbs:
   python: circleci/python@1.2
-  codecov: codecov/codecov@3.1.1
 
 jobs:
   build-and-test:
@@ -25,8 +24,11 @@ jobs:
           path: test-results
       - store_artifacts:
           path: test-results
-      - codecov/upload:
-          file: "coverage.xml"
+      - run:
+          name: Upload to Codecov
+          command: |
+            pip install codecov-cli
+            codecovcli upload-coverage --file coverage.xml
 
 workflows:
   run-tests:
diff --git a/codecov.yml b/codecov.yml
index 9fb3afac..ed4c8e12 100644
--- a/codecov.yml
+++ b/codecov.yml
@@ -19,9 +19,9 @@ coverage:
 # --------------
 # which folders/files to ignore
 ignore:
-  - *tests*
-  - *__init__.py
-  - *hyppo/_utils.py
+  - "*tests*"
+  - "*__init__.py"
+  - "*hyppo/_utils.py"
 
 # Pull request comments:
 # ----------------------
diff --git a/hyppo/conditional/FCIT.py b/hyppo/conditional/FCIT.py
index c84a7db5..c5a09454 100644
--- a/hyppo/conditional/FCIT.py
+++ b/hyppo/conditional/FCIT.py
@@ -3,10 +3,12 @@
 import joblib
 import numpy as np
 from scipy.stats import ttest_1samp
+from sklearn.base import clone
 from sklearn.metrics import mean_squared_error as mse
 from sklearn.model_selection import GridSearchCV, ShuffleSplit
 from sklearn.preprocessing import StandardScaler
 from sklearn.tree import DecisionTreeRegressor
+from sklearn.utils import check_random_state
 
 from .base import ConditionalIndependenceTest, ConditionalIndependenceTestOutput
 
@@ -30,6 +32,16 @@ class FCIT(ConditionalIndependenceTest):
         Proportion of data to evaluate test stat on.
     discrete: tuple of string
         Whether :math:`X` or :math:`Y` are discrete
+    random_state: int, RandomState instance, or None
+        Controls every source of randomness used internally (the data
+        permutations, the cross-validation splits, and, when ``model``
+        doesn't already pin its own ``random_state``, the regressor
+        fit on each permutation). Pass an int for reproducible results.
+        When ``None`` (the default), results are not reproducible and,
+        because the underlying regressor then draws from NumPy's global
+        random state, the exact statistic/p-value can also drift across
+        scikit-learn versions whenever they change how many times that
+        global state gets consumed internally.
 
     Notes
     -----
@@ -62,6 +74,7 @@ def __init__(
         num_perm=8,
         prop_test=0.1,
         discrete=(False, False),
+        random_state=None,
     ):
 
         self.model = model
@@ -69,6 +82,7 @@ def __init__(
         self.num_perm = num_perm
         self.prop_test = prop_test
         self.discrete = discrete
+        self.random_state = random_state
         ConditionalIndependenceTest.__init__(self)
 
     def statistic(self, x, y, z=None):
@@ -91,16 +105,29 @@ def statistic(self, x, y, z=None):
         n_samples = x.shape[0]
         n_test = int(n_samples * self.prop_test)
 
+        rng = check_random_state(self.random_state)
+
         data_permutations = [
-            np.random.permutation(x.shape[0]) for i in range(self.num_perm)
+            rng.permutation(x.shape[0]) for i in range(self.num_perm)
         ]
+        reshuffle_seeds = rng.randint(0, np.iinfo(np.int32).max, size=self.num_perm)
 
-        clf = _cross_val(x, y, z, self.cv_grid, self.model, prop_test=self.prop_test)
+        cv_seed = rng.randint(0, np.iinfo(np.int32).max)
+        clf = _cross_val(
+            x,
+            y,
+            z,
+            self.cv_grid,
+            self.model,
+            prop_test=self.prop_test,
+            random_state=cv_seed,
+        )
         datadict = {
             "x": x,
             "y": y,
             "z": z,
             "data_permutation": data_permutations,
+            "reshuffle_seeds": reshuffle_seeds,
             "n_test": n_test,
             "reshuffle": False,
             "clf": clf,
@@ -113,12 +140,19 @@ def statistic(self, x, y, z=None):
         )
 
         if z.shape[1] == 0:
-            x_indep_y = x[np.random.permutation(n_samples)]
+            x_indep_y = x[rng.permutation(n_samples)]
         else:
             x_indep_y = np.empty([x.shape[0], 0])
 
+        cv_seed = rng.randint(0, np.iinfo(np.int32).max)
         clf = _cross_val(
-            x_indep_y, y, z, self.cv_grid, self.model, prop_test=self.prop_test
+            x_indep_y,
+            y,
+            z,
+            self.cv_grid,
+            self.model,
+            prop_test=self.prop_test,
+            random_state=cv_seed,
         )
 
         datadict["reshuffle"] = True
@@ -192,23 +226,30 @@ def test(self, x, y, z=None):
         return ConditionalIndependenceTestOutput(stat, pvalue)
 
 
-def _cross_val(x, y, z, cv_grid, model, prop_test):
+def _cross_val(x, y, z, cv_grid, model, prop_test, random_state=None):
     """
     Choose the regression hyperparameters by
     cross-validation.
     """
 
-    splitter = ShuffleSplit(n_splits=3, test_size=prop_test)
+    splitter = ShuffleSplit(n_splits=3, test_size=prop_test, random_state=random_state)
     cv = GridSearchCV(estimator=model, cv=splitter, param_grid=cv_grid, n_jobs=-1)
     cv.fit(_interleave(x, z), y)
 
-    return type(model)(**cv.best_params_)
+    best_model = clone(model)
+    best_model.set_params(**cv.best_params_)
+
+    params = best_model.get_params()
+    if params.get("random_state") is None and random_state is not None:
+        best_model.set_params(random_state=random_state)
+
+    return best_model
 
 
 def _interleave(x, z, seed=None):
     """Interleave x and z dimension-wise."""
     state = np.random.get_state()
-    np.random.seed(seed or int(time.time()))
+    np.random.seed(seed if seed is not None else int(time.time()))
     total_ids = np.random.permutation(x.shape[1] + z.shape[1])
     np.random.set_state(state)
     out = np.zeros([x.shape[0], x.shape[1] + z.shape[1]])
@@ -227,7 +268,12 @@ def _obtain_error(data_and_i):
     y = data["y"]
     z = data["z"]
     if data["reshuffle"]:
+        reshuffle_seeds = data.get("reshuffle_seeds")
+        seed = None if reshuffle_seeds is None else int(reshuffle_seeds[i])
+        state = np.random.get_state()
+        np.random.seed(seed)
         perm_ids = np.random.permutation(x.shape[0])
+        np.random.set_state(state)
     else:
         perm_ids = np.arange(x.shape[0])
     data_permutation = data["data_permutation"][i]
diff --git a/hyppo/conditional/_utils.py b/hyppo/conditional/_utils.py
index 2a9c61e8..329705eb 100644
--- a/hyppo/conditional/_utils.py
+++ b/hyppo/conditional/_utils.py
@@ -6,13 +6,14 @@
 class _CheckInputs:
     """Checks inputs for all independence tests"""
 
-    def __init__(self, x, y, z, reps=None, max_dims=None, ignore_z_var=False):
+    def __init__(self, x, y, z=None, reps=None, max_dims=None, ignore_z_var=False):
         self.x = x
         self.y = y
         self.z = z
         self.reps = reps
         self.max_dims = max_dims
         self.ignore_z_var = ignore_z_var # to allow for constant z input
+        self.is_zero_variance = False
 
     def __call__(self):
         check_ndarray_xyz(self.x, self.y, self.z)
@@ -59,8 +60,8 @@ def check_dim_xyz(self, max_dims):
                 raise ValueError(
                     f"x, y, z must have be univariate and have shape [n,{max_dims}]"
                 )
-
         self._check_nd_indeptest()
+        self._check_variance()
 
         return self.x, self.y, self.z
 
@@ -86,10 +87,9 @@ def _check_min_samples(self):
 
     def _check_variance(self):
         if np.var(self.x) == 0:
-        # or np.var(self.y) == 0 or np.var(self.z) == 0:
-            raise ValueError("Test cannot be run. Input array x has 0 variance.")
+            self.is_zero_variance = True
         if np.var(self.y) == 0:
-            raise ValueError("Test cannot be run. Input array y has 0 variance")
-        if not self.ignore_z_var:
+            self.is_zero_variance = True
+        if not self.ignore_z_var and self.z is not None:
             if np.var(self.z) == 0:
-                raise ValueError("Test cannot be run. Input array z has 0 variance")
+                self.is_zero_variance = True
diff --git a/hyppo/conditional/cdcorr.py b/hyppo/conditional/cdcorr.py
index ce71195d..2071d8c8 100644
--- a/hyppo/conditional/cdcorr.py
+++ b/hyppo/conditional/cdcorr.py
@@ -172,6 +172,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps, ignore_z_var=True)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if not self.is_distance:
             x, y = compute_dist(x, y, metric=self.compute_distance, **self.kwargs)
             z = self._compute_kde(z)
diff --git a/hyppo/conditional/kci.py b/hyppo/conditional/kci.py
index 6b76d153..fec60e70 100644
--- a/hyppo/conditional/kci.py
+++ b/hyppo/conditional/kci.py
@@ -112,6 +112,10 @@ def test(self, x, y):
 
         T = len(y)
 
+        if np.var(x) == 0 or np.var(y) == 0:
+            self.stat = 0.0
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         Kx, Ky = self.compute_kern(x, y)
         stat = self.statistic(x, y)
 
diff --git a/hyppo/conditional/pcorr.py b/hyppo/conditional/pcorr.py
index 898e4a21..09bba286 100644
--- a/hyppo/conditional/pcorr.py
+++ b/hyppo/conditional/pcorr.py
@@ -148,6 +148,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps, max_dims=1)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if auto:  # run t-stat
             stat = self.statistic(x, y, z)
 
diff --git a/hyppo/conditional/pdcorr.py b/hyppo/conditional/pdcorr.py
index c6786b9b..17331ac1 100644
--- a/hyppo/conditional/pdcorr.py
+++ b/hyppo/conditional/pdcorr.py
@@ -181,6 +181,12 @@ def test(
         check_input = _CheckInputs(x, y, z, reps=reps)
         x, y, z = check_input()
 
+        if check_input.is_zero_variance:
+            self.stat = 0.0
+            self.pvalue = 1.0
+            self.null_dist = None
+            return ConditionalIndependenceTestOutput(0.0, 1.0)
+
         if not self.is_distance:
             x, y, z = compute_dist(x, y, z, metric=self.compute_distance, **self.kwargs)
             self.is_distance = True
diff --git a/hyppo/conditional/tests/test_FCIT.py b/hyppo/conditional/tests/test_FCIT.py
index 1ec22f9d..f5960497 100644
--- a/hyppo/conditional/tests/test_FCIT.py
+++ b/hyppo/conditional/tests/test_FCIT.py
@@ -10,22 +10,20 @@ class TestFCIT:
     @pytest.mark.parametrize(
         "n, obs_stat, obs_pvalue",
         [
-            (2000, 11.677197, 3.8168e-06),
+            (2000, 15.854040, 4.815497e-07),
         ],
     )
     def test_linear_oned(self, n, obs_stat, obs_pvalue):
         np.random.seed(123456789)
         x, y = rot_ksamp("linear", n, 1, k=2)
-        np.random.seed(123456789)
-        stat, pvalue = FCIT().test(x, y)
+        stat, pvalue = FCIT(random_state=0).test(x, y)
 
         assert_almost_equal(stat, obs_stat, decimal=-1)
         assert_almost_equal(pvalue, obs_pvalue, decimal=4)
 
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
-        # 0.56139, -0.16024
-        [(1, 100000, -0.06757, 0.52599), (2, 100000, -4.59882, 0.99876)],
+        [(1, 100000, -1.339294, 0.888837), (2, 100000, -6.985438, 0.999893)],
     )
     def test_null(self, dim, n, obs_stat, obs_pvalue):
         np.random.seed(12)
@@ -47,8 +45,9 @@ def test_null(self, dim, n, obs_stat, obs_pvalue):
             ).T
         )
 
-        np.random.seed(122)
-        stat, pvalue = FCIT().test(x1.T, y1.T, z1)
+        # random_state pins every source of randomness FCIT uses internally,
+        # so this is reproducible across scikit-learn versions (gh-427).
+        stat, pvalue = FCIT(random_state=0).test(x1.T, y1.T, z1)
 
         assert_almost_equal(pvalue, obs_pvalue, decimal=4)
         assert_almost_equal(stat, obs_stat, decimal=4)
@@ -56,9 +55,8 @@ def test_null(self, dim, n, obs_stat, obs_pvalue):
     @pytest.mark.parametrize(
         "dim, n, obs_stat, obs_pvalue",
         [
-            #89.271754, 161.35165
-            (1, 100000, 89.184784, 2.91447597e-12),
-            (2, 100000, 161.35105, 4.63412957e-14),
+            (1, 100000, 100.371119, 1.283942e-12),
+            (2, 100000, 187.502936, 1.619644e-14),
         ],
     )
     def test_alternative(self, dim, n, obs_stat, obs_pvalue):
@@ -83,8 +81,7 @@ def test_alternative(self, dim, n, obs_stat, obs_pvalue):
             ).T
         )
 
-        np.random.seed(122)
-        stat, pvalue = FCIT().test(x2.T, y2.T, z2)
+        stat, pvalue = FCIT(random_state=0).test(x2.T, y2.T, z2)
 
         assert_almost_equal(pvalue, obs_pvalue, decimal=12)
         assert_almost_equal(stat, obs_stat, decimal=4)
diff --git a/hyppo/conditional/tests/test_utils.py b/hyppo/conditional/tests/test_utils.py
index ef9f01d5..450f7d58 100644
--- a/hyppo/conditional/tests/test_utils.py
+++ b/hyppo/conditional/tests/test_utils.py
@@ -64,7 +64,9 @@ def test_constant_input(self):
         y = np.arange(20).reshape(-1, 1)
         z = np.ones(20).reshape(-1, 1)
 
-        assert_raises(ValueError, _CheckInputs(x, y, z, ignore_z_var=False))
+        check_input = _CheckInputs(x, y, z, ignore_z_var=False)
+        check_input()
+        assert check_input.is_zero_variance == True
 
         try:
             _CheckInputs(x, y, z, ignore_z_var=True)
diff --git a/hyppo/discrim/_utils.py b/hyppo/discrim/_utils.py
index 15f0f8d7..b5aca947 100644
--- a/hyppo/discrim/_utils.py
+++ b/hyppo/discrim/_utils.py
@@ -20,27 +20,18 @@ def __call__(self):
                 msg = "The input matrices do not have the same number of rows."
                 raise ValueError(msg)
 
-        tmp_ = []
-        for x1 in self.x:
+        # Pre-process validations for all inputs
+        for i, x1 in enumerate(self.x):
             check_ndarray_xy(x1, self.y)
             contains_nan(x1)
             contains_nan(self.y)
             check_min_samples(x1)
+            # convert_xy_float64 converts both x1 and self.y to float64
             x1, self.y = convert_xy_float64(x1, self.y)
-            tmp_.append(self._condition_input(x1))
-
-        self.x = tmp_
+            self.x[i] = x1
 
-        if self.reps:
-            check_reps(self.reps)
-
-        return self.x, self.y
-
-    def _condition_input(self, x1):
-        """Checks whether there is only one subject and removes
-        isolates and calculate distance."""
+        # Calculate uniqueness and isolation ONLY ONCE using the original sized y
         uniques, counts = np.unique(self.y, return_counts=True)
-
         if (counts != 1).sum() <= 1:
             msg = "You have passed a vector containing only a single unique sample id."
             raise ValueError(msg)
@@ -48,8 +39,24 @@ def _condition_input(self, x1):
         if self.remove_isolates:
             idx = np.isin(self.y, uniques[counts != 1])
             self.y = self.y[idx]
+        else:
+            idx = np.ones(len(self.y), dtype=bool)
 
-            x1 = np.asarray(x1)
+        tmp_ = []
+        for x1 in self.x:
+            tmp_.append(self._condition_input(x1, idx))
+
+        self.x = tmp_
+
+        if self.reps:
+            check_reps(self.reps)
+
+        return self.x, self.y
+
+    def _condition_input(self, x1, idx):
+        """Removes isolates and calculates distance based on given indices."""
+        x1 = np.asarray(x1)
+        if self.remove_isolates:
             if not self.is_distance:
                 x1 = x1[idx]
             else:
diff --git a/hyppo/discrim/tests/test_discrim_one_samp.py b/hyppo/discrim/tests/test_discrim_one_samp.py
index 680bffd6..5786f6b2 100644
--- a/hyppo/discrim/tests/test_discrim_one_samp.py
+++ b/hyppo/discrim/tests/test_discrim_one_samp.py
@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_almost_equal, assert_raises, assert_warns
+from numpy.testing import assert_almost_equal, assert_raises
 
 from .. import DiscrimOneSample
 
@@ -32,6 +32,13 @@ def test_diff_one(self):
         assert_almost_equal(stat, obs_stat, decimal=3)
         assert_almost_equal(p, obs_p, decimal=3)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        x = np.concatenate((np.zeros((50, 2)), np.ones((50, 2))), axis=0).astype(dtype)
+        y = np.concatenate((np.zeros(50), np.ones(50)), axis=0).astype(dtype)
+        stat, p, _ = DiscrimOneSample().test(x, y, reps=10)
+        assert_almost_equal(stat, 1.0, decimal=3)
+
 
 class TestOneSampleWarn:
     """Tests errors and warnings derived from one sample test."""
diff --git a/hyppo/discrim/tests/test_discrim_two_samp.py b/hyppo/discrim/tests/test_discrim_two_samp.py
index 36cfdc2e..b9ecf5de 100644
--- a/hyppo/discrim/tests/test_discrim_two_samp.py
+++ b/hyppo/discrim/tests/test_discrim_two_samp.py
@@ -61,3 +61,18 @@ def test_no_indeptest(self):
         x2 = np.arange(3, 23)
         y = np.arange(5, 25)
         assert_raises(ValueError, DiscrimTwoSample().test, x1, x2, y, alt="abcd")
+
+    def test_symmetry_isolates(self):
+        # test #422: ensure output is symmetrical when isolates are removed
+        x1 = np.ones((20, 2), dtype=float)
+        x2 = np.concatenate([np.zeros((10, 2)), np.ones((10, 2))], axis=0)
+        y = np.concatenate([np.zeros(10), np.ones(9), [2]], axis=0) # 2 is an isolate
+
+        discrim1, discrim2, pvalue = DiscrimTwoSample(remove_isolates=True).test(x1, x2, y, workers=1, reps=10)
+        
+        y2 = np.concatenate([np.zeros(10), np.ones(9), [2]], axis=0)
+        discrim1_rev, discrim2_rev, pvalue_rev = DiscrimTwoSample(remove_isolates=True).test(x2, x1, y2, workers=1, reps=10)
+
+        assert_almost_equal(discrim1, discrim2_rev, decimal=2)
+        assert_almost_equal(discrim2, discrim1_rev, decimal=2)
+
diff --git a/hyppo/independence/dcorr.py b/hyppo/independence/dcorr.py
index a64dfee6..27b18981 100644
--- a/hyppo/independence/dcorr.py
+++ b/hyppo/independence/dcorr.py
@@ -1,7 +1,8 @@
 import numpy as np
 from numba import jit
 
-from ..tools import check_perm_blocks_dim, chi2_approx, compute_dist
+from ..tools import check_perm_blocks_dim, chi2_approx, compute_dist, convert_xy_float64
+
 from ._utils import _CheckInputs
 from .base import IndependenceTest, IndependenceTestOutput
 
@@ -138,6 +139,11 @@ def statistic(self, x, y):
         stat : float
             The computed Dcorr statistic.
         """
+        # Only convert dtype when inputs are raw data, not precomputed distance
+        # matrices. Distance matrices produced by compute_dist are already float64;
+        # converting them on every permutation call causes large unnecessary copies.
+        if not (self.is_distance or self.is_fast):
+            x, y = convert_xy_float64(x, y)
         distx = x
         disty = y
 
@@ -302,15 +308,18 @@ def _fast_1d_dcov(x, y, bias=False):  # pragma: no cover
 
     # sort inputs
     x_orig = x.ravel()
+    y_orig = y.ravel()
     x = np.sort(x_orig)
-    y = y[np.argsort(x_orig)]
+    y = y_orig[np.argsort(x_orig)]
     x = x.reshape(-1, 1)  # for numba
+    y_col = y.reshape(-1, 1)
 
     # cumulative sum
     si = _cpu_cumsum(x)
     ax = (np.arange(-(n - 2), n + 1, 2) * x.ravel()).reshape(-1, 1) + (si[-1] - 2 * si)
 
-    v = np.hstack((x, y, x * y))
+    v = np.hstack((x, y_col, x * y_col))
+
     nw = v.shape[1]
 
     idx = np.vstack((np.arange(n), np.zeros(n))).astype(np.int64).T
@@ -370,7 +379,8 @@ def _fast_1d_dcov(x, y, bias=False):  # pragma: no cover
     c4 = np.sum(iv3.T @ x)
     d = 4 * ((c1 + c2) - (c3 + c4)) - 2 * covterm
 
-    y_sorted = y[idx[n::-1, r], :]
+    y_sorted = y_col[idx[n::-1, r], :]
+
     si = _cpu_cumsum(y_sorted)
     by = np.zeros((n, 1))
     by[idx[::-1, r]] = (np.arange(-(n - 2), n + 1, 2) * y_sorted.ravel()).reshape(
diff --git a/hyppo/independence/mgc.py b/hyppo/independence/mgc.py
index 1a715da5..03c7b5f5 100644
--- a/hyppo/independence/mgc.py
+++ b/hyppo/independence/mgc.py
@@ -4,7 +4,8 @@
 import numpy as np
 from scipy.stats import multiscale_graphcorr
 
-from ..tools import compute_dist
+from ..tools import compute_dist, convert_xy_float64
+
 from ._utils import _CheckInputs
 from .base import IndependenceTest
 
@@ -148,6 +149,11 @@ def statistic(self, x, y):
         stat : float
             The computed MGC statistic.
         """
+        # Only convert dtype when inputs are raw data, not precomputed distance
+        # matrices. Distance matrices from compute_dist are already float64;
+        # converting them on every permutation call causes large unnecessary copies.
+        if not self.is_distance:
+            x, y = convert_xy_float64(x, y)
         distx = x
         disty = y
 
@@ -158,9 +164,12 @@ def statistic(self, x, y):
 
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore")
-            mgc = multiscale_graphcorr(distx, disty, compute_distance=None, reps=0)
+            try:
+                mgc = multiscale_graphcorr(distx, disty, compute_distance=None, reps=0)
+                stat = mgc.stat
+            except (IndexError, ValueError):
+                stat = 0.0
 
-        stat = mgc.stat
         self.stat = stat
 
         return stat
@@ -251,8 +260,15 @@ def test(self, x, y, reps=1000, workers=1, random_state=None):
         # scipy gives significantly faster results
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore")
-            _, _, mgc_dict = multiscale_graphcorr(x, y, compute_distance=None, reps=0)
-        mgc_dict.pop("null_dist")
+            try:
+                _, _, mgc_dict = multiscale_graphcorr(x, y, compute_distance=None, reps=0)
+            except (IndexError, ValueError):
+                mgc_dict = {
+                    "stat_mgc_map": np.zeros((x.shape[0], y.shape[0])),
+                    "opt_scale": (x.shape[0], y.shape[0]),
+                    "null_dist": []
+                }
+        mgc_dict.pop("null_dist", None)
 
         stat, pvalue = super(MGC, self).test(
             x, y, reps, workers, random_state=random_state
diff --git a/hyppo/independence/tests/test_dcorr.py b/hyppo/independence/tests/test_dcorr.py
index 9f6aea54..28779b62 100644
--- a/hyppo/independence/tests/test_dcorr.py
+++ b/hyppo/independence/tests/test_dcorr.py
@@ -1,6 +1,6 @@
 import numpy as np
 import pytest
-from numpy.testing import assert_almost_equal, assert_raises, assert_warns
+from numpy.testing import assert_almost_equal
 
 from ...tools import linear, power
 from .. import Dcorr
@@ -30,12 +30,21 @@ def test_rep(self, n):
         assert pvalue1 == pvalue2
 
     def test_dcorr_sqrt_bug(self):
-        x = np.array([1,2,3,4,5])
-        y = np.array([1,2,9,4,4])
+        x = np.array([1, 2, 3, 4, 5])
+        y = np.array([1, 2, 9, 4, 4])
         stat = Dcorr(bias=True).test(x, y, reps=0)[0]
 
         assert_almost_equal(stat, 0.762676242417, decimal=2)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(123456789)
+        x, y = linear(100, 1)
+        x_cast = (x * 100).astype(dtype)
+        y_cast = (y * 100).astype(dtype)
+        stat = Dcorr().statistic(x_cast, y_cast)
+        assert_almost_equal(stat, 1.0, decimal=2)
+
 
 class TestDcorrTypeIError:
     def test_oned(self):
diff --git a/hyppo/independence/tests/test_mgc.py b/hyppo/independence/tests/test_mgc.py
index 97410047..6d0aaf07 100644
--- a/hyppo/independence/tests/test_mgc.py
+++ b/hyppo/independence/tests/test_mgc.py
@@ -3,11 +3,10 @@
 from numpy.testing import (
     assert_almost_equal,
     assert_approx_equal,
-    assert_equal,
     assert_warns,
 )
 
-from ...tools import linear, multimodal_independence, power, spiral
+from ...tools import linear, power, spiral
 from .. import MGC
 
 
@@ -68,6 +67,15 @@ def test_rep(self, sim, obs_stat, obs_pvalue):
         assert stat1 == stat2
         assert pvalue1 == pvalue2
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(12345678)
+        x, y = linear(n=100, p=1)
+        x_cast = (x * 100).astype(dtype)
+        y_cast = (y * 100).astype(dtype)
+        stat = MGC().statistic(x_cast, y_cast)
+        assert_approx_equal(stat, 0.97, significant=1)
+
 
 class TestMGCTypeIError:
     def test_oned(self):
diff --git a/hyppo/kgof/kernel.py b/hyppo/kgof/kernel.py
index 12628459..2cb2f741 100644
--- a/hyppo/kgof/kernel.py
+++ b/hyppo/kgof/kernel.py
@@ -86,7 +86,7 @@ def pair_gradY_X(self, X, Y):
 
     @abstractmethod
     def pair_gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(X, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: n x d numpy array.
@@ -125,7 +125,7 @@ def gradY_X(self, X, Y, dim):
 
     @abstractmethod
     def gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(x, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: nx x d numpy array.
@@ -260,7 +260,7 @@ def gradXY_sum(self, X, Y):
         return G
 
     def pair_gradXY_sum(self, X, Y):
-        """
+        r"""
         Compute \sum_{i=1}^d \frac{\partial^2 k(X, Y)}{\partial x_i \partial y_i}
         evaluated at each x_i in X, and y_i in Y.
         X: n x d numpy array.
diff --git a/hyppo/ksample/mean_embedding.py b/hyppo/ksample/mean_embedding.py
index fc855208..d426fd72 100644
--- a/hyppo/ksample/mean_embedding.py
+++ b/hyppo/ksample/mean_embedding.py
@@ -166,8 +166,8 @@ def mean_embed_distance(difference, num_randfeatures):
     mu = np.mean(difference, 0)
 
     if num_randfeatures == 1:
-        stat = float(num_samples * mu**2) / float(sigma)
+        stat = float(num_samples * mu.item() ** 2 / sigma.item())
     else:
-        stat = num_samples * mu.dot(np.linalg.solve(sigma, np.transpose(mu)))
+        stat = float(num_samples * mu.dot(np.linalg.solve(sigma, np.transpose(mu))))
 
     return stat
diff --git a/hyppo/ksample/tests/test_ksamp.py b/hyppo/ksample/tests/test_ksamp.py
index 867e8c00..263b6e96 100644
--- a/hyppo/ksample/tests/test_ksamp.py
+++ b/hyppo/ksample/tests/test_ksamp.py
@@ -45,6 +45,14 @@ def test_rep(self, n, obs_stat, obs_pvalue, indep_test):
         assert stat == stat2
         assert pvalue == pvalue2
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64, np.int32, np.int64])
+    def test_dtypes(self, dtype):
+        np.random.seed(123456789)
+        inputs = rot_ksamp("linear", 100, 1, k=2)
+        inputs_cast = [mat.astype(dtype) for mat in inputs]
+        stat, _ = KSample("Dcorr").test(*inputs_cast, reps=0)
+        assert_almost_equal(stat, 0.045646974150778084, decimal=1)
+
 
 class TestKSampleErrorWarn:
     """Tests errors and warnings derived from MGC."""
diff --git a/hyppo/time_series/_utils.py b/hyppo/time_series/_utils.py
index fa8e0557..3dfe3052 100644
--- a/hyppo/time_series/_utils.py
+++ b/hyppo/time_series/_utils.py
@@ -121,6 +121,9 @@ def compute_scale_at_lag(x, y, opt_lag, compute_distance, **kwargs):
     mgc = MGC()
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore")
-        opt_scale = mgc.test(slice_distx, slice_disty, reps=0)[2]["opt_scale"]
+        try:
+            opt_scale = mgc.test(slice_distx, slice_disty, reps=0)[2]["opt_scale"]
+        except ValueError:
+            opt_scale = (slice_distx.shape[0], slice_disty.shape[0])
 
     return opt_scale
diff --git a/hyppo/tools/common.py b/hyppo/tools/common.py
index b8c37e39..dec50a6b 100644
--- a/hyppo/tools/common.py
+++ b/hyppo/tools/common.py
@@ -171,20 +171,26 @@ def check_ndarray_or_dataframe(data, col_id):
 
 def convert_xy_float64(x, y):
     """Convert x or y to np.float64 (if not already done)"""
-    # convert x and y to floats
-    x = np.asarray(x).astype(np.float64)
-    y = np.asarray(y).astype(np.float64)
-
+    x = np.asarray(x)
+    y = np.asarray(y)
+    if x.dtype != np.float64:
+        x = x.astype(np.float64)
+    if y.dtype != np.float64:
+        y = y.astype(np.float64)
     return x, y
 
 
 def convert_xyz_float64(x, y, z):
     """Convert x or y or z to np.float64 (if not already done)"""
-    # convert x and y to floats
-    x = np.asarray(x).astype(np.float64)
-    y = np.asarray(y).astype(np.float64)
-    z = np.asarray(z).astype(np.float64)
-
+    x = np.asarray(x)
+    y = np.asarray(y)
+    z = np.asarray(z)
+    if x.dtype != np.float64:
+        x = x.astype(np.float64)
+    if y.dtype != np.float64:
+        y = y.astype(np.float64)
+    if z.dtype != np.float64:
+        z = z.astype(np.float64)
     return x, y, z