diff --git a/doc/source/tutorial/stats.rst b/doc/source/tutorial/stats.rst
index 858b70e26948..6786cde57269 100644
--- a/doc/source/tutorial/stats.rst
+++ b/doc/source/tutorial/stats.rst
@@ -99,7 +99,7 @@ introspection:
     >>> dist_discrete = [d for d in dir(stats) if
     ...                  isinstance(getattr(stats, d), stats.rv_discrete)]
     >>> print('number of continuous distributions: %d' % len(dist_continu))
-    number of continuous distributions: 103
+    number of continuous distributions: 104
     >>> print('number of discrete distributions:   %d' % len(dist_discrete))
     number of discrete distributions:   19
 
diff --git a/doc/source/tutorial/stats/continuous.rst b/doc/source/tutorial/stats/continuous.rst
index 4a6e7e275485..f072627b98a3 100644
--- a/doc/source/tutorial/stats/continuous.rst
+++ b/doc/source/tutorial/stats/continuous.rst
@@ -286,6 +286,7 @@ Continuous Distributions in `scipy.stats`
    continuous_rice
    continuous_recipinvgauss
    continuous_semicircular
+   continuous_studentized_range
    continuous_t
    continuous_trapezoid
    continuous_triang
diff --git a/doc/source/tutorial/stats/continuous_studentized_range.rst b/doc/source/tutorial/stats/continuous_studentized_range.rst
new file mode 100644
index 000000000000..95a76cf9b971
--- /dev/null
+++ b/doc/source/tutorial/stats/continuous_studentized_range.rst
@@ -0,0 +1,38 @@
+.. _continuous-studentized_range:
+
+Studentized Range Distribution
+==============================
+This distribution has two shape parameters, :math:`k>1` and :math:`\nu>0`, and the support is :math:`x \geq 0`.
+
+.. math::
+   :nowrap:
+
+    \begin{eqnarray*}
+    f(x; k, \nu) = \frac{k(k-1)\nu^{\nu/2}}{\Gamma(\nu/2)2^{\nu/2-1}}
+    \int_{0}^{\infty} \int_{-\infty}^{\infty} s^{\nu} e^{-\nu s^2/2} \phi(z) \phi(sx + z)
+    [\Phi(sx + z) - \Phi(z)]^{k-2} \,dz \,ds
+    \end{eqnarray*}
+
+.. math::
+   :nowrap:
+
+    \begin{eqnarray*}
+    F(q; k, \nu) = \frac{k\nu^{\nu/2}}{\Gamma(\nu/2)2^{\nu/2-1}}
+    \int_{0}^{\infty} \int_{-\infty}^{\infty} s^{\nu-1} e^{-\nu s^2/2} \phi(z)
+    [\Phi(sq + z) - \Phi(z)]^{k-1} \,dz \,ds
+    \end{eqnarray*}
+
+Note: :math:`\phi(z)` and :math:`\Phi(z)` represent the normal PDF and normal CDF, respectively.
+
+When :math:`\nu` exceeds 100,000, the asymptopic approximation of :math:`F(x; k, \nu=\infty)` is used:
+
+.. math::
+   :nowrap:
+
+    \begin{eqnarray*}
+    F(x; k, \nu=\infty) = k \int_{-\infty}^{\infty} \phi(z)
+    [\Phi(x + z) - \Phi(z)]^{k-1} \,dz
+    \end{eqnarray*}
+
+
+Implementation: `scipy.stats.studentized_range`
diff --git a/doc/sphinxext b/doc/sphinxext
index d2baf3a4cac5..278cf213ecd3 160000
--- a/doc/sphinxext
+++ b/doc/sphinxext
@@ -1 +1 @@
-Subproject commit d2baf3a4cac5895ed85fb94636d75dea59f3ac75
+Subproject commit 278cf213ecd308d4f77ab700d1efa8f01059b9e6
diff --git a/scipy/stats/__init__.py b/scipy/stats/__init__.py
index b42e6b38096b..86a4e96157bd 100644
--- a/scipy/stats/__init__.py
+++ b/scipy/stats/__init__.py
@@ -136,6 +136,7 @@
    semicircular      -- Semicircular
    skewcauchy        -- Skew Cauchy
    skewnorm          -- Skew normal
+   studentized_range    -- Studentized Range
    t                 -- Student's T
    trapezoid         -- Trapezoidal
    triang            -- Triangular
@@ -269,6 +270,7 @@
    siegelslopes
    theilslopes
    multiscale_graphcorr
+   tukeykramer
 
 Statistical tests
 =================
diff --git a/scipy/stats/_continuous_distns.py b/scipy/stats/_continuous_distns.py
index 83c0dca94971..3ae60e2fdd6e 100644
--- a/scipy/stats/_continuous_distns.py
+++ b/scipy/stats/_continuous_distns.py
@@ -9315,6 +9315,192 @@ def _updated_ctor_param(self):
         return dct
 
 
