Repository: cogent3/cogent3
Pull request #1170 (2022-10-12 03:36:35)

DEP: deprecating distribution stats, replacing usage with SciPy


Patch:
diff --git a/doc/cookbook/evo_modelling.rst b/doc/cookbook/evo_modelling.rst
index 30b00ab32..c562983be 100644
--- a/doc/cookbook/evo_modelling.rst
+++ b/doc/cookbook/evo_modelling.rst
@@ -433,11 +433,11 @@ We import the function for computing the probability of a chi-square test statis
 
 .. jupyter-execute::
 
-    from cogent3.maths.stats import chisqprob
+    from scipy.stats.distributions import chi2
 
     LR = 2 * (alt_lnL - null_lnL)  # the likelihood ratio statistic
     df = alt_nfp - null_nfp  # the test degrees of freedom
-    p = chisqprob(LR, df)
+    p = chi2.sf(LR, df)
     print(f"LR={LR:.4f} ; df={df}; p={df:.4f}")
 
 Testing Hypotheses - By parametric bootstrapping
diff --git a/doc/examples/codon_models.rst b/doc/examples/codon_models.rst
index 0138e59bd..e6fe9f02e 100644
--- a/doc/examples/codon_models.rst
+++ b/doc/examples/codon_models.rst
@@ -103,11 +103,11 @@ We then conduct a likelihood ratio test whether the MLE of omega significantly i
 
 .. jupyter-execute::
 
-    from cogent3.maths.stats import chisqprob
+    from scipy.stats.distributions import chi2
 
     LR = 2 * (non_neutral_lnL - neutral_lnL)
     df = non_neutral_nfp - neutral_nfp
-    print(chisqprob(LR, df))
+    print(chi2.sf(LR, df))
 
 Not surprisingly, this is significant. We then ask whether the Human and Chimpanzee edges have a value of omega that is significantly different from the rest of the tree.
 
@@ -125,7 +125,7 @@ Not surprisingly, this is significant. We then ask whether the Human and Chimpan
 
     LR = 2 * (chimp_human_clade_lnL - non_neutral_lnL)
     df = chimp_human_clade_nfp - non_neutral_nfp
-    print(chisqprob(LR, df))
+    print(chi2.sf(LR, df))
 
 This is basically a replication of the original Huttley et al (2000) result for *BRCA1*.
 
@@ -194,7 +194,7 @@ The above statement essentially assigns a probability of nearly 1 to the 'neutra
 
 .. jupyter-execute::
 
-    print(chisqprob(LR, df))
+    print(chi2.sf(LR, df))
 
 We can get the posterior probabilities of site-classifications out of this model as
 
diff --git a/doc/examples/hmm_par_heterogeneity.rst b/doc/examples/hmm_par_heterogeneity.rst
index bc254b1c0..c3432b215 100644
--- a/doc/examples/hmm_par_heterogeneity.rst
+++ b/doc/examples/hmm_par_heterogeneity.rst
@@ -25,7 +25,7 @@ First import standard components necessary for all of the following calculations
         TimeReversibleNucleotide,
         predicate,
     )
-    from cogent3.maths.stats import chisqprob
+    from scipy.stats.distributions import chi2
 
 Load the alignment and tree.
 
@@ -110,7 +110,7 @@ We conduct the test between the sequentially nested models.
 
 .. jupyter-execute::
 
-    print("%.4f" % chisqprob(lr, df_patches - df_bins))
+    print("%.4f" % chi2.sf(lr, df_patches - df_bins))
 
 The stationary bin probabilities are labelled as ``bprobs`` and can be obtained as follows.
 
diff --git a/doc/examples/relative_rate.rst b/doc/examples/relative_rate.rst
index 1a86f622e..49d4756c2 100644
--- a/doc/examples/relative_rate.rst
+++ b/doc/examples/relative_rate.rst
@@ -92,7 +92,7 @@ These two lnL's are now used to calculate the likelihood ratio statistic it's de
 
     LR = 2 * (lf.get_log_likelihood() - null_lnL)
     df = lf.get_num_free_params() - null_nfp
-    P = stats.chisqprob(LR, df)
+    P = stats.chi2.sf(LR, df)
 
 Print this and look up a :math:`\chi^2` with number of edges - 1 degrees of freedom.
 
diff --git a/src/cogent3/app/result.py b/src/cogent3/app/result.py
index 6377c3c30..871918d11 100644
--- a/src/cogent3/app/result.py
+++ b/src/cogent3/app/result.py
@@ -7,8 +7,9 @@ from pathlib import Path
 
 import numpy
 
+from scipy.stats.distributions import chi2
+
 from cogent3.app.data_store import get_data_source
-from cogent3.maths.stats import chisqprob
 from cogent3.util.misc import extend_docstring_from, get_object_provenance
 from cogent3.util.table import Table
 
@@ -605,13 +606,13 @@ class hypothesis_result(model_collection_result):
 
     @property
     def pvalue(self):
-        """returns p-value from chisqprob(LR, df)
+        """returns p-value from chi2.sf(LR, df)
 
         None if LR < 0"""
         if self.LR == 0:
             pvalue = 1
         elif self.LR > 0:
-            pvalue = chisqprob(self.LR, self.df)
+            pvalue = chi2.sf(self.LR, self.df)
         else:
             pvalue = None
         return pvalue
diff --git a/src/cogent3/maths/stats/__init__.py b/src/cogent3/maths/stats/__init__.py
index 1fc190e32..a12b15d4e 100644
--- a/src/cogent3/maths/stats/__init__.py
+++ b/src/cogent3/maths/stats/__init__.py
@@ -4,11 +4,9 @@
 Also provides NumberList and FrequencyDistribution, two classes for
 working with statistical data.
 """