+class studentized_range_gen(rv_continuous):
+    r"""A studentized range continuous random variable.
+
+    %(before_notes)s
+
+    See Also
+    --------
+    t: Student's t distribution
+
+    Notes
+    -----
+    The probability density function for `studentized_range` is:
+
+    .. math::
+
+         f(x; k, \nu) = \frac{k(k-1)\nu^{\nu/2}}{\Gamma(\nu/2)
+                        2^{\nu/2-1}} \int_{0}^{\infty} \int_{-\infty}^{\infty}
+                        s^{\nu} e^{-\nu s^2/2} \phi(z) \phi(sx + z)
+                        [\Phi(sx + z) - \Phi(z)]^{k-2} \,dz \,ds
+
+    for :math:`x ≥ 0`, :math:`k > 1`, and :math:`\nu > 0`.
+
+    `studentized_range` takes ``k`` and ``v`` as shape parameters.
+
+    When :math:`\nu` exceeds 100,000, an asymptotic approximation (infinite
+    degrees of freedom) is used to compute the cumulative distribution
+    function [4]_.
+
+    %(after_notes)s
+
+    References
+    ----------
+
+    .. [1] "Studentized range distribution",
+           https://en.wikipedia.org/wiki/Studentized_range_distribution
+    .. [2] Batista, Ben Dêivide, et al. "Externally Studentized Normal Midrange
+           Distribution." Ciência e Agrotecnologia, vol. 41, no. 4, 2017, pp.
+           378-389., doi:10.1590/1413-70542017414047716.
+    .. [3] Harter, H. Leon. "Tables of Range and Studentized Range." The Annals
+           of Mathematical Statistics, vol. 31, no. 4, 1960, pp. 1122-1147.
+           JSTOR, www.jstor.org/stable/2237810. Accessed 18 Feb. 2021.
+    .. [4] Lund, R. E., and J. R. Lund. "Algorithm AS 190: Probabilities and
+           Upper Quantiles for the Studentized Range." Journal of the Royal
+           Statistical Society. Series C (Applied Statistics), vol. 32, no. 2,
+           1983, pp. 204-210. JSTOR, www.jstor.org/stable/2347300. Accessed 18
+           Feb. 2021.
+
+    Examples
+    --------
+    >>> from scipy.stats import studentized_range
+    >>> import matplotlib.pyplot as plt
+    >>> fig, ax = plt.subplots(1, 1)
+
+    Calculate the first four moments:
+
+    >>> k, df = 3, 10
+    >>> mean, var, skew, kurt = studentized_range.stats(k, df, moments='mvsk')
+
+    Display the probability density function (``pdf``):
+
+    >>> x = np.linspace(studentized_range.ppf(0.01, k, df),
+    ...                 studentized_range.ppf(0.99, k, df), 100)
+    >>> ax.plot(x, studentized_range.pdf(x, k, df),
+    ...         'r-', lw=5, alpha=0.6, label='studentized_range pdf')
+
+    Alternatively, the distribution object can be called (as a function)
+    to fix the shape, location and scale parameters. This returns a "frozen"
+    RV object holding the given parameters fixed.
+
+    Freeze the distribution and display the frozen ``pdf``:
+
+    >>> rv = studentized_range(k, df)
+    >>> ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')
+
+    Check accuracy of ``cdf`` and ``ppf``:
+
+    >>> vals = studentized_range.ppf([0.001, 0.5, 0.999], k, df)
+    >>> np.allclose([0.001, 0.5, 0.999], studentized_range.cdf(vals, k, df))
+    True
+
+    Rather than using (``studentized_range.rvs``) to generate random variates,
+    which is very slow for this distribution, we can approximate the inverse
+    CDF using an interpolator, and then perform inverse transform sampling
+    with this approximate inverse CDF.
+
+    This distribution has an infinite but thin right tail, so we focus our
+    attention on the leftmost 99.9 percent.
+
+    >>> a, b = studentized_range.ppf([0, .999], k, df)
+    >>> a, b
+    0, 7.41058083802274
+
+    >>> from scipy.interpolate import interp1d
+    >>> random_state = np.random.default_rng(123)
+    >>> xs = np.linspace(a, b, 50)
+    >>> cdf = studentized_range.cdf(xs, k, df)
+    # Create an interpolant of the inverse CDF
+    >>> ppf = interp1d(cdf, xs, fill_value='extrapolate')
+    # Perform inverse transform sampling using the interpolant
+    >>> r = ppf(random_state.uniform(size=1000))
+
+    And compare the histogram:
+
+    >>> ax.hist(r, density=True, histtype='stepfilled', alpha=0.2)
+    >>> ax.legend(loc='best', frameon=False)
+    >>> plt.show()
+
+    """
+
+    def _argcheck(self, k, df):
+        return (k > 1) & (df > 0)
+
+    def _fitstart(self, data):
+        # Default is k=1, but that is not a valid value of the parameter.
+        return super(studentized_range_gen, self)._fitstart(data, args=(2, 1))
+
+    def _munp(self, K, k, df):
+        cython_symbol = '_studentized_range_moment'
+        _a, _b = self._get_support()
+        # all three of these are used to create a numpy array so they must
+        # be the same shape.
+
+        def _single_moment(K, k, df):
+            log_const = _stats._studentized_range_pdf_logconst(k, df)
+            arg = [K, k, df, log_const]
+            usr_data = np.array(arg, float).ctypes.data_as(ctypes.c_void_p)
+
+            llc = LowLevelCallable.from_cython(_stats, cython_symbol, usr_data)
+
+            ranges = [(-np.inf, np.inf), (0, np.inf), (_a, _b)]
+            opts = dict(epsabs=1e-11, epsrel=1e-12)
+
+            return integrate.nquad(llc, ranges=ranges, opts=opts)[0]
+
+        ufunc = np.frompyfunc(_single_moment, 3, 1)
+        return np.float64(ufunc(K, k, df))
+
+    def _pdf(self, x, k, df):
+        cython_symbol = '_studentized_range_pdf'
+
+        def _single_pdf(q, k, df):
+            log_const = _stats._studentized_range_pdf_logconst(k, df)
+            arg = [q, k, df, log_const]
+            usr_data = np.array(arg, float).ctypes.data_as(ctypes.c_void_p)
+
+            llc = LowLevelCallable.from_cython(_stats, cython_symbol, usr_data)
+
+            ranges = [(-np.inf, np.inf), (0, np.inf)]
+            opts = dict(epsabs=1e-11, epsrel=1e-12)
+
+            return integrate.nquad(llc, ranges=ranges, opts=opts)[0]
+
+        ufunc = np.frompyfunc(_single_pdf, 3, 1)
+        return np.float64(ufunc(x, k, df))
+
+    def _cdf(self, x, k, df):
+
+        def _single_cdf(q, k, df):
+            # "When the degrees of freedom V are infinite the probability
+            # integral takes [on a] simpler form," and a single asymptotic
+            # integral is evaluated rather than the standard double integral.
+            # (Lund, Lund, page 205)
+            if df < 100000:
+                cython_symbol = '_studentized_range_cdf'
+                log_const = _stats._studentized_range_cdf_logconst(k, df)
+                arg = [q, k, df, log_const]
+                usr_data = np.array(arg, float).ctypes.data_as(ctypes.c_void_p)
+                ranges = [(-np.inf, np.inf), (0, np.inf)]
+            else:
+                cython_symbol = '_studentized_range_cdf_asymptotic'
+                arg = [q, k]
+                usr_data = np.array(arg, float).ctypes.data_as(ctypes.c_void_p)
+                ranges = [(-np.inf, np.inf)]
+
+            llc = LowLevelCallable.from_cython(_stats, cython_symbol, usr_data)
+            opts = dict(epsabs=1e-11, epsrel=1e-12)
+            return integrate.nquad(llc, ranges=ranges, opts=opts)[0]
+
+        ufunc = np.frompyfunc(_single_cdf, 3, 1)
+        return np.float64(ufunc(x, k, df))
+
+
+studentized_range = studentized_range_gen(name='studentized_range', a=0,
+                                          b=np.inf)
+
+
 # Collect names of classes and objects in this module.
 pairs = list(globals().copy().items())
 _distn_names, _distn_gen_names = get_distribution_names(pairs, rv_continuous)
diff --git a/scipy/stats/_distr_params.py b/scipy/stats/_distr_params.py
index cddde758c2c8..437be86fcb65 100644
--- a/scipy/stats/_distr_params.py
+++ b/scipy/stats/_distr_params.py
@@ -98,6 +98,7 @@
     ['semicircular', ()],
     ['skewcauchy', (0.5,)],
     ['skewnorm', (4.0,)],
+    ['studentized_range', (3.0, 10.0)],
     ['t', (2.7433514990818093,)],
     ['trapezoid', (0.2, 0.8)],
     ['triang', (0.15785029824528218,)],
@@ -254,6 +255,7 @@
     ['recipinvgauss', (-1, )],
     ['semicircular', ()],
     ['skewnorm', (np.inf, )],
+    ['studentized_range', (-1, 1)],
     ['t', (-1, )],
     ['trapezoid', (0, 2)],
     ['triang', (2, )],
diff --git a/scipy/stats/_hypotests.py b/scipy/stats/_hypotests.py
index fd5485212c9c..64d6c65238ea 100644
--- a/scipy/stats/_hypotests.py
+++ b/scipy/stats/_hypotests.py
@@ -6,12 +6,13 @@
 import scipy.stats
 from scipy.optimize import shgo
 from . import distributions
+from ._common import ConfidenceInterval
 from ._continuous_distns import chi2, norm
 from scipy.special import gamma, kv, gammaln
 from . import _wilcoxon_data
 
 __all__ = ['epps_singleton_2samp', 'cramervonmises', 'somersd',
-           'barnard_exact', 'cramervonmises_2samp']
+           'barnard_exact', 'cramervonmises_2samp', 'tukeykramer']
 
 Epps_Singleton_2sampResult = namedtuple('Epps_Singleton_2sampResult',
                                         ('statistic', 'pvalue'))
@@ -1216,3 +1217,231 @@ def cramervonmises_2samp(x, y, method='auto'):
             p = max(0, 1. - _cdf_cvm_inf(tn))
 
     return CramerVonMisesResult(statistic=t, pvalue=p)