-from .distribution import chi_high as chisqprob
 
 
 __all__ = [
-    "chisqprob",
     "contingency",
     "distribution",
     "information_criteria",
diff --git a/src/cogent3/maths/stats/contingency.py b/src/cogent3/maths/stats/contingency.py
index 92dd84318..3a7853bdd 100644
--- a/src/cogent3/maths/stats/contingency.py
+++ b/src/cogent3/maths/stats/contingency.py
@@ -2,7 +2,8 @@ from numpy import log, outer, sqrt, zeros
 from numpy.random import shuffle
 from numpy.testing import assert_allclose
 
-from cogent3.maths.stats import chisqprob
+from scipy.stats.distributions import chi2
+
 from cogent3.maths.stats.test import G_fit
 from cogent3.util.dict_array import DictArray
 
@@ -275,7 +276,7 @@ class CategoryCounts:
         """
         stat = calc_chisq(self.observed.array, self.expected.array)
         if not shuffled:
-            pval = chisqprob(stat, self.df)
+            pval = chi2.sf(stat, self.df)
         else:
             pval = estimate_pval(self.observed.array, calc_chisq, num_reps=shuffled)
         title = "Chisq-test for independence"
@@ -317,7 +318,7 @@ class CategoryCounts:
             williams=williams,
         )
         if not shuffled:
-            pval = chisqprob(G, self.df)
+            pval = chi2.sf(G, self.df)
         else:
             pval = estimate_pval(obs.array, calc_G, num_reps=shuffled)
 
diff --git a/src/cogent3/maths/stats/distribution.py b/src/cogent3/maths/stats/distribution.py
index be3e6977e..21265e5b6 100644
--- a/src/cogent3/maths/stats/distribution.py
+++ b/src/cogent3/maths/stats/distribution.py
@@ -6,6 +6,8 @@ which is (c) Stephen L. Moshier 1984, 1995.
 from numpy import arctan as atan
 from numpy import array, exp, sqrt
 
+from scipy.stats import norm, f, t
+
 from cogent3.maths.stats.special import (
     MACHEP,
     MAXNUM,
@@ -43,12 +45,17 @@ incbet = betai  # shouldn't have renamed it...
 # Probability integrals: low gives left-hand tail, high gives right-hand tail.
 
 
-def z_low(x):
+def z_low(x):  # pragma: no cover
     """Returns left-hand tail of z distribution (0 to x).
 
     x ranges from -infinity to +infinity; result ranges from 0 to 1
 
     See Cephes docs for details."""
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "z_low", "2022.12", "use scipy.stats.norm.cdf")
+
     y = x * SQRTH
     z = abs(y)  # distribution is symmetric
     if z < SQRTH:
@@ -60,12 +67,17 @@ def z_low(x):
             return 0.5 * erfc(z)
 
 
-def z_high(x):
+def z_high(x):  # pragma: no cover
     """Returns right-hand tail of z distribution (0 to x).
 
     x ranges from -infinity to +infinity; result ranges from 0 to 1
 
     See Cephes docs for details."""
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "z_high", "2022.12", "use scipy.stats.norm.sf")
+
     y = x * SQRTH
     z = abs(y)
     if z < SQRTH:
@@ -79,10 +91,10 @@ def z_high(x):
 
 def zprob(x):
     """Returns both tails of z distribution (-inf to -x, inf to x)."""
-    return 2 * z_high(abs(x))
+    return 2 * norm.sf(abs(x))
 
 
-def chi_low(x, df):
+def chi_low(x, df):  # pragma: no cover
     """Returns left-hand tail of chi-square distribution (0 to x), given df.
 
     x ranges from 0 to infinity.
@@ -94,6 +106,13 @@ def chi_low(x, df):
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued(
+        "function", "chi_low", "2022.12", "use scipy.stats.distributions.chi2.cdf"
+    )
+
     x = fix_rounding_error(x)
     if x < 0:
         raise ValueError(f"chi_low: x must be >= 0 (got {x}).")
@@ -102,7 +121,7 @@ def chi_low(x, df):
     return igam(df / 2, x / 2)
 
 
-def chi_high(x, df):
+def chi_high(x, df):  # pragma: no cover
     """Returns right-hand tail of chi-square distribution (x to infinity).
 
     df, the degrees of freedom, ranges from 1 to infinity (assume integers).
@@ -112,6 +131,13 @@ def chi_high(x, df):
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued(
+        "function", "chi_high", "2022.12", "use scipy.stats.distributions.chi2.sf"
+    )
+
     x = fix_rounding_error(x)
 
     if x < 0:
@@ -121,7 +147,7 @@ def chi_high(x, df):
     return igamc(df / 2, x / 2)
 
 
-def t_low(t, df):
+def t_low(t, df):  # pragma: no cover
     """Returns left-hand tail of Student's t distribution (-infinity to x).
 
     df, the degrees of freedom, ranges from 1 to infinity.
@@ -131,12 +157,17 @@ def t_low(t, df):
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "t_low", "2022.12", "use scipy.stats.t.cdf")
+
     if df < 1:
         raise ValueError(f"t_low: df must be >= 1 (got {df}).")
     return stdtr(df, t)
 
 
-def t_high(t, df):
+def t_high(t, df):  # pragma: no cover
     """Returns right-hand tail of Student's t distribution (x to infinity).
 
     df, the degrees of freedom, ranges from 1 to infinity.
@@ -146,14 +177,19 @@ def t_high(t, df):
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "t_high", "2022.12", "use scipy.stats.t.sf")
+
     if df < 1:
         raise ValueError(f"t_high: df must be >= 1 (got {df}).")
     return stdtr(df, -t)  # distribution is symmetric
 
 
-def tprob(t, df):
+def tprob(x, df):
     """Returns both tails of t distribution (-infinity to -x, infinity to x)"""
-    return 2 * t_high(abs(t), df)
+    return 2 * t.sf(abs(x), df)
 
 
 def poisson_high(successes, mean):
@@ -185,15 +221,25 @@ def poisson_exact(successes, mean):
         return pdtrc(successes - 1, mean) - pdtrc(successes, mean)
 
 
-def binomial_high(successes, trials, prob):
+def binomial_high(successes, trials, prob):  # pragma: no cover
     """Returns right-hand binomial tail (X > successes) given prob(success)."""
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "binomial_high", "2022.12", "use scipy.stats.binom.sf")
+
     if -1 <= successes < 0:
         return 1
     return bdtrc(successes, trials, prob)
 
 
-def binomial_low(successes, trials, prob):
+def binomial_low(successes, trials, prob):  # pragma: no cover
     """Returns left-hand binomial tail (X <= successes) given prob(success)."""
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "binomial_low", "2022.12", "use scipy.stats.binom.cdf")
+
     return bdtr(successes, trials, prob)
 
 
@@ -213,7 +259,7 @@ def binomial_exact(successes, trials, prob):
     return exp(ln_binomial(successes, trials, prob))
 
 
-def f_low(df1, df2, x):
+def f_low(df1, df2, x):  # pragma: no cover
     """Returns left-hand tail of f distribution (0 to x).
 
     x ranges from 0 to infinity.