+
+
+TukeyKramerResult = make_dataclass("TukeyKramerResult",
+                                   ("statistic", "pvalue", "lower_conf",
+                                    "upper_conf", "significant"))
+
+
+class TukeyKramerResult:
+    """
+    Result of `scipy.stats.tukeykramer`.
+
+    Attributes
+    ----------
+    sig_level : float
+            The desired significance level (copied from `tukeykramer` input).
+    statistic : float ndarray
+        The computed statistic of the test for each comparison. The
+        ``(i, j)`` index is the p-value for the ith - jth group comparison.
+    pvalue : float ndarray
+        The associated p-value from the studentized range distribution. The
+        ``(i, j)`` index is the p-value for the ith - jth group comparison.
+
+    Methods
+    -------
+    confidence_inteval() :
+        Compute the confidence interval for the significance level
+        specified in the test.
+    """
+
+    def __init__(self, sig_level, statistic, pvalue, _nargs, _nsamples,
+                 _stand_err):
+        self.sig_level = sig_level
+        self.statistic = statistic
+        self.pvalue = pvalue
+        self._nargs = _nargs
+        self._nsamples = _nsamples
+        self._stand_err = _stand_err
+
+    def confidence_interval(self):
+        # determine the critical value of the studentized range using the
+        # appropriate significance level, number of treatments, and degrees
+        # of freedom as determined by the number of data less the number of
+        # treatments. ("Confidence limits for Tukey's method")[1]. Note that
+        # in the cases of unequal sample sizes there will be a criterion for
+        # each group comparison
+        params = (
+        1 - self.sig_level, self._nargs, self._nsamples - self._nargs)
+        srd = distributions.studentized_range.ppf(*params)
+        # also called maximum critical value, the tukey criterion is the
+        # studentized range critical value * the square root of mean square
+        # error over the sample size.
+        tukey_criterion = srd * self._stand_err
+        # the confidence levels are determined by the
+        # `mean_differences` +- `tukey_criterion`
+        upper_conf = self.statistic + tukey_criterion
+        lower_conf = self.statistic - tukey_criterion
+        return ConfidenceInterval(low=lower_conf, high=upper_conf)
+
+
+def _tukeykramer_iv(args, sig_level):
+    args = [np.asarray(arg) for arg in args]
+    for arg in args:
+        if arg.ndim != 1:
+            raise ValueError("Input samples must be one-dimensional")
+        if arg.size <= 1:
+            raise ValueError("Input sample size must be greater than one.")
+        if np.isinf(arg).any():
+            raise ValueError("Input samples must be finite.")
+    if not 0 < sig_level < 1:
+        raise ValueError("Significance level should be between 0 and 1.")
+    return args
+
+
+def tukeykramer(*args, sig_level=.05):
+    """
+    Perform Tukey-Kramer.
+
+    Under the assumption of a rejected null hypothesis that two or more samples
+    have the same population mean, the Tukey Kramer test compares the absolute
+    mean difference between each input group to the Tukey criterion for
+    significance. For inputs of differing sample sizes, the Tukey Kramer
+    method is used, albeit with a higher confidence coefficient. [1]_
+
+    Parameters
+    ----------
+    sample1, sample2, ... : array_like
+        The sample measurements for each group. There must be at least
+        two arguments.
+    sig_level : float, optional
+        Significance level for which to determine the appropriate studentized
+        range distribution value.
+        Default is .05.
+
+    Returns
+    -------
+    result : `~scipy.stats._result_classes.TukeyKramerTestResult` instance
+        The return value is an object with the following attributes:
+
+        sig_level : float
+            The desired significance level (copied from `tukeykramer` input).
+        statistic : float ndarray
+            The computed statistic of the test for each comparison. The
+            (i, j) index is the p-value for the ith - jth group comparison.
+        pvalue : float ndarray
+            The associated p-value from the studentized range distribution. The
+            (i, j) index is the p-value for the ith - jth group comparison.
+
+
+        The object has the following methods:
+
+        confidence_ci() :
+            Compute the confidence interval for the significance level
+            specified in the test.
+    Notes
+    -----
+    The use of this test relies on several assumptions.
+
+    1. There are equal variances between the samples.
+    2. Each sample is from a normally distributed population.
+
+    References
+    ----------
+    .. [1] NIST/SEMATECH e-Handbook of Statistical Methods, "7.4.7.1. Tukey's
+           Method."
+           https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm,
+           28 November 2020.
+    .. [2] Abdi, Herve & Williams, Lynne. (2021). Tukey's Honestly Signiflcant
+           Difierence (HSD) Test. 
+    .. [3] "One-Way ANOVA Using SAS PROC ANOVA & PROC GLM." SAS
+           Tutorials, 2007, www.stattutorials.com/SAS/TUTORIAL-PROC-GLM.htm. 
+    .. [4] Kramer, Clyde Young. "Extension of Multiple Range Tests to Group
+           Means with Unequal Numbers of Replications." Biometrics, vol. 12,
+           no. 3, 1956, pp. 307-310. JSTOR, www.jstor.org/stable/3001469.
+           Accessed 25 May 2021.       
+    .. [5] 
+
+
+    Examples
+    --------
+    Here are some data comparing the time to relief of three brands of
+    headache medicine, reported in minutes. Data adapted from [3]_.
+
+    >>> from scipy.stats import tukeykramer
+    >>> group0 = [24.5, 23.5, 26.4, 27.1, 29.9]
+    >>> group1 = [28.4, 34.2, 29.5, 32.2, 30.1]
+    >>> group2 = [26.1, 28.3, 24.3, 26.2, 27.8]
+
+    We would like to see if the means between any of the groups are
+    significantly different. First, visually examine a box and whisker plot.
+
+    >>> import matplotlib.pyplot as plt
+    >>> fig, ax = plt.subplots(1, 1)
+    >>> ax.boxplot([group0, group1, group2])
+    >>> ax.set_xticklabels(["group0", "group1", "group2"])
+    >>> ax.set_ylabel("mean")
+    >>> plt.show()
+
+    From a box and whisker plot we can see overlap in the interquartile ranges
+    group 1 to group 2 and group 3, but we can apply the ``tukeykramer`` test
+    to determine if the difference between means is significant.
+
+    >>> res = tukeykramer(group0, group1, group2)
+    >>> for ((i, j), p) in np.ndenumerate(res.pvalue):
+    >>>     # filter out self comparisons
+    >>>     if i != j:
+    >>>         print(f"({i} - {j}) {p:.03f}")
+    (0 - 1) 0.012
+    (0 - 2) 0.980
+    (1 - 0) 0.012
+    (1 - 2) 0.017
+    (2 - 0) 0.980
+    (2 - 1) 0.017
+
+    The null hypothesis is that each group has the same mean. We choose a
+    significance level of .05. The p-value for comparisons (``group0`` -
+    ``group1``), and (``group1`` - ``group2``), do not exceed .05, so we reject
+    the null hypothesis that they have the same means. The comparison of
+    (``group0`` - ``group2``) exceeds .05, so we accept the null hypothesis
+    that they have the same mean.
+
+    There are also confidence intervals associated with the test statistic.
+    They are applicable for the 95% confidence level becuase the test was
+    applied with a .05 significance level.
+
+    >>> conf = res.confidence_interval()
+    >>> for ((i, j), l) in np.ndenumerate(conf.low):
+    >>>     # filter out self comparisons
+    >>>     if i != j:
+    >>>         print(f"({i} - {j}) {p:.03f}")
+    """
+    args = _tukeykramer_iv(args, sig_level)
+    nargs = len(args)
+    means = np.asarray([np.mean(arg) for arg in args])
+    nsamples_args = np.asarray([a.size for a in args])
+    nsamples = np.sum(nsamples_args)
+
+    # determine mean square error
+    mse = (np.sum([np.var(arg, ddof=1) for arg in args] *
+                  (nsamples_args - 1)) / (nsamples - len(args)))
+
+    # The calculation of the standard error differs when treatments differ in
+    # size. See ("Unequal sample sizes")[1].
+    if np.unique(nsamples_args).size == 1:
+        # all input groups are the same length, so only one value needs to be
+        # calculated [1].
+        normalize = 2 / nsamples_args[0]
+    else:
+        # to compare groups of differing sizes, we must compute a variance
+        # value for each individual comparison. Use broadcasting to get the
+        # resulting matrix. [3], verified against [4] (page 308).
+        normalize = 1 / nsamples_args + 1 / nsamples_args.reshape((nargs, 1))
+
+    # the standard error is used in the computation of the tukey criterion and
+    # finding the p-values.
+    stand_err = np.sqrt(normalize * mse / 2)
+
+    # the mean difference is the test statistic.
+    mean_differences = means.reshape((nargs, 1)) - means
+
+    # Calculate the t-statistic to use within the survival function of the
+    # studentized range to get the p-value.
+    t_stat = np.abs(mean_differences) / stand_err
+
+    params = t_stat, nargs, nsamples - 1
+    pvalues = distributions.studentized_range.sf(*params)
+
+    return TukeyKramerResult(sig_level, mean_differences, pvalues, nargs,
+                             nsamples, stand_err)
\ No newline at end of file
diff --git a/scipy/stats/_stats.pxd b/scipy/stats/_stats.pxd
index 10c3de6d075a..bd7e6879e6b4 100644
--- a/scipy/stats/_stats.pxd
+++ b/scipy/stats/_stats.pxd
@@ -1,4 +1,8 @@
 # destined to be used in a LowLevelCallable
 cdef double _geninvgauss_pdf(double x, void *user_data) nogil except *
+cdef double _studentized_range_cdf(int n, double[2] x, void *user_data) nogil
+cdef double _studentized_range_cdf_asymptotic(double z, void *user_data) nogil
+cdef double _studentized_range_pdf(int n, double[2] x, void *user_data) nogil
+cdef double _studentized_range_moment(int n, double[3] x_arg, void *user_data) nogil
 cdef double _genhyperbolic_pdf(double x, void *user_data) nogil except *
 cdef double _genhyperbolic_logpdf(double x, void *user_data) nogil except *
diff --git a/scipy/stats/_stats.pyx b/scipy/stats/_stats.pyx
index 8697ea172202..528531e7a22d 100644
--- a/scipy/stats/_stats.pyx
+++ b/scipy/stats/_stats.pyx
@@ -513,6 +513,116 @@ cdef double _geninvgauss_pdf(double x, void *user_data) nogil except *:
     return math.exp(_geninvgauss_logpdf_kernel(x, p, b))
 
 
+cdef double _phi(double z) nogil:
+    """evaluates the normal PDF. Used in `studentized_range`"""
+    cdef double inv_sqrt_2pi = 0.3989422804014327
+    return inv_sqrt_2pi * math.exp(-0.5 * z * z)
+
+
+cdef double _logphi(double z) nogil:
+    """evaluates the log of the normal PDF. Used in `studentized_range`"""
+    cdef double log_inv_sqrt_2pi = -0.9189385332046727
+    return log_inv_sqrt_2pi - 0.5 * z * z
+
+
+cdef double _Phi(double z) nogil:
+    """evaluates the normal CDF. Used in `studentized_range`"""
+    # use a custom function because using cs.ndtr results in incorrect PDF at
+    # q=0 on 32bit systems. Use a hardcoded 1/sqrt(2) constant rather than
+    # math constants because they're not availible on all systems.
+    cdef double inv_sqrt_2 = 0.7071067811865475
+    return 0.5 * math.erfc(-z * inv_sqrt_2)
+
+
+cpdef double _studentized_range_cdf_logconst(double k, double df):
+    """Evaluates log of constant terms in the cdf integrand"""
+    cdef double log_2 = 0.6931471805599453
+    return (math.log(k) + (df / 2) * math.log(df)
+            - (math.lgamma(df / 2) + (df / 2 - 1) * log_2))
+
+
+cpdef double _studentized_range_pdf_logconst(double k, double df):
+    """Evaluates log of constant terms in the pdf integrand"""
+    cdef double log_2 = 0.6931471805599453
+    return (math.log(k) + math.log(k - 1) + (df / 2) * math.log(df)
+            - (math.lgamma(df / 2) + (df / 2 - 1) * log_2))
+
+
+cdef double _studentized_range_cdf(int n, double[2] integration_var,
+                                   void *user_data) nogil:
+    # evaluates the integrand of Equation (3) by Batista, et al [2]
+    # destined to be used in a LowLevelCallable
+    q = (<double *> user_data)[0]
+    k = (<double *> user_data)[1]
+    df = (<double *> user_data)[2]
+    log_cdf_const = (<double *> user_data)[3]
+
+    s = integration_var[1]
+    z = integration_var[0]
+
+    # suitable terms are evaluated within logarithms to avoid under/overflows
+    log_terms = (log_cdf_const
+                 + (df - 1) * math.log(s)
+                 - (df * s * s / 2)
+                 + _logphi(z))
+
+    # multiply remaining term outside of log because it can be 0
+    return math.exp(log_terms) * math.pow(_Phi(z + q * s) - _Phi(z), k - 1)
+
+
+cdef double _studentized_range_cdf_asymptotic(double z, void *user_data) nogil:
+    # evaluates the integrand of equation (2) by Lund, Lund, page 205. [4]
+    # destined to be used in a LowLevelCallable
+    q = (<double *> user_data)[0]
+    k = (<double *> user_data)[1]
+
+    return k * _phi(z) * math.pow(_Phi(z + q) - _Phi(z), k - 1)
+
+
+cdef double _studentized_range_pdf(int n, double[2] integration_var,
+                                   void *user_data) nogil:
+    # evaluates the integrand of equation (4) by Batista, et al [2]
+    # destined to be used in a LowLevelCallable
+    q = (<double *> user_data)[0]
+    k = (<double *> user_data)[1]
+    df = (<double *> user_data)[2]
+    log_pdf_const = (<double *> user_data)[3]
+
+    z = integration_var[0]
+    s = integration_var[1]
+
+    # suitable terms are evaluated within logarithms to avoid under/overflows
+    log_terms = (log_pdf_const
+                 + df * math.log(s)
+                 - df * s * s / 2
+                 + _logphi(z)
+                 + _logphi(s * q + z))
+
+    # multiply remaining term outside of log because it can be 0
+    return math.exp(log_terms) * math.pow(_Phi(s * q + z) - _Phi(z), k - 2)
+
+
+cdef double _studentized_range_moment(int n, double[3] integration_var,
+                                      void *user_data) nogil:
+    # destined to be used in a LowLevelCallable
+    K = (<double *> user_data)[0]  # the Kth moment to calc.
+    k = (<double *> user_data)[1]
+    df = (<double *> user_data)[2]
+    log_pdf_const = (<double *> user_data)[3]
+
+    # Pull outermost integration variable out to pass as q to PDF
+    q = integration_var[2]
+
+    cdef double pdf_data[4]
+    pdf_data[0] = q
+    pdf_data[1] = k
+    pdf_data[2] = df
+    pdf_data[3] = log_pdf_const
+
+    return (math.pow(q, K) *
+            _studentized_range_pdf(4, integration_var, pdf_data))
+
+
 cpdef double genhyperbolic_pdf(double x, double p, double a, double b) nogil:
     return math.exp(_genhyperbolic_logpdf_kernel(x, p, a, b))
 
diff --git a/scipy/stats/statlib.pyf b/scipy/stats/statlib.pyf
index eb539783c73d..cfa4e5e60da5 100644
--- a/scipy/stats/statlib.pyf
+++ b/scipy/stats/statlib.pyf
@@ -29,8 +29,8 @@ python module statlib ! in
            integer intent(in) :: is
            integer intent(out) :: ifault
            double precision intent(out) :: prho
-         end function prho         
-         
+         end function prho
+
     end interface 
 end python module statlib
 
diff --git a/scipy/stats/stats.py b/scipy/stats/stats.py
index ae611dd051ff..62439e510a65 100644
--- a/scipy/stats/stats.py
+++ b/scipy/stats/stats.py
@@ -48,6 +48,7 @@
 from ._stats import (_kendall_dis, _toint64, _weightedrankedtau,
                      _local_correlations)
 from dataclasses import make_dataclass
+from itertools import permutations
 
 
 # Functions/classes in other files should be added in `__init__.py`, not here
@@ -3691,6 +3692,7 @@ def f_oneway(*args, axis=0):
     return F_onewayResult(f, prob)
 
 