@@ -222,16 +268,26 @@ def f_low(df1, df2, x):
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "f_low", "2022.12", "use scipy.stats.f.cdf(x, df1, df2)")
+
     return fdtr(df1, df2, x)
 
 
-def f_high(df1, df2, x):
+def f_high(df1, df2, x):  # pragma: no cover
     """Returns right-hand tail of f distribution (x to infinity).
 
     Result ranges from 0 to 1.
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "f_high", "2022.12", "use scipy.stats.f.sf(x, df1, df2)")
+
     return fdtrc(df1, df2, x)
 
 
@@ -249,9 +305,9 @@ def fprob(dfn, dfd, F, side="right"):
     if F < 0:
         raise ValueError(f"fprob: F must be >= 0 (got {F}).")
     if side == "right":
-        return 2 * f_high(dfn, dfd, F)
+        return 2 * f.sf(F, dfn, dfd)
     elif side == "left":
-        return 2 * f_low(dfn, dfd, F)
+        return 2 * f.cdf(F, dfn, dfd)
     else:
         raise ValueError(f"Not a valid value for side {side}")
 
@@ -380,11 +436,16 @@ def pdtrc(k, m):
     return igam(k + 1, m)
 
 
-def fdtr(a, b, x):
+def fdtr(a, b, x):  # pragma: no cover
     """Returns left tail of F distribution, 0 to x.
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "fdtr", "2022.12", "use scipy.stats.f.cdf(x, df1, df2)")
+
     if min(a, b) < 1:
         raise ValueError("F a and b (degrees of freedom) must both be >= 1.")
     if x < 0:
@@ -394,11 +455,16 @@ def fdtr(a, b, x):
     return betai(0.5 * a, 0.5 * b, w)
 
 
-def fdtrc(a, b, x):
+def fdtrc(a, b, x):  # pragma: no cover
     """Returns right tail of F distribution, x to infinity.
 
     See Cephes docs for details.
     """
+
+    from cogent3.util.warning import discontinued
+
+    discontinued("function", "fdtrc", "2022.12", "use scipy.stats.f.sf(x, df1, df2)")
+
     if min(a, b) < 1:
         raise ValueError("F a and b (degrees of freedom) must both be >= 1.")
     if x < 0:
@@ -424,8 +490,15 @@ def gdtrc(a, b, x):
 # note: ndtri for the normal distribution is already imported
 
 
-def chdtri(df, y):
+def chdtri(df, y):  # pragma: no cover
     """Returns inverse of chi-squared distribution."""
+
+    from cogent3.util.warning import discontinued
+
+    discontinued(
+        "function", "chdtri", "2022.12", "use scipy.stats.distributions.chi2.isf"
+    )
+
     y = fix_rounding_error(y)
     if y < 0.0 or y > 1.0 or df < 1.0:
         raise ZeroDivisionError("y must be between 0 and 1; df >= 1")
diff --git a/src/cogent3/recalculation/scope.py b/src/cogent3/recalculation/scope.py
index d3e583cb2..e122e906c 100644
--- a/src/cogent3/recalculation/scope.py
+++ b/src/cogent3/recalculation/scope.py
@@ -6,8 +6,9 @@ from contextlib import contextmanager
 
 import numpy
 