+    
 def alexandergovern(*args, nan_policy='propagate'):
     """Performs the Alexander Govern test.
 
diff --git a/scipy/stats/tests/data/studentized_range_mpmath_ref.json b/scipy/stats/tests/data/studentized_range_mpmath_ref.json
new file mode 100644
index 000000000000..ecfa9ab4fcae
--- /dev/null
+++ b/scipy/stats/tests/data/studentized_range_mpmath_ref.json
@@ -0,0 +1,329 @@
+{
+  "COMMENT": "!!!!!!!!!! THIS FILE WAS AUTOGENERATED BY RUNNING `python gen_studentized_range_mpmath_ref.py` !!!!!!!!!!",
+  "cdf_data": [
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.0027520560662338336
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.23510918942128056
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.999908236635449
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 2.8544145010066327e-12
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.001300816146573152
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.9991738325963548
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.0027527430186246545
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.2401356460422021
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.9999999993640496
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 9.39089126131273e-13
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.0005841098057517027
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.9999999905199205
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 1000,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.2406547688726472
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 1000,
+        "expected_atol": 1e-10,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.999999999991348
+    }
+  ],
+  "pdf_data": [
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.05494947290360002
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.398772020275752
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 7.728665739003835e-05
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 2.564099684606509e-10
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.009782573637596617
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 10,
+        "v": 10,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.0006784051704217357
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.05497690690332341
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.4157583991350054
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 2.159522630228393e-09
+    },
+    {
+      "src_case": {
+        "q": 0.1,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 8.442593793786411e-11
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.0047089182958790715
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 10,
+        "v": 100,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 3.215093171597309e-08
+    },
+    {
+      "src_case": {
+        "q": 1,
+        "k": 3,
+        "v": 1000,
+        "expected_atol": 1e-11,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 0.41754151596262973
+    },
+    {
+      "src_case": {
+        "q": 10,
+        "k": 3,
+        "v": 1000,
+        "expected_atol": 1e-10,
+        "expected_rtol": 1e-11
+      },
+      "mp_result": 1.4993650577943627e-29
+    }
+  ],
+  "moment_data": [
+    {
+      "src_case": {
+        "m": 1,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-9,
+        "expected_rtol": 1e-9
+      },
+      "mp_result": 1.8342745127927962
+    },
+    {
+      "src_case": {
+        "m": 2,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-9,
+        "expected_rtol": 1e-9
+      },
+      "mp_result": 4.567483357831711
+    },
+    {
+      "src_case": {
+        "m": 3,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-9,
+        "expected_rtol": 1e-9
+      },
+      "mp_result": 14.412156886227011
+    },
+    {
+      "src_case": {
+        "m": 4,
+        "k": 3,
+        "v": 10,
+        "expected_atol": 1e-9,
+        "expected_rtol": 1e-9
+      },
+      "mp_result": 56.012250366720444
+    }
+  ]
+}
\ No newline at end of file
diff --git a/scipy/stats/tests/studentized_range_mpmath_ref.py b/scipy/stats/tests/studentized_range_mpmath_ref.py
new file mode 100644
index 000000000000..a391e89d9d86
--- /dev/null
+++ b/scipy/stats/tests/studentized_range_mpmath_ref.py
@@ -0,0 +1,254 @@
+# To run this script, run
+# `python gen_studentized_range_mpmath_ref.py`
+# in the "scipy/stats/tests/" directory
+
+# This script generates a JSON file "./data/studentized_range_mpmath_ref.json"
+# that is used to compare the accuracy of `studentized_range` functions against
+# precise (20 DOP) results generated using `mpmath`.
+
+# Equations in this file have been taken from
+# https://en.wikipedia.org/wiki/Studentized_range_distribution
+# and have been checked against the following reference:
+# Lund, R. E., and J. R. Lund. "Algorithm AS 190: Probabilities and
+# Upper Quantiles for the Studentized Range." Journal of the Royal
+# Statistical Society. Series C (Applied Statistics), vol. 32, no. 2,
+# 1983, pp. 204-210. JSTOR, www.jstor.org/stable/2347300. Accessed 18
+# Feb. 2021.
+
+# Note: I would have prefered to use pickle rather than JSON, but -
+# due to security concerns - decided against it.
+import itertools
+from collections import namedtuple
+import json
+import time
+
+import os
+from multiprocessing import Pool, cpu_count
+
+import numpy as np
+from mpmath import gamma, pi, sqrt, quad, inf, mpf, mp
+from mpmath import npdf as phi
+from mpmath import ncdf as Phi
+
+results_filepath = "data/studentized_range_mpmath_ref.json"
+num_pools = max(cpu_count() - 1, 1)
+
+MPResult = namedtuple("MPResult", ["src_case", "mp_result"])
+
+CdfCase = namedtuple("CdfCase",
+                     ["q", "k", "v", "expected_atol", "expected_rtol"])
+
+MomentCase = namedtuple("MomentCase",
+                        ["m", "k", "v", "expected_atol", "expected_rtol"])
+
+# Load previously generated JSON results, or init a new dict if none exist
+if os.path.isfile(results_filepath):
+    res_dict = json.load(open(results_filepath, mode="r"))
+else:
+    res_dict = dict()
+
+# Frame out data structure. Store data with the function type as a top level
+# key to allow future expansion
+res_dict["COMMENT"] = ("!!!!!! THIS FILE WAS AUTOGENERATED BY RUNNING "
+                       "`python gen_studentized_range_mpmath_ref.py` !!!!!!")
+res_dict.setdefault("cdf_data", [])
+res_dict.setdefault("pdf_data", [])
+res_dict.setdefault("moment_data", [])
+
+general_atol, general_rtol = 1e-11, 1e-11
+large_nu_atol = 1e-9  # Increase tolorances for large values of nu
+
+mp.dps = 24
+
+cp_q = [0.1, 1, 10]
+cp_k = [3, 10]
+cp_nu = [10, 100, 1000, np.inf]
+
+cdf_pdf_cases = [
+    CdfCase(*case,
+            general_atol if case[2] < 120 else large_nu_atol,
+            general_rtol)
+    for case in
+    itertools.product(cp_q, cp_k, cp_nu)
+]
+
+mom_atol, mom_rtol = 1e-9, 1e-9
+# These are EXTREMELY slow - Multiple days each in worst case.
+moment_cases = [
+    MomentCase(i, 3, 10, mom_atol, mom_rtol)
+    for i in range(5)
+]
+
+
+def write_data():
+    """Writes the current res_dict to the target JSON file"""
+    with open(results_filepath, mode="w") as f:
+        json.dump(res_dict, f, indent=2)
+
+
+def to_dict(named_tuple):
+    """Converts a namedtuple to a dict"""
+    return dict(named_tuple._asdict())
+
+
+def mp_res_to_dict(mp_result):
+    """Formats an MPResult namedtuple into a dict for JSON dumping"""
+    return {
+        "src_case": to_dict(mp_result.src_case),
+
+        # np assert can't handle mpf, so take the accuracy hit here.
+        "mp_result": float(mp_result.mp_result)
+    }
+
+
+def cdf_mp(q, k, nu):
+    """Straightforward implementation of studentized range CDF"""
+    q, k, nu = mpf(q), mpf(k), mpf(nu)
+
+    def inner(s, z):
+        return phi(z) * (Phi(z + q * s) - Phi(z)) ** (k - 1)
+
+    def outer(s, z):
+        return s ** (nu - 1) * phi(sqrt(nu) * s) * inner(s, z)
+
+    def whole(s, z):
+        return (sqrt(2 * pi) * k * nu ** (nu / 2)
+                / (gamma(nu / 2) * 2 ** (nu / 2 - 1)) * outer(s, z))
+
+    res = quad(whole, [0, inf], [-inf, inf],
+               method="gauss-legendre", maxdegree=10)
+    return res
+
+
+def pdf_mp(q, k, nu):
+    """Straightforward implementation of studentized range PDF"""
+    q, k, nu = mpf(q), mpf(k), mpf(nu)
+
+    def inner(s, z):
+        return phi(z + q * s) * phi(z) * (Phi(z + q * s) - Phi(z)) ** (k - 2)
+
+    def outer(s, z):
+        return s ** nu * phi(sqrt(nu) * s) * inner(s, z)
+
+    def whole(s, z):
+        return (sqrt(2 * pi) * k * (k - 1) * nu ** (nu / 2)
+                / (gamma(nu / 2) * 2 ** (nu / 2 - 1)) * outer(s, z))
+
+    res = quad(whole, [0, inf], [-inf, inf],
+               method="gauss-legendre", maxdegree=10)
+    return res
+
+
+def moment_mp(m, k, nu):
+    """Implementation of the studentized range moment"""
+    m, k, nu = mpf(m), mpf(k), mpf(nu)
+
+    def inner(q, s, z):
+        return phi(z + q * s) * phi(z) * (Phi(z + q * s) - Phi(z)) ** (k - 2)
+
+    def outer(q, s, z):
+        return s ** nu * phi(sqrt(nu) * s) * inner(q, s, z)
+
+    def pdf(q, s, z):
+        return (sqrt(2 * pi) * k * (k - 1) * nu ** (nu / 2)
+                / (gamma(nu / 2) * 2 ** (nu / 2 - 1)) * outer(q, s, z))
+
+    def whole(q, s, z):
+        return q ** m * pdf(q, s, z)
+
+    res = quad(whole, [0, inf], [0, inf], [-inf, inf],
+               method="gauss-legendre", maxdegree=10)
+    return res
+
+
+def result_exists(set_key, case):
+    """Searches the results dict for a result in the set that matches a case.
+    Returns True if such a case exists."""
+    if set_key not in res_dict:
+        raise ValueError(f"{set_key} not present in data structure!")
+
+    case_dict = to_dict(case)
+    existing_res = list(filter(
+        lambda res: res["src_case"] == case_dict,  # dict comparison
+        res_dict[set_key]))
+
+    return len(existing_res) > 0
+
+
+def run(case, run_lambda, set_key, index=0, total_cases=0):
+    """Runs the single passed case, returning an mp dictionary and index"""
+    t_start = time.perf_counter()
+
+    res = run_lambda(case)
+
+    print(f"Finished {index + 1}/{total_cases} in batch. "
+          f"(Took {time.perf_counter() - t_start}s)")
+
+    return index, set_key, mp_res_to_dict(MPResult(case, res))
+
+
+def write_result(res):
+    """A callback for completed jobs. Inserts and writes a calculated result
+     to file."""
+    index, set_key, result_dict = res
+    res_dict[set_key].insert(index, result_dict)
+    write_data()
+
+
+def run_cases(cases, run_lambda, set_key):
+    """Runs an array of cases and writes to file"""
+    # Generate jobs to run from cases that do not have a result in
+    # the previously loaded JSON.
+    job_arg = [(case, run_lambda, set_key, index, len(cases))
+               for index, case in enumerate(cases)
+               if not result_exists(set_key, case)]
+
+    print(f"{len(cases) - len(job_arg)}/{len(cases)} cases won't be "
+          f"calculated because their results already exist.")
+
+    jobs = []
+    pool = Pool(num_pools)
+
+    # Run all using multiprocess
+    for case in job_arg:
+        jobs.append(pool.apply_async(run, args=case, callback=write_result))
+
+    pool.close()
+    pool.join()
+
+
+def run_pdf(case):
+    return pdf_mp(case.q, case.k, case.v)
+
+
+def run_cdf(case):
+    return cdf_mp(case.q, case.k, case.v)
+
+
+def run_moment(case):
+    return moment_mp(case.m, case.k, case.v)
+
+
+def main():
+    t_start = time.perf_counter()
+
+    total_cases = 2 * len(cdf_pdf_cases) + len(moment_cases)
+    print(f"Processing {total_cases} test cases")
+
+    print(f"Running 1st batch ({len(cdf_pdf_cases)} PDF cases). "
+          f"These take about 30s each.")
+    run_cases(cdf_pdf_cases, run_pdf, "pdf_data")
+
+    print(f"Running 2nd batch ({len(cdf_pdf_cases)} CDF cases). "
+          f"These take about 30s each.")
+    run_cases(cdf_pdf_cases, run_cdf, "cdf_data")
+
+    print(f"Running 3rd batch ({len(moment_cases)} moment cases). "
+          f"These take about anywhere from a few hours to days each.")
+    run_cases(moment_cases, run_moment, "moment_data")
+
+    print(f"Test data generated in {time.perf_counter() - t_start}s")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/scipy/stats/tests/test_continuous_basic.py b/scipy/stats/tests/test_continuous_basic.py
index 3666a4a8472e..aaa5692aeff5 100644
--- a/scipy/stats/tests/test_continuous_basic.py
+++ b/scipy/stats/tests/test_continuous_basic.py
@@ -38,20 +38,24 @@
             'geninvgauss', 'genhyperbolic']
 # distslow are sorted by speed (very slow to slow)
 
+distxslow = ['studentized_range']
+# distxslow are sorted by speed (very slow to slow)
+
 # skip check_fit_args (test is slow)
 skip_fit_test_mle = ['exponpow', 'exponweib', 'gausshyper', 'genexpon',
                      'halfgennorm', 'gompertz', 'johnsonsb', 'johnsonsu',
                      'kappa4', 'ksone', 'kstwo', 'kstwobign', 'mielke', 'ncf',
                      'nct', 'powerlognorm', 'powernorm', 'recipinvgauss',
                      'trapezoid', 'vonmises', 'vonmises_line', 'levy_stable',
-                     'rv_histogram_instance']
+                     'rv_histogram_instance', 'studentized_range']
 
 # these were really slow in `test_fit`.py.
 # note that this list is used to skip both fit_test and fit_fix tests
 slow_fit_test_mm = ['argus', 'exponpow', 'exponweib', 'gausshyper', 'genexpon',
                     'genhalflogistic', 'halfgennorm', 'gompertz', 'johnsonsb',
                     'kappa4', 'kstwobign', 'recipinvgauss', 'skewnorm',
-                    'trapezoid', 'truncexpon', 'vonmises', 'vonmises_line']
+                    'trapezoid', 'truncexpon', 'vonmises', 'vonmises_line',
+                    'studentized_range']
 # pearson3 fails due to something weird
 # the first list fails due to non-finite distribution moments encountered
 # most of the rest fail due to integration warnings
@@ -74,7 +78,8 @@
                          'johnsonsu', 'kappa4', 'ksone', 'kstwo', 'kstwobign',
                          'levy_stable', 'mielke', 'ncf', 'ncx2',
                          'powerlognorm', 'powernorm', 'rdist', 'recipinvgauss',
-                         'trapezoid', 'vonmises', 'vonmises_line']
+                         'trapezoid', 'vonmises', 'vonmises_line',
+                         'studentized_range']
 # the first list fails due to non-finite distribution moments encountered
 # most of the rest fail due to integration warnings
 # pearson3 is overriden as not implemented due to gh-11746
@@ -103,7 +108,8 @@
                    'logistic', 'loguniform', 'maxwell', 'nakagami',
                    'ncf', 'nct', 'ncx2', 'norminvgauss', 'pearson3', 'rdist',
                    'reciprocal', 'rice', 'skewnorm', 't', 'tukeylambda',
-                   'vonmises', 'vonmises_line', 'rv_histogram_instance'])
+                   'vonmises', 'vonmises_line', 'rv_histogram_instance',
+                   'studentized_range'])
 
 _h = np.histogram([1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6,
                    6, 6, 6, 7, 7, 7, 8, 8, 9], bins=8)
@@ -116,6 +122,8 @@ def cases_test_cont_basic():
             continue
         elif distname in distslow:
             yield pytest.param(distname, arg, marks=pytest.mark.slow)
+        elif distname in distxslow:
+            yield pytest.param(distname, arg, marks=pytest.mark.xslow)
         else:
             yield distname, arg
 
@@ -242,6 +250,13 @@ def cases_test_moments():
         if distname == 'levy_stable':
             continue
 
+        if distname == 'studentized_range':
+            msg = ("studentized_range is far too slow for this test and it is "
+                   "redundant with test_distributions::TestStudentizedRange::"
+                   "test_moment_against_mp")
+            yield pytest.param(distname, arg, True, True, True,
+                               marks=pytest.mark.xslow(reason=msg))
+            continue
         cond1 = distname not in fail_normalization
         cond2 = distname not in fail_higher
 
@@ -287,7 +302,7 @@ def test_moments(distname, arg, normalization_ok, higher_ok, is_xfailing):
 
 @pytest.mark.parametrize('dist,shape_args', distcont)
 def test_rvs_broadcast(dist, shape_args):
-    if dist in ['gausshyper', 'genexpon']:
+    if dist in ['gausshyper', 'genexpon', 'studentized_range']:
         pytest.skip("too slow")
 
     # If shape_only is True, it means the _rvs method of the
@@ -543,6 +558,7 @@ def check_pdf_logpdf_at_endpoints(distfn, args, msg):
             "invalid value encountered in add",  # genextreme
             "invalid value encountered in subtract",  # gengamma
             "invalid value encountered in multiply"  # recipinvgauss
+            # "invalid value encountered in log"  # studentized_range
             ]
         for msg in suppress_messsages:
             sup.filter(category=RuntimeWarning, message=msg)
diff --git a/scipy/stats/tests/test_distributions.py b/scipy/stats/tests/test_distributions.py
index 9850def7efe7..2750c8506392 100644
--- a/scipy/stats/tests/test_distributions.py
+++ b/scipy/stats/tests/test_distributions.py
@@ -6,6 +6,7 @@
 import sys
 import pickle
 import os
+import json
 
 from numpy.testing import (assert_equal, assert_array_equal,
                            assert_almost_equal, assert_array_almost_equal,
@@ -20,7 +21,8 @@
 from numpy.lib.recfunctions import rec_append_fields
 from scipy import special
 from scipy._lib._util import check_random_state
-from scipy.integrate import IntegrationWarning, quad, trapezoid
+from scipy.integrate import (IntegrationWarning, quad, trapezoid,
+                             cumulative_trapezoid)
 import scipy.stats as stats
 from scipy.stats._distn_infrastructure import argsreduce
 import scipy.stats.distributions
@@ -30,6 +32,7 @@
 from .test_discrete_basic import distdiscrete, invdistdiscrete
 from scipy.stats._continuous_distns import FitDataError
 from scipy.optimize import root, fmin
+from itertools import product
 
 # python -OO strips docstrings
 DOCSTRINGS_STRIPPED = sys.flags.optimize > 1
@@ -4490,6 +4493,131 @@ def test_burr_nan_mean_var_9544(self):
         assert_(np.isfinite(e4))
 
 
+class TestStudentizedRange:
+    # For alpha = .05, .01, and .001, and for each value of
+    # v = [1, 3, 10, 20, 120, inf], a Q was picked from each table for
+    # k = [2, 8, 14, 20].
+
+    # these arrays are written with `k` as column, and `v` as rows.
+    # Q values are taken from table 3:
+    # https://www.jstor.org/stable/2237810
+    q05 = [17.97, 45.40, 54.33, 59.56,
+           4.501, 8.853, 10.35, 11.24,
+           3.151, 5.305, 6.028, 6.467,
+           2.950, 4.768, 5.357, 5.714,
+           2.800, 4.363, 4.842, 5.126,
+           2.772, 4.286, 4.743, 5.012]
+    q01 = [90.03, 227.2, 271.8, 298.0,
+           8.261, 15.64, 18.22, 19.77,
+           4.482, 6.875, 7.712, 8.226,
+           4.024, 5.839, 6.450, 6.823,
+           3.702, 5.118, 5.562, 5.827,
+           3.643, 4.987, 5.400, 5.645]
+    q001 = [900.3, 2272, 2718, 2980,
+            18.28, 34.12, 39.69, 43.05,
+            6.487, 9.352, 10.39, 11.03,
+            5.444, 7.313, 7.966, 8.370,
+            4.772, 6.039, 6.448, 6.695,
+            4.654, 5.823, 6.191, 6.411]
+    qs = np.concatenate((q05, q01, q001))
+    ps = [.95, .99, .999]
+    vs = [1, 3, 10, 20, 120, np.inf]
+    ks = [2, 8, 14, 20]
+
+    data = zip(product(ps, vs, ks), qs)
+
+    def test_cdf_against_tables(self):
+        for pvk, q in self.data:
+            p_expected, v, k = pvk
+            res_p = stats.studentized_range.cdf(q, k, v)
+            assert_allclose(res_p, p_expected, rtol=1e-4)
+
+    @pytest.mark.slow
+    def test_ppf_against_tables(self):
+        for pvk, q_expected in self.data:
+            res_q = stats.studentized_range.ppf(*pvk)
+            assert_allclose(res_q, q_expected, rtol=1e-4)
+
+    path_prefix = os.path.dirname(__file__)
+    relative_path = "data/studentized_range_mpmath_ref.json"
+    with open(os.path.join(path_prefix, relative_path), "r") as file:
+        pregenerated_data = json.load(file)
+
+    @pytest.mark.parametrize("case_result", pregenerated_data["cdf_data"])
+    def test_cdf_against_mp(self, case_result):
+        src_case = case_result["src_case"]
+        mp_result = case_result["mp_result"]
+        qkv = src_case["q"], src_case["k"], src_case["v"]
+        res = stats.studentized_range.cdf(*qkv)
+
+        assert_allclose(res, mp_result,
+                        atol=src_case["expected_atol"],
+                        rtol=src_case["expected_rtol"])
+
+    @pytest.mark.parametrize("case_result", pregenerated_data["pdf_data"])
+    def test_pdf_against_mp(self, case_result):
+        src_case = case_result["src_case"]
+        mp_result = case_result["mp_result"]
+        qkv = src_case["q"], src_case["k"], src_case["v"]
+        res = stats.studentized_range.pdf(*qkv)
+
+        assert_allclose(res, mp_result,
+                        atol=src_case["expected_atol"],
+                        rtol=src_case["expected_rtol"])
+
+    @pytest.mark.slow
+    @pytest.mark.parametrize("case_result", pregenerated_data["moment_data"])
+    def test_moment_against_mp(self, case_result):
+        src_case = case_result["src_case"]
+        mp_result = case_result["mp_result"]
+        mkv = src_case["m"], src_case["k"], src_case["v"]
+        res = stats.studentized_range.moment(*mkv)
+
+        assert_allclose(res, mp_result,
+                        atol=src_case["expected_atol"],
+                        rtol=src_case["expected_rtol"])
+
+    def test_pdf_integration(self):
+        k, v = 3, 10
+        # Test whether PDF integration is 1 like it should be.
+        res = quad(stats.studentized_range.pdf, 0, np.inf, args=(k, v))
+        assert_allclose(res[0], 1)
+
+    @pytest.mark.xslow
+    def test_pdf_against_cdf(self):
+        k, v = 3, 10
+
+        # Test whether the integrated PDF matches the CDF using cumulative
+        # integration. Use a small step size to reduce error due to the
+        # summation. This is slow, but tests the results well.
+        x = np.arange(0, 10, step=0.01)
+
+        y_cdf = stats.distributions.studentized_range.cdf(x, k, v)[1:]
+        y_pdf_raw = stats.distributions.studentized_range.pdf(x, k, v)
+        y_pdf_cumulative = cumulative_trapezoid(y_pdf_raw, x)
+
+        # Because of error caused by the summation, use a relatively large rtol
+        assert_allclose(y_pdf_cumulative, y_cdf, rtol=1e-4)
+
+    @pytest.mark.slow
+    def test_moment_vectorization(self):
+        # Test moment broadcasting. Calls `_munp` directly because
+        # `rv_continuous.moment` is broken at time of writing. See gh-12192
+        m = stats.studentized_range._munp([1, 2], [4, 5], [10, 11])
+        assert_allclose(m.shape, (2,))
+
+        with pytest.raises(ValueError, match="...could not be broadcast..."):
+            stats.studentized_range._munp(1, [4, 5], [10, 11, 12])
+
+    @pytest.mark.xslow
+    def test_fitstart_valid(self):
+        with suppress_warnings() as sup, np.errstate(invalid="ignore"):
+            # the integration warning message may differ
+            sup.filter(IntegrationWarning)
+            k, df, _, _ = stats.studentized_range._fitstart([1, 2, 3])
+        assert_(stats.studentized_range._argcheck(k, df))
+
+
 def test_540_567():
     # test for nan returned in tickets 540, 567
     assert_almost_equal(stats.norm.cdf(-1.7624320982), 0.03899815971089126,
diff --git a/scipy/stats/tests/test_fit.py b/scipy/stats/tests/test_fit.py
index 832763c1d15b..964cc8238550 100644
--- a/scipy/stats/tests/test_fit.py
+++ b/scipy/stats/tests/test_fit.py
@@ -34,6 +34,7 @@
         'vonmises',
         'levy_stable',
         'trapezoid',
+        'studentized_range'
 ]
 
 mm_failing_fits = ['alpha', 'betaprime', 'burr', 'burr12', 'cauchy', 'chi',
@@ -45,7 +46,7 @@
                    'levy_stable', 'loglaplace', 'lomax', 'mielke', 'nakagami',
                    'ncf', 'nct', 'ncx2', 'pareto', 'powerlognorm', 'powernorm',
                    'skewcauchy', 't',
-                   'trapezoid', 'triang', 'tukeylambda']
+                   'trapezoid', 'triang', 'tukeylambda', 'studentized_range']
 
 # not sure if these fail, but they caused my patience to fail
 mm_slow_fits = ['argus', 'exponpow', 'exponweib', 'gausshyper', 'genexpon',
diff --git a/scipy/stats/tests/test_stats.py b/scipy/stats/tests/test_stats.py
index beb26c4de0c1..3a823811b904 100644
--- a/scipy/stats/tests/test_stats.py
+++ b/scipy/stats/tests/test_stats.py
@@ -30,6 +30,7 @@
 from .common_tests import check_named_results
 from scipy.sparse.sputils import matrix
 from scipy.spatial.distance import cdist
+from itertools import permutations
 from numpy.lib import NumpyVersion
 from scipy.stats.stats import (_broadcast_concatenate,
                                AlexanderGovernConstantInputWarning)
@@ -6089,6 +6090,212 @@ def test_bad_shapes(self):
             stats.f_oneway(a, b, axis=1)
 
 
+class TestTukeykramer:
+
+    data_same_size = ([24.5, 23.5, 26.4, 27.1, 29.9],
+                      [28.4, 34.2, 29.5, 32.2, 30.1],
+                      [26.1, 28.3, 24.3, 26.2, 27.8])
+    data_diff_size = ([24.5, 23.5, 26.28, 26.4, 27.1, 29.9, 30.1, 30.1],
+                      [28.4, 34.2, 29.5, 32.2, 30.1],
+                      [26.1, 28.3, 24.3, 26.2, 27.8])
+    extreme_size = ([24.5, 23.5, 26.4],
+                    [28.4, 34.2, 29.5, 32.2, 30.1, 28.4, 34.2, 29.5, 32.2,
+                     30.1],
+                    [26.1, 28.3, 24.3, 26.2, 27.8])
+    sas_same_size = [
+        ['2', '3', '4.340', '.6908830567505320', '7.989116943249460', '1'],
+        ['2', '1', '4.600', '.9508830567505330', '8.249116943249470', '1'],
+        ['3', '2', '-4.340', '-7.98911694324947', '-.690883056750532', '1'],
+        ['3', '1', '0.260', '-3.38911694324947', '3.909116943249470', '0'],
+        ['1', '2', '-4.600', '-8.24911694324947', '-.950883056750534', '1'],
+        ['1', '3', '-0.260', '-3.90911694324947', '3.389116943249460', '0']]
+
+    sas_diff_size = [
+        ['2', '1', '3.645', '0.2679292644988780', '7.022070735501120', '1'],
+        ['2', '3', '4.340', '0.5934764007020940', '8.086523599297900', '1'],
+        ['1', '2', '-3.645', '-7.02207073550113', '-.267929264498878', '1'],
+        ['1', '3', '.6949999999999960', '-2.68207073550113',
+         '4.072070735501120', '0'],
+        ['3', '2', '-4.340', '-8.08652359929791', '-.593476400702094', '1'],
+        ['3', '1', '-.694999999999997', ' -4.07207073550112',
+        '2.682070735501120', '0']]
+
+    sas_extreme = [
+        ['2', '3', '4.34', '1.561605075444930',  '7.118394924555070', '1'],
+        ['2', '1', '6.08', '2.740784878675930', '9.419215121324070', '1'],
+        ['3', '2', '-4.34', '-7.11839492455507',  '-1.56160507544493', '1'],
+        ['3', '1', '1.74', '-1.96452656607342',  '5.444526566073420', '0'],
+        ['1', '2', '-6.08', '-9.41921512132407', '-2.74078487867593', '1'],
+        ['1', '3', '-1.74', '-5.44452656607343',  '1.964526566073420', '0']
+    ]
+
+    @pytest.mark.parametrize("data,SAS_res", ((data_same_size, sas_same_size),
+                                              (data_diff_size, sas_diff_size),
+                                              (extreme_size, sas_extreme)),
+                             ids=["equal size sample", "unequal sample size",
+                                  "extreme sample size differences"])
+    def test_compare_sas(self, data, SAS_res):
+        '''
+        SAS code used to generate results for each sample:
+
+        DATA ACHE;
+        INPUT BRAND RELIEF;
+        CARDS;
+        1 24.5
+        1 23.5
+        1 26.28
+        1 26.4
+        1 27.1
+        1 29.9
+        1 30.1
+        1 30.1
+        2 28.4
+        2 34.2
+        2 29.5
+        2 32.2
+        2 30.1
+        3 26.1
+        3 28.3
+        3 24.3
+        3 26.2
+        3 27.8
+        ;
+        ods graphics on;   ODS RTF;ODS LISTING CLOSE;
+           PROC ANOVA DATA=ACHE;
+           CLASS BRAND;
+           MODEL RELIEF=BRAND;
+           MEANS BRAND/TUKEY CLDIFF;
+           TITLE 'COMPARE RELIEF ACROSS MEDICINES  - ANOVA EXAMPLE';
+           ods output  CLDiffs =tc;
+        proc print data=tc;
+            format LowerCL 17.16 UpperCL 17.16 Difference 17.16;
+            title “Output with many digits”;
+        RUN;
+        QUIT;
+        ODS RTF close;
+        ODS LISTING;
+
+        # The standard output, not showing the additional digits. The code
+        # example above outputs a table with more digits.
+
+        Brand Comparison    Mean Difference  Simultanious 95% Confidence Limits
+        2 - 1               3.645               0.268           7.022   1
+        2 - 3               4.340               0.593           8.087   1
+        1 - 2              -3.645              -7.022          -0.268   1
+        1 - 3               0.695              -2.682           4.072   0
+        3 - 2              -4.340              -8.087          -0.593   1
+        3 - 1              -0.695              -4.072           2.682   0
+        '''
+
+        res = stats.tukeykramer(*data)
+        conf = res.confidence_interval()
+        for comp in SAS_res:
+            i, j, stat, lower, upper, sig = np.asarray(comp, dtype=float)
+            # make int for slicing, zero index
+            i, j = int(i) - 1, int(j) - 1
+            assert_allclose(stat, res.statistic[i, j])
+            assert_allclose(lower, conf.low[i, j], atol=1e-4)
+            assert_allclose(upper, conf.high[i, j], atol=1e-4)
+            assert ((sig == 1 and res.pvalue[i, j] < .05) or
+                    (sig == 0 and res.pvalue[i, j] > .05))
+
+
+    def test_comp_R(self):
+        '''  
+        Testing against results and p-values from R:
+        from: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/TukeyHSD
+
+        > require(graphics)
+
+        > summary(fm1 <- aov(breaks ~ tension, data = warpbreaks))
+        > TukeyHSD(fm1, "tension", ordered = TRUE)
+        > plot(TukeyHSD(fm1, "tension"))
+
+        Tukey multiple comparisons of means
+        95% family-wise confidence level
+        factor levels have been ordered
+
+        Fit: aov(formula = breaks ~ tension, data = warpbreaks)
+
+        $tension
+                 diff        lwr      upr     p adj
+        M-H  4.722222 -4.8376022 14.28205 0.4630831
+        L-H 14.722222  5.1623978 24.28205 0.0014315
+        L-M 10.000000  0.4401756 19.55982 0.0384598
+
+        '''
+        l = [26, 30, 54, 25, 70, 52, 51, 26, 67,
+             27, 14, 29, 19, 29, 31, 41, 20, 44]
+        m = [18, 21, 29, 17, 12, 18, 35, 30, 36,
+             42, 26, 19, 16, 39, 28, 21, 39, 29]
+        h = [36, 21, 24, 18, 10, 43, 28, 15, 26,
+             20, 21, 24, 17, 13, 15, 15, 16, 28]
+        res = stats.tukeykramer(l, m, h)
+
+        r_pvalue = np.asarray(
+            [[0, 0.0384598, 0.0014315],
+             [0, 0, 0.4630831],
+             [0, 0, 0]])
+        r_stat = np.asarray(
+            [[0, 10.000000, 14.722222],
+             [0, 0, 4.722222],
+             [0, 0, 0]])
+        r_lower = np.asarray(
+            [[0, 0.4401756, 5.1623978],
+             [0, 0, -4.8376022],
+             [0, 0, 0]])
+        r_upper = np.asarray(
+            [[0, 19.55982, 24.28205],
+             [0, 0, 14.28205],
+             [0, 0, 0]])
+
+        conf = res.confidence_interval()
+        atol = 1e-3
+        for (i, j) in [(0, 1), (0, 2), (1, 2)]:
+            assert_allclose(r_pvalue[i, j], res.pvalue[i, j], atol=atol)
+            assert_allclose(r_stat[i, j], res.statistic[i, j], atol=atol)
+            assert_allclose(r_lower[i, j], conf.low[i, j], atol=atol)
+            assert_allclose(r_upper[i, j], conf.high[i, j], atol=atol)
+
+
+    def test_engineering_stat_handbook(self):
+        '''
+        Example sourced from:
+        https://www.itl.nist.gov/div898/handbook/prc/section4/prc471.htm
+        '''
+        group1 = [6.9, 5.4, 5.8, 4.6, 4.0]
+        group2 = [8.3, 6.8,	7.8, 9.2, 6.5]
+        group3 = [8.0, 10.5, 8.1, 6.9, 9.3]
+        group4 = [5.8, 3.8,	6.1, 5.6, 6.2]
+        res = stats.tukeykramer(group1, group2, group3, group4)
+        conf = res.confidence_interval()
+        lower = np.asarray([
+            [0, 0, 0, -2.25],
+            [.29, 0, -2.93, .13],
+            [1.13, 0, 0, .97],
+            [0, 0, 0, 0]])
+        upper = np.asarray([
+            [0, 0, 0, 1.93],
+            [4.47, 0,  1.25, 4.31],
+            [5.31, 0, 0, 5.15],
+            [0, 0, 0, 0]])
+
+        for (i, j) in [(1, 0), (2, 0), (0, 3), (1, 2), (2, 3)]:
+            assert_allclose(lower[i, j], conf.low[i, j], atol=1e-2)
+            assert_allclose(upper[i, j], conf.high[i, j], atol=1e-2)
+
+    def test_no_inf(self):
+        with assert_raises(ValueError, match="...must be finite."):
+            stats.tukeykramer([1, 2, 3], [2, np.inf])
+
+    def test_is_1d(self):
+        with assert_raises(ValueError, match="...must be one-dimensional"):
+            stats.tukeykramer([[1, 2,], [2, 3]], [2, 5])
+
+    def test_no_empty(self):
+        with assert_raises(ValueError, match="...must be greater than one"):
+            stats.tukeykramer([], [2, 5])
+
 class TestKruskal:
     def test_simple(self):
         x = [1]
diff --git a/tools/unicode-check.py b/tools/unicode-check.py
index 0a0182b634e7..6f4f1d981aa6 100755
--- a/tools/unicode-check.py
+++ b/tools/unicode-check.py
@@ -10,7 +10,7 @@
 # The set of Unicode code points greater than 127 that we
 # allow in the source code.
 box_drawing_chars = set(chr(cp) for cp in range(0x2500, 0x2580))
-allowed = (set(['®', 'é', 'ó', 'ö', 'ü', 'λ', 'π', 'ω', '∫', '≠']) |
+allowed = (set(['®', 'é', 'ê', 'ó', 'ö', 'ü', 'λ', 'π', 'ω', '∫', '≠', '≥']) |
            box_drawing_chars)