+from scipy.stats.distributions import chi2
+
 from cogent3.maths.optimisers import MaximumEvaluationsReached
-from cogent3.maths.stats.distribution import chdtri
 
 from .calculation import Calculator
 from .setting import ConstVal, Var
@@ -728,7 +729,7 @@ class ParameterController(object):
     def get_param_interval(self, par_name, *args, **kw):
         """Confidence interval for 'par_name' found by adjusting the
         single parameter until the final result falls by 'dropoff', which
-        can be specified directly or via 'p' as chdtri(1, p).  Additional
+        can be specified directly or via 'p' as chi2.isf(1, p).  Additional
         arguments are taken to specify the scope."""
         dropoff = kw.pop("dropoff", None)
         p = kw.pop("p", None)
@@ -763,7 +764,7 @@ class ParameterController(object):
         """Make a setting -> value function"""
         if p is not None:
             assert dropoff is None, (p, dropoff)
-            dropoff = chdtri(1, p) / 2.0
+            dropoff = chi2.isf(1, p) / 2.0
         if dropoff is None:
 
             def callback(defn, posn):


Test patch:
diff --git a/doc/examples/neutral_test.rst b/doc/examples/neutral_test.rst
index dd19b7431..d26c087fe 100644
--- a/doc/examples/neutral_test.rst
+++ b/doc/examples/neutral_test.rst
@@ -90,7 +90,7 @@ The lnL's from the two models are now used to calculate the likelihood ratio sta
 
     LR = 2 * (lf.get_log_likelihood() - null_lnL)
     df = lf.get_num_free_params() - null_nfp
-    P = stats.chisqprob(LR, df)
+    P = stats.chi2.sf(LR, df)
 
 Print this and look up a chi-sq with number of edges - 1 degrees of freedom.
 
diff --git a/doc/examples/testing_multi_loci.rst b/doc/examples/testing_multi_loci.rst
index 15a3cd291..da4cb8cae 100644
--- a/doc/examples/testing_multi_loci.rst
+++ b/doc/examples/testing_multi_loci.rst
@@ -14,7 +14,7 @@ We want to know whether an exchangeability parameter is different between alignm
 
     from cogent3 import load_aligned_seqs, make_table, make_tree
     from cogent3.evolve.models import HKY85
-    from cogent3.maths.stats import chisqprob
+    from scipy.stats.distributions import chi2
     from cogent3.recalculation.scope import ALL, EACH
 
     aln = load_aligned_seqs("data/long_testseqs.fasta")
@@ -62,5 +62,5 @@ Just to pretty up the result display, I'll print(a table consisting of the test
 .. jupyter-execute::
 
     make_table(
-        header=["LR", "df", "p"], rows=[[LR, df, chisqprob(LR, df)]], digits=2, space=3,
+        header=["LR", "df", "p"], rows=[[LR, df, chi2.sf(LR, df)]], digits=2, space=3,
     )
diff --git a/src/cogent3/maths/stats/test.py b/src/cogent3/maths/stats/test.py
index a1483f792..9d30e888e 100644
--- a/src/cogent3/maths/stats/test.py
+++ b/src/cogent3/maths/stats/test.py
@@ -28,19 +28,13 @@ from numpy import sum as npsum
 from numpy import take, tanh, trace, zeros
 from numpy.random import permutation, randint
 
+from scipy.stats.distributions import chi2
+from scipy.stats import norm, binom, f, t
+
 from cogent3.maths.stats.distribution import (
-    binomial_high,
-    binomial_low,
-    chi_high,
-    f_high,
-    f_low,
     fprob,
     ndtri,
-    t_high,
-    t_low,
     tprob,
-    z_high,
-    z_low,
     zprob,
 )
 from cogent3.maths.stats.kendall import kendalls_tau, pkendall
@@ -133,7 +127,7 @@ def var(x, axis=None):
     else:
         n = x.shape[axis]
     # compute the sum of squares from the mean(s)
-    sample_SS = npsum(x ** 2, axis) - npsum(x, axis) ** 2 / n
+    sample_SS = npsum(x**2, axis) - npsum(x, axis) ** 2 / n
     return sample_SS / (n - 1)
 
 
@@ -228,7 +222,7 @@ def G_2_by_2(a, b, c, d, williams=1, directional=1):
         q = 1 + ((((n / ab) + (n / cd)) - 1) * (((n / ac) + (n / bd)) - 1)) / (6 * n)
         G /= q
 
-    p = chi_high(max(G, 0), 1)
+    p = chi2.sf(max(G, 0), 1)
 
     # find which tail we were in if the test was directional
     if directional:
@@ -271,7 +265,7 @@ def G_ind(m, williams=False):
             / (6 * tot * df)
         )
         G = G / q
-    return G, chi_high(max(G, 0), df)
+    return G, chi2.sf(max(G, 0), df)
 
 
 def calc_contingency_expected(matrix):
@@ -334,7 +328,7 @@ def G_fit(obs, exp, williams=1):
         q = 1 + (k + 1) / (6 * obs.sum())
         G /= q
 
-    return G, chi_high(G, k - 1)
+    return G, chi2.sf(G, k - 1)
 
 
 def likelihoods(d_given_h, priors):
@@ -849,7 +843,7 @@ def spearman(x_items, y_items):
     if ties1 == 0 and ties2 == 0:
         n = len(rank1)
         sum_sqr = npsum([(x - y) ** 2 for x, y in zip(rank1, rank2)])
-        rho = 1 - (6 * sum_sqr / (n * (n ** 2 - 1)))
+        rho = 1 - (6 * sum_sqr / (n * (n**2 - 1)))
     else:
         avg = lambda x: npsum(x) / len(x)
 
@@ -1265,24 +1259,24 @@ def z_test(a, popmean=0, popstdev=1, tails=None):
 def z_tailed_prob(z, tails):
     """Returns appropriate p-value for given z, depending on tails."""
     if tails == "high":
-        return z_high(z)
+        return norm.sf(z)
     elif tails == "low":
-        return z_low(z)
+        return norm.cdf(z)
     else:
         return zprob(z)
 
 
-def t_tailed_prob(t, df, tails):
+def t_tailed_prob(x, df, tails):
     """Return appropriate p-value for given t and df, depending on tails."""
     tails = tails or "2"
     tails = _get_alternate(str(tails))
 
     if tails == ALT_HIGH:
-        return t_high(t, df)
+        return t.sf(x, df)
     elif tails == ALT_LOW:
-        return t_low(t, df)
+        return t.cdf(x, df)
     else:
-        return tprob(t, df)
+        return tprob(x, df)
 
 
 def reverse_tails(tails):
@@ -1356,7 +1350,7 @@ def fisher(probs):
     -2 * SUM(ln(P)) gives chi-squared distribution with 2n degrees of freedom.
     """
     try:
-        return chi_high(-2 * npsum(list(map(log, probs))), 2 * len(probs))
+        return chi2.sf(-2 * npsum(list(map(log, probs))), 2 * len(probs))
     except OverflowError:
         return 0.0
 
@@ -1395,9 +1389,9 @@ def f_two_sample(a, b, tails=None):
 
     dfn, dfd, F = f_value(a, b)
     if tails == ALT_LOW:
-        return dfn, dfd, F, f_low(dfn, dfd, F)
+        return dfn, dfd, F, f.cdf(F, dfn, dfd)
     elif tails == ALT_HIGH:
-        return dfn, dfd, F, f_high(dfn, dfd, F)
+        return dfn, dfd, F, f.sf(F, dfn, dfd)
     else:
         if var(a) >= var(b):
             side = "right"
@@ -1442,7 +1436,7 @@ def ANOVA_one_way(a):
     dfn = len(group_means) - 1
     between_MS = between_MS / dfn
     F = between_MS / within_MS
-    return dfn, dfd, F, between_MS, within_MS, group_means, f_high(dfn, dfd, F)
+    return dfn, dfd, F, between_MS, within_MS, group_means, f.sf(F, dfn, dfd)
 
 
 def MonteCarloP(value, rand_values, tail="high"):
@@ -1490,14 +1484,14 @@ def sign_test(success, trials, alt="two sided"):
     """
     alt = _get_alternate(str(alt))
     if alt == ALT_LOW:
-        p = binomial_low(success, trials, 0.5)
+        p = binom.cdf(success, trials, 0.5)
     elif alt == ALT_HIGH:
         success -= 1
-        p = binomial_high(success, trials, 0.5)
+        p = binom.sf(success, trials, 0.5)
     else:
         success = min(success, trials - success)
-        hi = 1 - binomial_high(success, trials, 0.5)
-        lo = binomial_low(success, trials, 0.5)
+        hi = 1 - binom.sf(success, trials, 0.5)
+        lo = binom.cdf(success, trials, 0.5)
         p = hi + lo
 
     return p
@@ -1560,7 +1554,7 @@ def ks_test(x, y=None, alt="two sided", exact=None, warn_for_ties=True):
         Pval = (
             1 - pkstwo(sqrt(n) * stat)
             if alt == ALT_TWO_SIDED and not ties
-            else exp(-2 * n * stat ** 2)
+            else exp(-2 * n * stat**2)
         )
     if ties and warn_for_ties:
         warnings.warn("Cannot compute correct KS probability with ties")
@@ -1651,7 +1645,7 @@ def mw_test(x, y):
         if value != prev and start is not None:
             ave_rank = _average_rank(start, index)
             num_tied = index - start + 1
-            T += num_tied ** 3 - num_tied
+            T += num_tied**3 - num_tied
             for i in range(start - 1, index):
                 combined["rank"][i] = ave_rank
             start = None
@@ -1661,7 +1655,7 @@ def mw_test(x, y):
     if start is not None:
         ave_rank = _average_rank(start, index)
         num_tied = index - start + 2
-        T += num_tied ** 3 - num_tied
+        T += num_tied**3 - num_tied
         for i in range(start - 1, index + 1):
             combined["rank"][i] = ave_rank
 
@@ -1674,7 +1668,7 @@ def mw_test(x, y):
     U2 = prod - U1
     U = max([U1, U2])
     numerator = U - prod / 2
-    denominator = sqrt((prod / (total * (total - 1))) * ((total ** 3 - total - T) / 12))
+    denominator = sqrt((prod / (total * (total - 1))) * ((total**3 - total - T) / 12))
     z = numerator / denominator
     p = zprob(z)
     return U, p
diff --git a/tests/test_maths/test_stats/test_test.py b/tests/test_maths/test_stats/test_test.py
index 4336624f4..cde0a2d27 100644
--- a/tests/test_maths/test_stats/test_test.py
+++ b/tests/test_maths/test_stats/test_test.py
@@ -15,6 +15,7 @@ from numpy import (
     reshape,
     testing,
     tril,
+    isnan,
 )
 
 from cogent3.maths.stats.number import NumberCounter
@@ -771,7 +772,7 @@ class StatTests(TestsHelper):
         assert_allclose(t_two_sample(sample, x), (1.5637254, 0.1929248))
 
         # can't do the test if both samples have single item
-        self.assertEqual(t_two_sample(x, x), (None, None))
+        self.assertTrue(isnan(t_two_sample(x, x)).all())
 
         # Test special case if t=0.
         assert_allclose(t_two_sample([2], [1, 2, 3]), (0.0, 1.0))
diff --git a/tests/test_recalculation.py b/tests/test_recalculation.py
index 9a514fcde..17bd49381 100644
--- a/tests/test_recalculation.py
+++ b/tests/test_recalculation.py
@@ -157,7 +157,7 @@ class RecalculationTest(TestCase):
         )
 
         # For likelihood functions it is more convenient to provide 'p' rather than
-        # 'dropoff', dropoff = chdtri(1, p) / 2.0.  Also in general you won't need ultra precise answers,
+        # 'dropoff', dropoff = chi2.isf(1, p) / 2.0.  Also in general you won't need ultra precise answers,
         # so don't use 'xtol=0.0', that's just to make the doctest work.
         pc.graphviz()