Repository: cjekel/tolerance_interval_py

Issue:
code exact solution for two sided normal distribution

Discussion:

Patch:
diff --git a/.gitignore b/.gitignore
index a2cc631..b4aec76 100644
--- a/.gitignore
+++ b/.gitignore
@@ -103,6 +103,10 @@ venv.bak/
 .spyderproject
 .spyproject
 
+# Sublime project settings
+*.sublime-project
+*.sublime-workspace
+
 # Rope project settings
 .ropeproject
 
diff --git a/environment.yml b/environment.yml
new file mode 100644
index 0000000..9236bce
--- /dev/null
+++ b/environment.yml
@@ -0,0 +1,15 @@
+# Simple environment for developing tolerance_interval_py
+# To use:
+#   $ conda env create -f environment.yml  # `mamba` works too for this command
+#   $ conda activate tolerance-interval-py
+#
+name: tolerance-interval-py
+channels:
+  - conda-forge
+dependencies:
+  - numpy
+  - scipy
+  - sympy
+  - setuptools
+  # Avoid pulling in large MKL libraries.
+  - nomkl
diff --git a/toleranceinterval/twoside/__init__.py b/toleranceinterval/twoside/__init__.py
index 3a60577..05d4a95 100644
--- a/toleranceinterval/twoside/__init__.py
+++ b/toleranceinterval/twoside/__init__.py
@@ -1,2 +1,3 @@
+from .twoside import normal_factor  # noqa F401
 from .twoside import normal  # noqa F401
 from .twoside import lognormal  # noqa F401
diff --git a/toleranceinterval/twoside/_normal_approx.py b/toleranceinterval/twoside/_normal_approx.py
new file mode 100644
index 0000000..88063e0
--- /dev/null
+++ b/toleranceinterval/twoside/_normal_approx.py
@@ -0,0 +1,130 @@
+# -- coding: utf-8 --
+# MIT License
+#
+# Copyright (c) 2019 Charles Jekel
+#
+# Permission is hereby granted, free of charge, to any person obtaining a copy
+# of this software and associated documentation files (the "Software"), to deal
+# in the Software without restriction, including without limitation the rights
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+# copies of the Software, and to permit persons to whom the Software is
+# furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in
+# all copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+r"""
+Algorithms for computing approximate two-sided statistical tolerance interval
+factors under the assumption of a normal distribution.
+
+"""
+
+import numpy as np
+from scipy.stats import norm, chi2
+
+
+def tolerance_factor_howe(n, p, g, m=None, nu=None):
+    r"""
+    Compute two-side central tolerance interval factor using Howe's method.
+
+    Computes the two-sided tolerance interval (TI) factor under a normal
+    distribution assumption using Howe's method. This follows the derivation
+    in [1]. This is an approximation, and does not represent the exact TI.
+
+    Parameters
+    ----------
+    n : scalar
+        Sample size.
+    p : float
+        Percentile for central TI to estimate.
+    g : float
+        Confidence level where g > 0. and g < 1.
+    m : scalar
+        Number of independent random samples (of size n). If None,
+        default value is m = 1.
+    nu : scalar
+        Degrees of freedom for distribution of the (pooled) sample
+        variance. If None, default value is nu = m*(n-1).
+
+    Returns
+    -------
+    float
+        The calculated tolerance factor for the tolerance interval.
+
+    References
+    ----------
+    .. [1] Howe, W. G. (1969). "Two-sided Tolerance Limits for Normal
+        Populations - Some Improvements", Journal of the American Statistical
+        Association, 64 , pages 610-620.
+
+    """
+    # Handle defaults for keyword inputs.
+    if m is None:
+        m = 1
+    if nu is None:
+        nu = m * (n - 1)
+
+    alpha = 1.0 - g
+    zp = norm.ppf((1.0 + p) / 2.0)
+    u = zp * np.sqrt(1.0 + (1.0 / n))
+    chi2_nu = chi2.ppf(alpha, df=nu)
+    v = np.sqrt(nu / chi2_nu)
+    k = u * v
+    return k
+
+
+def tolerance_factor_guenther(n, p, g, m=None, nu=None):
+    r"""
+    Compute two-side central tolerance interval factor using Guenther's method.
+
+    Computes the two-sided tolerance interval (TI) factor under a normal
+    distribution assumption using Guenthers's method. This follows the
+    derivation in [1]. This is an approximation, and does not represent the
+    exact TI.
+
+    Parameters
+    ----------
+    n : scalar
+        Sample size.
+    p : float
+        Percentile for central TI to estimate.
+    g : float
+        Confidence level where g > 0. and g < 1.
+    m : scalar
+        Number of independent random samples (of size n). If None,
+        default value is m = 1.
+    nu : scalar
+        Degrees of freedom for distribution of the (pooled) sample
+        variance. If None, default value is nu = m*(n-1).
+
+    Returns
+    -------
+    float
+        The calculated tolerance factor for the tolerance interval.
+
+    References
+    ----------
+    .. [1] Guenther, W. C. (1977). "Sampling Inspection in Statistical Quality
+        Control", Griffin's Statistical Monographs, Number 37, London.
+
+    """
+    # Handle defaults for keyword inputs.
+    if m is None:
+        m = 1
+    if nu is None:
+        nu = m * (n - 1)
+
+    k = tolerance_factor_howe(n, p, g, m, nu)
+    alpha = 1.0 - g
+    chi2_nu = chi2.ppf(alpha, df=nu)
+    w = np.sqrt(1.0 + ((n - 3.0 - chi2_nu) / (2.0 * (n + 1.0) ** 2)))
+    k *= w
+    return k
diff --git a/toleranceinterval/twoside/_normal_exact.py b/toleranceinterval/twoside/_normal_exact.py
new file mode 100644
index 0000000..6f67e84
--- /dev/null
+++ b/toleranceinterval/twoside/_normal_exact.py
@@ -0,0 +1,356 @@
+# Author: Copyright (c) 2021 Jed Ludlow
+# License: MIT License
+
+r"""
+Algorithm for computing the exact two-sided statistical tolerance interval
+factor under the assumption of a normal distribution.
+
+This module is a Python port of an algorithm written in MATLAB by
+Viktor Witkovsky and posted to the MATLAB Central File Exchange, retrieved
+2021-03-12 at this URL:
+
+https://www.mathworks.com/matlabcentral/fileexchange/24135-tolerancefactor
+
+Here is Witkovsky's original copyright and disclaimer:
+
+-------------------------------------------------------------------------------
+Copyright (c) 2020, Viktor Witkovsky
+Copyright (c) 2013, Viktor Witkovsky
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+* Redistributions of source code must retain the above copyright notice, this
+  list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright notice,
+  this list of conditions and the following disclaimer in the documentation
+  and/or other materials provided with the distribution
+* Neither the name of Institute of Measurement Science, Slovak Academy of
+  Sciences, Bratislava, Slovakia nor the names of its
+  contributors may be used to endorse or promote products derived from this
+  software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+-------------------------------------------------------------------------------
+"""
+
+import numpy as np
+import scipy.stats as stats
+import scipy.special as spec
+import scipy.integrate as integ
+import scipy.optimize as optim
+
+SQRT_2 = 1.4142135623730950488
+SQRT_2PI = 2.5066282746310005024
+
+
+def tolerance_factor(n, coverage, confidence, m=None, nu=None, d2=None,
+                     simultaneous=False, tailprob=False):
+    r"""
+    Compute the exact tolerance factor k for the two-sided statistical
+    tolerance interval to contain a proportion of a population with given
+    confidence under normal distribution assumptions.
+
+    Computes (by Gauss-Kronod quadrature) the exact tolerance factor k for the
+    two-sided coverage-content and (1-alpha)-confidence tolerance interval
+
+        TI = [Xmean - k * S, Xmean + k * S]
+
+    where Xmean = mean(X), S = std(X), X = [X_1,...,X_n] is a random sample
+    of size n from the distribution N(mu,sig2) with unknown mean mu and
+    variance sig2.
+
+    The tolerance factor k is determined such that the tolerance intervals
+    with confidence (1-alpha) cover at least the coverage fraction
+    of the distribution N(mu,sigma^2), i.e.
+    Prob[ Prob( Xmean - k * S < X < Xmean + k * S ) >= p ]= 1-alpha,
+    for X ~ N(mu,sig2) which is independent with Xmean and S. For more
+    details see e.g. Krishnamoorthy and Mathew (2009).
+
+    REMARK:
+    If S is a pooled estimator of sig, based on m random samples of size n,
+    tolerance_factor computes the simultaneous (or non-simultaneous) exact
+    tolerance factor k for the two-sided coverage-content and
+    (1-alpha)-confidence tolerance intervals
+
+        TI = [Xmean_i - k * S, Xmean_i + k * S], for i = 1,...,m
+
+    where Xmean_i = mean(X_i), X_i = [X_i1,...,X_in] is a random sample of
+    size n from the distribution N(mu,sig2) with unknown mean mu and
+    variance sig2, and S = sqrt(S2), S2 is the pooled estimator of sig2,
+
+        S2 = (1/nu) * sum_i=1:m ( sum_j=1:n (X_ij - Xmean_i)^2 )
+
+    with nu degrees of freedom, nu = m * (n-1).
+
+    Parameters
+    ----------
+    n : scalar
+        Sample size
+    coverage : scalar in the interval [0.0, 1.0]
+        Coverage (or content) probability,
+        Prob( Xmean - k * S < X < Xmean + k * S ) >= coverage
+    confidence : scalar in the interval [0.0, 1.0]
+        Confidence probability,
+        Prob[ Prob( Xmean-k*S < X < Xmean+k*S ) >= coverage ] = confidence.
+    m : scalar
+        Number of independent random samples (of size n). If None,
+        default value is m = 1.
+    nu : scalar
+        Degrees of freedom for distribution of the (pooled) sample
+        variance S2. If None, default value is nu = m*(n-1).
+    d2 : scalar
+        Normalizing constant. For computing the factors of the
+        non-simultaneous tolerance limits (xx'*betaHat +/- k * S)
+        for the linear regression y = XX*beta +epsilon, set d2 =
+        xx'*inv(XX'*XX)*xx.
+        Typically, in simple linear regression the estimator S2 has
+        nu = n-2 degrees of freedom. If None, default value is d2 = 1/n.
+    simultaneous : boolean
+        Logical flag for calculating the factor for
+        simultaneous tolerance intervals. If False, tolerance_factor will
+        calculate the factor for the non-simultaneous tolerance interval.
+        Default value is False.
+    tailprob : boolean
+        Logical flag for representing the input probabilities
+        'coverage' and 'confidence'. If True, the input parameters are
+        represented as the tail coverage (i.e. 1 - coverage) and tail confidence
+        (i.e. 1 - confidence). This option is useful if the interest is to
+        calculate the tolerance factor for extremely large values
+        of coverage and/or confidence, close to 1, as e.g. coverage = 1 - 1e-18.
+        Default value is False.
+
+    Returns
+    -------
+    float
+        The calculated tolerance factor for the tolerance interval.
+
+    References
+    ----------
+    [1] Krishnamoorthy K, Mathew T. (2009). Statistical Tolerance Regions:
+        Theory, Applications, and Computation. John Wiley & Sons, Inc.,
+        Hoboken, New Jersey. ISBN: 978-0-470-38026-0, 512 pages.
+    [2] Witkovsky V. On the exact two-sided tolerance intervals for
+        univariate normal distribution and linear regression. Austrian
+        Journal of Statistics 43(4), 2014, 279-92. http://
+        ajs.data-analysis.at/index.php/ajs/article/viewFile/vol43-4-6/35
+    [3] ISO 16269-6:2014: Statistical interpretation of data - Part 6:
+        Determination of statistical tolerance intervals.
+    [4] Janiga I., Garaj I.: Two-sided tolerance limits of normal
+        distributions with unknown means and unknown common variability.
+        MEASUREMENT SCIENCE REVIEW, Volume 3, Section 1, 2003, 75-78.
+    [5] Robert W. Mee (1990) Simultaneous Tolerance Intervals for Normal
+        Populations With Common Variance, Technometrics, 32:1, 83-92,
+        DOI: 10.1080/00401706.1990.10484595
+
+    """
+    # Coverage and confidence must be within [0.0, 1.0].
+    if (coverage < 0.0) or (coverage > 1.0):
+        raise ValueError('coverage must be within the interval [0.0, 1.0].')
+    if (confidence < 0.0) or (confidence > 1.0):
+        raise ValueError('confidence must be within the interval [0.0, 1.0].')
+
+    # Handle defaults for keyword inputs.
+    if m is None:
+        m = 1
+    if nu is None:
+        nu = m * (n - 1)
+    if d2 is None:
+        d2 = 1.0 / n
+
+    # Set values for limiting cases of coverage and confidence.
+    if confidence == 0:
+        k = np.nan
+        return k
+    if confidence == 1:
+        k = np.inf
+        return k
+    if coverage == 1:
+        k = np.inf
+        return k
+    if coverage == 0:
+        k = 0.0
+        return k
+
+    # Set the constants.
+    if tailprob:
+        tailconfidence = confidence
+        tailcoverage = coverage
+    else:
+        tailconfidence = round(1.0e+16 * (1.0 - confidence)) / 1.0e+16
+        tailcoverage = round(1.0e+16 * (1.0 - coverage)) / 1.0e+16
+
+    # Return the result for nu = Inf.
+    if nu == np.inf:
+        k = SQRT_2 * spec.erfcinv(tailcoverage)
+        return k
+
+    # Handle the case of nu not being positive.
+    if nu <= 0:
+        raise(ValueError, 'Degrees of freedom should be positive.')
+
+    # Compute the two-sided tolerance factor
+    tol_high_precision = np.spacing(tailconfidence)
+
+    # Integration limits
+    A = 0.0
+    B = 10.0
+
+    if m > 1:
+        k0, _ = _approx_tol_factor_witkovsky(
+            tailcoverage, tailconfidence, d2, m, nu, A, B)
+    else:
+        k0, _ = _approx_tol_factor_wald_wolfowitz(
+            tailcoverage, tailconfidence, d2, nu)
+
+    # Compute the tolerance factor.
+    sol = optim.newton(lambda k: (_integral_gauss_kronod(
+        k, nu, m, d2, tailcoverage, simultaneous, A, B, tol_high_precision)
+        - tailconfidence), k0)
+    k = sol
+
+    return k
+
+
+def _integral_gauss_kronod(k, nu, m, c, tailcoverage, simultaneous, A, B, tol):
+    r"""
+    Evaluates the integral defined by eqs. (1.2.4) and (2.5.8) in
+    Krishnamoorthy and Mathew: Statistical Tolerance Regions, Wiley, 2009,
+    (pp.7 and 52), by adaptive Gauss-Kronod quadrature.
+
+    (The tol argument is currently ignored but may be used at some point. It
+    was used in Witkovsky's integration, but we ignore it here and use
+    SciPy default tolerances, which seem to be adequate.)
+
+    """
+    val, _ = integ.quad(lambda z: _integrand(
+        z, k, nu, m, c, tailcoverage, simultaneous), A, B)
+    val *= 2
+    return val
+
+
+def _integrand(z, k, nu, m, c, tailcoverage, simultaneous):
+    r"""
+    Integrand for Gauss-Kronod quadrature.
+
+    """
+    root = _find_root(np.sqrt(c) * z, tailcoverage)
+    ncx2pts = nu * root**2
+    factor = np.exp(-0.5 * z**2) / SQRT_2PI
+
+    if simultaneous:
+        factor = factor * (m * (1 - (spec.erfc(z / SQRT_2))) ** (m - 1))
+
+    x = ncx2pts / k ** 2
+    fun = spec.gammainc(nu / 2.0, x / 2.0) * factor
+    return fun
+
+
+def _find_root(x, tailcoverage):
+    r"""
+    Numerically finds the solution (root), of the equation
+
+    normcdf(x+root) - normcdf(x-root) = coverage = 1 - tailcoverage
+
+    by Halley's method for finding the root of the function
+
+    fun(r) = fun(r|x,tailcoverage)
+
+    based on two derivatives, funD1(r|x,tailcoverage)
+    and funD2(r|x,tailcoverage) of the fun(r|x,tailcoverage), (for given x
+    and tailcoverage).
+
+    Note that r = sqrt(ncx2inv(1-tailcoverage,1,x^2)), where ncx2inv is the
+    inverse of the noncentral chi-square cdf.
+
+    """
+    # Set the constants
+    max_iterations = 100
+    iteration = 0
+
+    # Set the appropriate tolerance
+    if np.spacing(tailcoverage) < np.spacing(1):
+        tol = min(10 * np.spacing(tailcoverage), np.spacing(1))
+
+    # Set the starting value of the root r: r0 = x + norminv(coverage)
+    r = x + SQRT_2 * spec.erfcinv(2 * tailcoverage)
+
+    # Main loop (Halley's method)
+    while True:
+        iteration += 1
+        fun, fun_d1, fun_d2 = _complementary_content(r, x, tailcoverage)
+        # Halley's method
+        r = r - 2 * fun * fun_d1 / (2 * fun_d1 ** 2 - fun * fun_d2)
+        if iteration > max_iterations:
+            break
+        if np.all(np.abs(fun) < tol):
+            break
+    return r
+
+
+def _complementary_content(r, x, tailcoverage):
+    r"""
+    Complementary content function.
+
+    Calculates difference between the complementary content and given
+    the tailcoverage
+
+    fun(r|x,tailcoverage) = 1 - (normcdf(x+r) - normcdf(x-r)) - tailcoverage
+
+    and the first (fun_d1) and the second (fun_d2) derivative of the function
+
+    fun(r|x,tailcoverage)
+
+    """
+    fun = 0.5 * (
+        spec.erfc((x + r) / SQRT_2)
+        + spec.erfc(-(x - r) / SQRT_2)
+    ) - tailcoverage
+    aux1 = np.exp(-0.5 * (x + r)**2)
+    aux2 = np.exp(-0.5 * (x - r)**2)
+    fun_d1 = -(aux1 + aux2) / SQRT_2PI
+    fun_d2 = -((x - r) * aux2 - (x + r) * aux1) / SQRT_2PI
+    return (fun, fun_d1, fun_d2)
+
+
+def _approx_tol_factor_wald_wolfowitz(tailcoverage, tailconfidence, c, nu):
+    r"""
+    Compute the approximate tolerance factor (Wald-Wolfowitz).
+
+    """
+    r = _find_root(np.sqrt(c), tailcoverage)
+    k = r * np.sqrt(nu / stats.chi2.ppf(tailconfidence, nu))
+    return (k, r)
+
+
+def _approx_tol_factor_witkovsky(tailcoverage, tailconfidence, c, m, nu, A, B):
+    r"""
+    Compute the approximate tolerance factor (Witkovsky).
+
+    """
+    val, _ = integ.quad(lambda z: _expect_fun(z, c, tailcoverage, m), A, B)
+    r = np.sqrt(2 * val)
+    k = r * np.sqrt(nu / stats.chi2.ppf(tailconfidence, nu))
+    return (k, r)
+
+
+def _expect_fun(z, c, tailcoverage, m):
+    r"""
+    Expectation function.
+
+    """
+    r = _find_root(np.sqrt(c) * z, tailcoverage)
+    f = r ** 2 * stats.norm.pdf(z)
+    if m > 1:
+        f = f * (m * (1 - (spec.erfc(z / SQRT_2))) ** (m - 1))
+    return f
diff --git a/toleranceinterval/twoside/twoside.py b/toleranceinterval/twoside/twoside.py
index bd48025..b4aa376 100644
--- a/toleranceinterval/twoside/twoside.py
+++ b/toleranceinterval/twoside/twoside.py
@@ -22,30 +22,194 @@
 # SOFTWARE.
 
 import numpy as np
-from scipy.stats import norm, chi2
 from ..checks import numpy_array, assert_2d_sort
+from . import _normal_exact as exact
+from . import _normal_approx as approx
 
 
-def normal(x, p, g):
+def normal_factor(n, p, g, method=None, m=None, nu=None, d2=None,
+                  simultaneous=False, tailprob=False):
     r"""
-    Compute two-side central tolerance interval using the normal distribution.
+    Compute two-sided central tolerance factor using the normal distribution.
 
-    Computes the two-sided tolerance interval (TI) using the normal
-    distribution. This follows the derivation in [1] to calculate the interval
-    as a factor of sample standard deviations away from the sample mean. This
-    is an approximation, and does not represent the exact TI. For exact TI see
-    the eqution 2.3.4 from [2].
+    Computes the tolerance factor k for the two-sided central tolerance
+    interval (TI) to cover a proportion p of the population with confidence g:
 
-    Parameters
+        TI = [Xmean - k * S, Xmean + k * S]
+
+    where Xmean = mean(X), S = std(X), X = [X_1,...,X_n] is a random sample
+    of size n from the distribution N(mu,sig2) with unknown mean mu and
+    variance sig2.
+
+    The tolerance factor k is determined such that the tolerance intervals
+    with confidence g cover at least the coverage fraction
+    of the distribution N(mu,sigma^2), i.e.
+    Prob[ Prob( Xmean - k * S < X < Xmean + k * S ) >= p ] = g,
+    for X ~ N(mu,sig2) which is independent of Xmean and S.
+
+    By default, this function uses an 'exact' method for computing the factor
+    by Gauss-Kronod quadrature as described in the references [1,2,4]. There
+    are also two approximate methods implemented: the 'howe' method as
+    described in [5], and the 'guenther' method as described in [6]. A brief
+    overview of both approximate methods can be found at NIST's website:
+    https://www.itl.nist.gov/div898/handbook/prc/section2/prc263.htm
+
+    Additional optional parameters are available to consider pooled variance
+    studies when m random samples of size n are available. Furthermore,
+    for the 'exact' method, optional parameters are available to
+    consider simultaneous tolerance intervals as described in [7,8].
+    If S is a pooled estimator of sig, based on m random samples of size n,
+    normal_factor computes the tolerance factor k for the two-sided p-content
+    and g-confidence tolerance intervals
+
+        TI = [Xmean_i - k * S, Xmean_i + k * S], for i = 1,...,m
+
+    where Xmean_i = mean(X_i), X_i = [X_i1,...,X_in] is a random sample of
+    size n from the distribution N(mu_i,sig2) with unknown mean mu_i and
+    variance sig2, and S = sqrt(S2), S2 is the pooled estimator of sig2,
+
+        S2 = (1/nu) * sum_i=1:m ( sum_j=1:n (X_ij - Xmean_i)^2 )
+
+    with nu degrees of freedom, nu = m * (n-1). For the 'exact' method, both
+    the simultaneous and non-simultaneous cases can be considered.
+
+    Standard Parameters
+    -------------------
+    n : scalar
+        Sample size
+    p : scalar in the interval [0.0, 1.0]
+        Coverage (or content) probability,
+        Prob( Xmean - k * S < X < Xmean + k * S ) >= p
+    g : scalar in the interval [0.0, 1.0]
+        Confidence probability,
+        Prob[ Prob( Xmean-k*S < X < Xmean+k*S ) >= p ] = g.
+    method : str
+        Method to use for computing the factor. Available methods are 'exact',
+        'howe', and 'guenther'. If None, the default method is 'exact'.
+
+    Optional Parameters for Pooled Variance Studies
+    -----------------------------------------------
+    m : scalar
+        Number of independent random samples (of size n). If None,
+        default value is m = 1.
+    nu : scalar
+        Degrees of freedom for distribution of the (pooled) sample
+        variance S2. If None, default value is nu = m*(n-1).
+
+    Advanced Optional Parameters Available for 'exact' Method
+    ---------------------------------------------------------
+    d2 : scalar
+        Normalizing constant. For computing the factors of the
+        non-simultaneous tolerance limits (xx'*betaHat +/- k * S)
+        for the linear regression y = XX*beta +epsilon, set d2 =
+        xx'*inv(XX'*XX)*xx.
+        Typically, in simple linear regression the estimator S2 has
+        nu = n-2 degrees of freedom. If None, default value is d2 = 1/n.
+    simultaneous : boolean
+        Logical flag for calculating the factor for
+        simultaneous tolerance intervals. If False, normal_factor will
+        calculate the factor for the non-simultaneous tolerance interval.
+        Default value is False.
+    tailprob : boolean
+        Logical flag for representing the input probabilities
+        'p' and 'g'. If True, the input parameters are
+        represented as the tail coverage (i.e. 1 - p) and tail confidence
+        (i.e. 1 - g). This option is useful if the interest is to
+        calculate the tolerance factor for extremely large values
+        of coverage and/or confidence, close to 1, as e.g. coverage = 1 - 1e-18.
+        Default value is False.
+
+    Returns
+    -------
+    float
+        The calculated tolerance factor for the tolerance interval.
+
+    References
     ----------
+    .. [1] Krishnamoorthy K, Mathew T. (2009). Statistical Tolerance Regions:
+        Theory, Applications, and Computation. John Wiley & Sons, Inc.,
+        Hoboken, New Jersey. ISBN: 978-0-470-38026-0, 512 pages.
+    .. [2] Witkovsky V. On the exact two-sided tolerance intervals for
+        univariate normal distribution and linear regression. Austrian
+        Journal of Statistics 43(4), 2014, 279-92. http://
+        ajs.data-analysis.at/index.php/ajs/article/viewFile/vol43-4-6/35
+    .. [3] ISO 16269-6:2014: Statistical interpretation of data - Part 6:
+        Determination of statistical tolerance intervals.
+    .. [4] Janiga I., Garaj I.: Two-sided tolerance limits of normal
+        distributions with unknown means and unknown common variability.
+        MEASUREMENT SCIENCE REVIEW, Volume 3, Section 1, 2003, 75-78.
+    .. [5] Howe, W. G. “Two-Sided Tolerance Limits for Normal Populations,
+        Some Improvements.” Journal of the American Statistical Association,
+        vol. 64, no. 326, [American Statistical Association, Taylor & Francis,
+        Ltd.], 1969, pp. 610–20, https://doi.org/10.2307/2283644.
+    .. [6] Guenther, W. C. (1977). "Sampling Inspection in Statistical Quality
+        Control",, Griffin's Statistical Monographs, Number 37, London.
+    .. [7] Robert W. Mee (1990) Simultaneous Tolerance Intervals for Normal
+        Populations With Common Variance, Technometrics, 32:1, 83-92,
+        DOI: 10.1080/00401706.1990.10484595
+    .. [8] K. Krishnamoorthy & Saptarshi Chakraberty (2022). Construction of
+        simultaneous tolerance intervals for several normal distributions,
+        Journal of Statistical Computation and Simulation, 92:1, 101-114,
+        DOI: 10.1080/00949655.2021.1932885
+
+    """
+    # Handle default method:
+    if method is None:
+        method = 'exact'
+
+    if method == 'exact':
+        k = exact.tolerance_factor(n, p, g, m, nu, d2, simultaneous, tailprob)
+    elif method == 'howe':
+        k = approx.tolerance_factor_howe(n, p, g, m, nu)
+    elif method == 'guenther':
+        k = approx.tolerance_factor_guenther(n, p, g, m, nu)
+    else:
+        raise ValueError(
+            "Invalid method requested. Valid methods are 'exact', 'howe', or "
+            "'guenther'."
+        )
+    return k
+
+
+def normal(x, p, g, method=None, pool_variance=False):
+    r"""
+    Compute two-sided central tolerance interval using the normal distribution.
+
+    Computes the two-sided tolerance interval (TI) to cover a proportion p of
+    the population with confidence g using the normal distribution. This follows
+    the standard approach to calculate the interval as a factor of sample
+    standard deviations away from the sample mean.
+
+        TI = [Xmean - k * S, Xmean + k * S]
+
+    where Xmean = mean(X), S = std(X), X = [X_1,...,X_n] is a random sample
+    of size n from the distribution N(mu,sig2) with unknown mean mu and
+    variance sig2.
+
+    By default, this function uses an 'exact' method for computing the TI
+    by Gauss-Kronod quadrature. There are also two approximate methods
+    implemented: the 'howe' method, and the 'guenther' method. See the
+    documentation for normal_factor for more details on the methods.
+
+    Standard Parameters
+    -------------------
     x : ndarray (1-D, or 2-D)
         Numpy array of samples to compute the tolerance interval. Assumed data
         type is np.float. Shape of (m, n) is assumed for 2-D arrays with m
         number of sets of sample size n.
     p : float
-        Percentile for central TI to estimate.
+        Percentile for central TI to cover.
     g : float
         Confidence level where g > 0. and g < 1.
+    method : str
+        Method to use for computing the TI. Available methods are 'exact',
+        'howe', and 'guenther'. If None, the default method is 'exact'.
+
+    Optional Parameters for Pooled Variance Studies
+    -----------------------------------------------
+    pool_variance : boolean
+        Consider the m random samples to share the same variance such that
+        the degrees of freedom are nu = m*(n-1). Default is False.
 
     Returns
     -------
@@ -56,11 +220,7 @@ def normal(x, p, g):
 
     References
     ----------
-    .. [1] Young, D. S. (2010). tolerance: An R Package for Estimating
-        Tolerance Intervals. Journal of Statistical Software; Vol 1, Issue 5
-        (2010). Retrieved from http://dx.doi.org/10.18637/jss.v036.i05
-    .. [2] Krishnamoorthy, K. and Mathew, T. (2009). Statistical Tolerance
-        Regions: Theory, Applications, and Computation. Wiley, Hoboken, NJ.
+    See the documentation for normal_factor for a complete list of references.
 
     Examples
     --------
@@ -83,24 +243,25 @@ def normal(x, p, g):
     x = numpy_array(x)  # check if numpy array, if not make numpy array
     x = assert_2d_sort(x)
     m, n = x.shape
-    alpha = 1.0 - g
-    zp = norm.ppf((1.0+p)/2.)
-    u = zp * np.sqrt(1 + (1.0/n))
-    chi2v = chi2.ppf(alpha, df=n-1.0)
-    v = np.sqrt((n-1.0)/chi2v)
-    w = np.sqrt(1.0 + ((n-3.0-chi2v)/(2.0*(n+1.0)**2)))
-    k = u * v * w
+
+    # Handle pooled variance case
+    if pool_variance:
+        _m = m
+    else:
+        _m = 1
+
+    k = normal_factor(n, p, g, method, _m)
     bound = np.zeros((m, 2))
     xmu = x.mean(axis=1)
-    kstd = k*x.std(axis=1, ddof=1)
+    kstd = k * x.std(axis=1, ddof=1)
     bound[:, 0] = xmu - kstd
     bound[:, 1] = xmu + kstd
     return bound
 
 
-def lognormal(x, p, g):
+def lognormal(x, p, g, method=None, pool_variance=False):
     r"""
-    Two-side central tolerance interval using the lognormal distribution.
+    Two-sided central tolerance interval using the lognormal distribution.
 
     Computes the two-sided tolerance interval using the lognormal distribution.
     This just performs a ln and exp transformations of the normal distribution.
@@ -115,6 +276,15 @@ def lognormal(x, p, g):
         Percentile for central TI to estimate.
     g : float
         Confidence level where g > 0. and g < 1.
+    method : str
+        Method to use for computing the TI. Available methods are 'exact',
+        'howe', and 'guenther'. If None, the default method is 'exact'.
+
+    Optional Parameters for Pooled Variance Studies
+    -----------------------------------------------
+    pool_variance : boolean
+        Consider the m random samples to share the same variance such that
+        the degrees of freedom are nu = m*(n-1). Default is False.
 
     Returns
     -------
@@ -143,4 +313,4 @@ def lognormal(x, p, g):
     """
     x = numpy_array(x)  # check if numpy array, if not make numpy array
     x = assert_2d_sort(x)
-    return np.exp(normal(np.log(x), p, g))
+    return np.exp(normal(np.log(x), p, g, method, pool_variance))


Test patch:
diff --git a/tests/test_oneside_nonparmetric.py b/tests/test_oneside_nonparmetric.py
index 95768d5..03c6994 100644
--- a/tests/test_oneside_nonparmetric.py
+++ b/tests/test_oneside_nonparmetric.py
@@ -16,44 +16,6 @@ class TestEverything(unittest.TestCase):
                   0.99, 0.99])
     G = np.array([0.90, 0.95, 0.99, 0.90, 0.95, 0.99, 0.90, 0.95, 0.99, 0.90,
                   0.95, 0.99])
-    K = np.array([10.253, 13.090, 18.500, 20.581, 26.260, 37.094, 103.029,
-                  131.426, 185.617, 4.258, 5.311, 7.340, 6.155, 7.656, 10.553,
-                  13.995, 17.370, 23.896, 3.188, 3.957, 5.438, 4.162, 5.144,
-                  7.042, 7.380, 9.083, 12.387, 2.742, 3.400, 4.666, 3.407,
-                  4.203, 5.741, 5.362, 6.578, 8.939, 2.494, 3.092, 4.243,
-                  3.006, 3.708, 5.062, 4.411, 5.406, 7.335, 2.333, 2.894,
-                  3.972, 2.755, 3.399, 4.642, 3.859, 4.728, 6.412, 2.219,
-                  2.754, 3.783, 2.582, 3.187, 4.354, 3.497, 4.285, 5.812,
-                  2.133, 2.650, 3.641, 2.454, 3.031, 4.143, 3.240, 3.972,
-                  5.389, 2.066, 2.568, 3.532, 2.355, 2.911, 3.981, 3.048,
-                  3.738, 5.074, 2.011, 2.503, 3.443, 2.275, 2.815, 3.852,
-                  2.898, 3.556, 4.829, 1.966, 2.448, 3.371, 2.210, 2.736,
-                  3.747, 2.777, 3.410, 4.633, 1.928, 2.402, 3.309, 2.155,
-                  2.671, 3.659, 2.677, 3.290, 4.472, 1.895, 2.363, 3.257,
-                  2.109, 2.614, 3.585, 2.593, 3.189, 4.337, 1.867, 2.329,
-                  3.212, 2.068, 2.566, 3.520, 2.521, 3.102, 4.222, 1.842,
-                  2.299, 3.172, 2.033, 2.524, 3.464, 2.459, 3.028, 4.123,
-                  1.819, 2.272, 3.137, 2.002, 2.486, 3.414, 2.405, 2.963,
-                  4.037, 1.800, 2.249, 3.105, 1.974, 2.453, 3.370, 2.357,
-                  2.905, 3.960, 1.782, 2.227, 3.077, 1.949, 2.423, 3.331,
-                  2.314, 2.854, 3.892, 1.765, 2.208, 3.052, 1.926, 2.396,
-                  3.295, 2.276, 2.808, 3.832, 1.750, 2.190, 3.028, 1.905,
-                  2.371, 3.263, 2.241, 2.766, 3.777, 1.737, 2.174, 3.007,
-                  1.886, 2.349, 3.233, 2.209, 2.729, 3.727, 1.724, 2.159,
-                  2.987, 1.869, 2.328, 3.206, 2.180, 2.694, 3.681, 1.712,
-                  2.145, 2.969, 1.853, 2.309, 3.181, 2.154, 2.662, 3.640,
-                  1.702, 2.132, 2.952, 1.838, 2.292, 3.158, 2.129, 2.633,
-                  3.601, 1.657, 2.080, 2.884, 1.777, 2.220, 3.064, 2.030,
-                  2.515, 3.447, 1.598, 2.010, 2.793, 1.697, 2.125, 2.941,
-                  1.902, 2.364, 3.249, 1.559, 1.965, 2.735, 1.646, 2.065,
-                  2.862, 1.821, 2.269, 3.125, 1.532, 1.933, 2.694, 1.609,
-                  2.022, 2.807, 1.764, 2.202, 3.038, 1.511, 1.909, 2.662,
-                  1.581, 1.990, 2.765, 1.722, 2.153, 2.974, 1.495, 1.890,
-                  2.638, 1.559, 1.964, 2.733, 1.688, 2.114, 2.924, 1.481,
-                  1.874, 2.618, 1.542, 1.944, 2.706, 1.661, 2.082, 2.883,
-                  1.470, 1.861, 2.601, 1.527, 1.927, 2.684, 1.639, 2.056,
-                  2.850])
-
     K = np.array([1, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan,
                   np.nan, np.nan, np.nan, np.nan, 2, 1, np.nan, np.nan, np.nan,
                   np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, 3, 2,
diff --git a/tests/test_twoside_lognormal.py b/tests/test_twoside_lognormal.py
index c75a160..5ac09a8 100644
--- a/tests/test_twoside_lognormal.py
+++ b/tests/test_twoside_lognormal.py
@@ -7,28 +7,53 @@ import unittest
 
 class TestEverything(unittest.TestCase):
 
-    def test_nist_approx(self):
-        # values from:
-        # https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/PASS/Tolerance_Intervals_for_Normal_Data.pdf
-        G = [0.9, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95,
-             0.95]
-        N = [26, 200, 400, 600, 800, 1000, 1200, 1400, 1600, 5910, 866, 179]
-        P = [0.8, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]
-        K = [1.6124, 1.7984, 1.7493, 1.7287, 1.7168, 1.7088, 1.7029, 1.6984,
-             1.6948, 1.6703, 1.7138, 1.8084]
+    def _run_test_cases(self, N, P, G, K, method):
         for i, k in enumerate(K):
             n = N[i]
-            x = np.random.random(n)*10
+            x = np.random.random(n) * 10
             xmu = np.mean(np.log(x))
             xstd = np.std(np.log(x), ddof=1)
             p = P[i]
             g = G[i]
-            bound = lognormal(x, p, g)
-            k_hat_l = (xmu - np.log(bound[0, 0]))/xstd
-            k_hat_u = (np.log(bound[0, 1]) - xmu)/xstd
+            bound = lognormal(x, p, g, method=method)
+            k_hat_l = (xmu - np.log(bound[0, 0])) / xstd
+            k_hat_u = (np.log(bound[0, 1]) - xmu) / xstd
             self.assertTrue(np.isclose(k, k_hat_l, rtol=1e-4, atol=1e-5))
             self.assertTrue(np.isclose(k, k_hat_u, rtol=1e-4, atol=1e-5))
 
+    def test_exact(self):
+        # We test only one case here mostly as a sanity check to make sure
+        # the functions are hooked up correctly. See two-sided tolerance
+        # factor unit tests for exhaustive checking of exact tolerance factors.
+        # Test case is drawn from ISO Table F1.
+        G = [0.90]
+        N = [35]
+        P = [0.90]
+        K = [1.9906]
+        self._run_test_cases(N, P, G, K, 'exact')
+
+    def test_guenther_approx(self):
+        # values from:
+        # https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/PASS/Tolerance_Intervals_for_Normal_Data.pdf
+        G = [0.9, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95,
+             0.95]
+        N = [26, 200, 400, 600, 800, 1000, 1200, 1400, 1600, 5910, 866, 179]
+        P = [0.8, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]
+        K = [1.6124, 1.7984, 1.7493, 1.7287, 1.7168, 1.7088, 1.7029, 1.6984,
+             1.6948, 1.6703, 1.7138, 1.8084]
+        self._run_test_cases(N, P, G, K, 'guenther')
+
+    def test_howe_approx(self):
+        # values from:
+        # https://www.itl.nist.gov/div898/handbook/prc/section2/prc263.htm
+        # Howe's method is implicitly tested to some extent by Guenther's
+        # method.
+        G = [0.99]
+        N = [43]
+        P = [0.9]
+        K = [2.2173]
+        self._run_test_cases(N, P, G, K, 'howe')
+
     def test_random_shapes(self):
         M = [3, 10, 20]
         N = [5, 10, 20]
diff --git a/tests/test_twoside_normal.py b/tests/test_twoside_normal.py
index 2edb795..e401cb0 100644
--- a/tests/test_twoside_normal.py
+++ b/tests/test_twoside_normal.py
@@ -7,80 +7,32 @@ import unittest
 
 class TestEverything(unittest.TestCase):
 
-    # Tabuluar values from Table XII in from:
-    # Montgomery, D. C., & Runger, G. C. (2018). Chapter 8. Statistical
-    # Intervals for a Single Sample. In Applied Statistics and Probability
-    # for Engineers, 7th Edition.
-
-    # Control",, Griffin's Statistical Monographs, Number 37, London.
-    # sample_sizes = np.array([2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
-    #                         16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 30, 40,
-    # 50,
-    #                         60, 70, 80, 90, 100])
-    # P = np.array([0.90, 0.95, 0.99, 0.90, 0.95, 0.99, 0.90, 0.95, 0.99])
-    # G = np.array([0.90, 0.90, 0.90, 0.95, 0.95, 0.95, 0.99, 0.99, 0.99])
-    # K = np.array([15.978, 18.800, 24.167, 32.019, 37.674, 48.430, 160.193,
-    #               188.491, 242.300, 5.847, 6.919, 8.974, 8.380, 9.916,
-    # 12.861,
-    #               18.930, 22.401, 29.055, 4.166, 4.943, 6.440, 5.369, 6.370,
-    #               8.299, 9.398, 11.150, 14.527, 3.949, 4.152, 5.423, 4.275,
-    #               5.079, 6.634, 6.612, 7.855, 10.260, 3.131, 3.723, 4.870,
-    #               3.712, 4.414, 5.775, 5.337, 6.345, 8.301, 2.902, 3.452,
-    #               4.521, 3.369, 4.007, 5.248, 4.613, 5.488, 7.187, 2.743,
-    #               3.264, 4.278, 3.136, 3.732, 4.891, 4.147, 4.936, 6.468,
-    #               2.626, 3.125, 4.098, 2.967, 3.532, 4.631, 3.822, 4.550,
-    #               5.966, 2.535, 3.018, 3.959, 2.839, 3.379, 4.433, 3.582,
-    #               4.265, 5.594, 2.463, 2.933, 3.849, 2.737, 3.259, 4.277,
-    #               3.397, 4.045, 5.308, 2.404, 2.863, 3.758, 2.655, 3.162,
-    #               4.150, 3.250, 3.870, 5.079, 2.355, 2.805, 3.682, 2.587,
-    #               3.081, 4.044, 3.130, 3.727, 4.893, 2.314, 2.756, 3.618,
-    #               2.529, 3.012, 3.955, 3.029, 3.608, 4.737, 2.278, 2.713,
-    #               3.562, 2.480, 2.954, 3.878, 2.945, 3.507, 4.605, 2.246,
-    #               2.676, 3.514, 2.437, 2.903, 3.812, 2.872, 3.421, 4.492,
-    #               2.219, 2.643, 3.471, 2.400, 2.858, 3.754, 2.808, 3.345,
-    #               4.393, 2.194, 2.614, 3.433, 2.366, 2.819, 3.702, 2.753,
-    #               3.279, 4.307, 2.172, 2.588, 3.399, 2.337, 2.784, 3.656,
-    #               2.703, 3.221, 4.230, 2.152, 2.564, 3.368, 2.310, 2.752,
-    #               3.615, 2.659, 3.168, 4.161, 2.135, 2.543, 3.340, 2.286,
-    #               2.723, 3.577, 2.620, 3.121, 4.100, 2.118, 2.524, 3.315,
-    #               2.264, 2.697, 3.543, 2.584, 3.078, 4.044, 2.103, 2.506,
-    #               3.292, 2.244, 2.673, 3.512, 2.551, 3.040, 3.993, 2.089,
-    #               2.489, 3.270, 2.225, 2.651, 3.483, 2.522, 3.004, 3.947,
-    #               2.077, 2.474, 3.251, 2.208, 2.631, 3.457, 2.494, 2.972,
-    #               3.904, 2.025, 2.413, 3.170, 2.140, 2.529, 3.350, 2.385,
-    #               2.841, 3.733, 1.959, 2.334, 3.066, 2.052, 2.445, 3.213,
-    #               2.247, 2.677, 3.518, 1.916, 2.284, 3.001, 1.996, 2.379,
-    #               3.126, 2.162, 2.576, 3.385, 1.887, 2.248, 2.955, 1.958,
-    #               2.333, 3.066, 2.103, 2.506, 3.293, 1.865, 2.222, 2.920,
-    #               1.929, 2.299, 3.021, 2.060, 2.454, 3.225, 1.848, 2.202,
-    #               2.894, 1.907, 2.272, 2.986, 2.026, 2.414, 3.173, 1.834,
-    #               2.185, 2.872, 1.889, 2.251, 2.958, 1.999, 2.382, 3.130,
-    #               1.822, 2.172, 2.854, 1.874, 2.233, 2.934, 1.977, 2.355,
-    #               3.096])
-
-    # K = K.reshape(sample_sizes.size, P.size)
-
-    # def test_montgomery_bounds(self):
-    #     for i, row in enumerate(self.K):
-    #         n = self.sample_sizes[i]
-    #         x = np.random.random(n)*10
-    #         xmu = x.mean()
-    #         xstd = x.std(ddof=1)
+    def _run_test_cases(self, N, P, G, K, method):
+        for i, k in enumerate(K):
+            n = N[i]
+            x = np.random.random(n) * 10
+            xmu = x.mean()
+            xstd = x.std(ddof=1)
+            p = P[i]
+            g = G[i]
+            bound = normal(x, p, g, method=method)
+            k_hat_l = (xmu - bound[0, 0]) / xstd
+            k_hat_u = (bound[0, 1] - xmu) / xstd
+            self.assertTrue(np.isclose(k, k_hat_l, rtol=1e-4, atol=1e-5))
+            self.assertTrue(np.isclose(k, k_hat_u, rtol=1e-4, atol=1e-5))
 
-    #         for j, k in enumerate(row):
-    #             p = self.P[j]
-    #             g = self.G[j]
-    #             chi2v = chi2.ppf(1-g, df=n-1.0)
-    #             w = np.sqrt(1.0 + ((n-3.0-chi2v)/(2.0*(n+1.0)**2)))
-    #             k = k*w
-    #             bound = normal(x, p, g)
-    #             k_hat_l = (xmu - bound[0, 0])/xstd
-    #             k_hat_u = (bound[0, 1] - xmu)/xstd
-    #             print(k, k_hat_l, k_hat_u, bound)
-    #             self.assertTrue(np.isclose(k, k_hat_l, rtol=1e-3, atol=1e-4))
-    #             self.assertTrue(np.isclose(k, k_hat_u, rtol=1e-3, atol=1e-4))
+    def test_exact(self):
+        # We test only one case here mostly as a sanity check to make sure
+        # the functions are hooked up correctly. See two-sided tolerance
+        # factor unit tests for exhaustive checking of exact tolerance factors.
+        # Test case is drawn from ISO Table F1.
+        G = [0.90]
+        N = [35]
+        P = [0.90]
+        K = [1.9906]
+        self._run_test_cases(N, P, G, K, 'exact')
 
-    def test_nist_approx(self):
+    def test_guenther_approx(self):
         # values from:
         # https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/PASS/Tolerance_Intervals_for_Normal_Data.pdf
         G = [0.9, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95, 0.95,
@@ -89,18 +41,18 @@ class TestEverything(unittest.TestCase):
         P = [0.8, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9, 0.9]
         K = [1.6124, 1.7984, 1.7493, 1.7287, 1.7168, 1.7088, 1.7029, 1.6984,
              1.6948, 1.6703, 1.7138, 1.8084]
-        for i, k in enumerate(K):
-            n = N[i]
-            x = np.random.random(n)*10
-            xmu = x.mean()
-            xstd = x.std(ddof=1)
-            p = P[i]
-            g = G[i]
-            bound = normal(x, p, g)
-            k_hat_l = (xmu - bound[0, 0])/xstd
-            k_hat_u = (bound[0, 1] - xmu)/xstd
-            self.assertTrue(np.isclose(k, k_hat_l, rtol=1e-4, atol=1e-5))
-            self.assertTrue(np.isclose(k, k_hat_u, rtol=1e-4, atol=1e-5))
+        self._run_test_cases(N, P, G, K, 'guenther')
+
+    def test_howe_approx(self):
+        # values from:
+        # https://www.itl.nist.gov/div898/handbook/prc/section2/prc263.htm
+        # Howe's method is implicitly tested to some extent by Guenther's
+        # method.
+        G = [0.99]
+        N = [43]
+        P = [0.9]
+        K = [2.2173]
+        self._run_test_cases(N, P, G, K, 'howe')
 
     def test_random_shapes(self):
         M = [3, 10, 20]
diff --git a/tests/test_twoside_normal_factor_iso.py b/tests/test_twoside_normal_factor_iso.py
new file mode 100644
index 0000000..05ee6f7
--- /dev/null
+++ b/tests/test_twoside_normal_factor_iso.py
@@ -0,0 +1,434 @@
+# Author: Copyright (c) 2021 Jed Ludlow
+# License: MIT License
+
+"""
+Test normla_factor against standard tables of tolerance factors
+as published in ISO 16269-6:2014 Annex F.
+
+A sampling of values from the tables is included here for brevity.
+
+"""
+
+import numpy as np
+import toleranceinterval.twoside as ts
+import unittest
+
+
+def decimal_ceil(x, places):
+    """
+    Apply ceiling function at a decimal place.
+
+    The tables of tolerance factors in ISO 16269-6:2014 provide
+    the tolerance factors to a certain number of decimal places. The values
+    at that final decimal place reflect the application of the ceiling function
+    at that decimal place.
+
+    """
+    x *= 10 ** places
+    x = np.ceil(x)
+    x /= 10 ** places
+    return x
+
+
+class BaseTestIso:
+
+    class TestIsoTableF(unittest.TestCase):
+
+        def test_tolerance_factor(self):
+            for row_idx, row in enumerate(self.factor_k5):
+                for col_idx, k5 in enumerate(row):
+                    k = ts.normal_factor(
+                        self.sample_size[row_idx],
+                        self.coverage,
+                        self.confidence,
+                        method='exact',
+                        m=self.number_of_samples[col_idx])
+                    k = decimal_ceil(k, places=4)
+                    self.assertAlmostEqual(k, k5, places=4)
+
+
+class TestIsoF1(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.90
+
+    confidence = 0.90
+
+    # This is n from the table.
+    sample_size = np.array([
+        2, 8, 16, 35, 100, 300, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 2
+        15.5124, 6.0755, 4.5088, 3.8875, 3.5544,
+        3.3461, 3.2032, 3.0989, 3.0193, 2.9565,
+        # n = 8
+        2.7542, 2.3600, 2.2244, 2.1530, 2.1081,
+        2.0769, 2.0539, 2.0361, 2.0220, 2.0104,
+        # n = 16
+        2.2537, 2.0574, 1.9833, 1.9426, 1.9163,
+        1.8977, 1.8837, 1.8727, 1.8638, 1.8564,
+        # n = 35
+        1.9906, 1.8843, 1.8417, 1.8176, 1.8017,
+        1.7902, 1.7815, 1.7747, 1.7690, 1.7643,
+        # n = 100
+        1.8232, 1.7697, 1.7473, 1.7343, 1.7256,
+        1.7193, 1.7144, 1.7105, 1.7073, 1.7047,
+        # n = 300
+        1.7401, 1.7118, 1.6997, 1.6925, 1.6877,
+        1.6842, 1.6815, 1.6793, 1.6775, 1.6760,
+        # n = 1000
+        1.6947, 1.6800, 1.6736, 1.6698, 1.6672,
+        1.6653, 1.6639, 1.6627, 1.6617, 1.6609,
+        # n = infinity
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF2(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.95
+
+    confidence = 0.90
+
+    # This is n from the table.
+    sample_size = np.array([
+        3, 9, 15, 30, 90, 150, 400, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 3
+        6.8233, 4.3320, 3.7087, 3.4207, 3.2528,
+        3.1420, 3.0630, 3.0038, 2.9575, 2.9205,
+        # n = 9
+        3.1323, 2.7216, 2.5773, 2.5006, 2.4521,
+        2.4182, 2.3931, 2.3737, 2.3581, 2.3454,
+        # n = 15
+        2.7196, 2.4718, 2.3789, 2.3280, 2.2951,
+        2.2719, 2.2545, 2.2408, 2.2298, 2.2206,
+        # n = 30
+        2.4166, 2.2749, 2.2187, 2.1870, 2.1662,
+        2.1513, 2.1399, 2.1309, 2.1236, 2.1175,
+        # n = 90
+        2.1862, 2.1182, 2.0898, 2.0733, 2.0624,
+        2.0544, 2.0482, 2.0433, 2.0393, 2.0360,
+        # n = 150
+        2.1276, 2.0775, 2.0563, 2.0439, 2.0356,
+        2.0296, 2.0249, 2.0212, 2.0181, 2.0155,
+        # n = 400
+        2.0569, 2.0282, 2.0158, 2.0085, 2.0035,
+        1.9999, 1.9971, 1.9949, 1.9930, 1.9915,
+        # n = 1000
+        2.0193, 2.0018, 1.9942, 1.9897, 1.9866,
+        1.9844, 1.9826, 1.9812, 1.9800, 1.9791,
+        # n = infinity
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF3(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.99
+
+    confidence = 0.90
+
+    # This is n from the table.
+    sample_size = np.array([
+        4, 8, 17, 28, 100, 300, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 4
+        6.3722, 4.6643, 4.1701, 3.9277, 3.7814,
+        3.6825, 3.6108, 3.5562, 3.5131, 3.4782,
+        # n = 8
+        4.2707, 3.6541, 3.4408, 3.3281, 3.2572,
+        3.2078, 3.1712, 3.1428, 3.1202, 3.1016,
+        # n = 17
+        3.4741, 3.1819, 3.0708, 3.0095, 2.9698,
+        2.9416, 2.9204, 2.9037, 2.8902, 2.8791,
+        # n = 28
+        3.2023, 3.0062, 2.9286, 2.8850, 2.8564,
+        2.8358, 2.8203, 2.8080, 2.7980, 2.7896,
+        # n = 100
+        2.8548, 2.7710, 2.7358, 2.7155, 2.7018,
+        2.6919, 2.6843, 2.6782, 2.6732, 2.6690,
+        # n = 300
+        2.7249, 2.6806, 2.6616, 2.6504, 2.6429,
+        2.6374, 2.6331, 2.6297, 2.6269, 2.6245,
+        # n = 1000
+        2.6538, 2.6308, 2.6208, 2.6148, 2.6108,
+        2.6079, 2.6056, 2.6037, 2.6022, 2.6009,
+        # n = infinity
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF4(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.90
+
+    confidence = 0.95
+
+    # This is n from the table.
+    sample_size = np.array([
+        2, 8, 16, 35, 150, 500, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 2
+        31.0923, 8.7252, 5.8380, 4.7912, 4.2571,
+        3.9341, 3.7179, 3.5630, 3.4468, 3.3565,
+        # n = 8
+        3.1561, 2.5818, 2.3937, 2.2974, 2.2381,
+        2.1978, 2.1685, 2.1463, 2.1289, 2.1149,
+        # n = 16
+        2.4486, 2.1771, 2.0777, 2.0241, 1.9899,
+        1.9661, 1.9483, 1.9346, 1.9237, 1.9147,
+        # n = 35
+        2.0943, 1.9515, 1.8953, 1.8638, 1.8432,
+        1.8285, 1.8174, 1.8087, 1.8016, 1.7957,
+        # n = 150
+        1.8260, 1.7710, 1.7478, 1.7344, 1.7254,
+        1.7188, 1.7137, 1.7097, 1.7064, 1.7036,
+        # n = 500
+        1.7374, 1.7098, 1.6979, 1.6908, 1.6861,
+        1.6826, 1.6799, 1.6777, 1.6760, 1.6744,
+        # n = 1000
+        1.7088, 1.6898, 1.6816, 1.6767, 1.6734,
+        1.6709, 1.6690, 1.6675, 1.6663, 1.6652,
+        # n = infinity
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF5(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.95
+
+    confidence = 0.95
+
+    # This is n from the table.
+    sample_size = np.array([
+        5, 10, 26, 90, 200, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 5
+        5.0769, 3.6939, 3.2936, 3.0986, 2.9820,
+        2.9041, 2.8482, 2.8062, 2.7734, 2.7472,
+        # n = 10
+        3.3935, 2.8700, 2.6904, 2.5964, 2.5377,
+        2.4973, 2.4677, 2.4450, 2.4271, 2.4125,
+        # n = 26
+        2.6188, 2.4051, 2.3227, 2.2771, 2.2476,
+        2.2266, 2.2108, 2.1985, 2.1886, 2.1803,
+        # n = 90
+        2.2519, 2.1622, 2.1251, 2.1037, 2.0895,
+        2.0792, 2.0713, 2.0650, 2.0598, 2.0555,
+        # n = 200
+        2.1430, 2.0877, 2.0642, 2.0505, 2.0413,
+        2.0346, 2.0294, 2.0253, 2.0219, 2.0190,
+        # n = 1000
+        2.0362, 2.0135, 2.0037, 1.9979, 1.9939,
+        1.9910, 1.9888, 1.9870, 1.9855, 1.9842,
+        # n = infinity
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF6(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.99
+
+    confidence = 0.95
+
+    # This is n from the table.
+    sample_size = np.array([
+        3, 9, 17, 35, 100, 500, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 3
+        12.6472, 6.8474, 5.5623, 4.9943, 4.6711,
+        4.4612, 4.3133, 4.2032, 4.1180, 4.0500,
+        # n = 9
+        4.6329, 3.8544, 3.5909, 3.4534, 3.3677,
+        3.3085, 3.2651, 3.2317, 3.2052, 3.1837,
+        # n = 17
+        3.7606, 3.3572, 3.2077, 3.1264, 3.0743,
+        3.0377, 3.0104, 2.9892, 2.9722, 2.9582,
+        # n = 35
+        3.2762, 3.0522, 2.9638, 2.9143, 2.8818,
+        2.8586, 2.8411, 2.8273, 2.8161, 2.8068,
+        # n = 100
+        2.9356, 2.8253, 2.7794, 2.7529, 2.7352,
+        2.7224, 2.7126, 2.7048, 2.6984, 2.6930,
+        # n = 500
+        2.7208, 2.6775, 2.6588, 2.6478, 2.6403,
+        2.6349, 2.6307, 2.6273, 2.6245, 2.6221,
+        # n = infinity
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF7(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.90
+
+    confidence = 0.99
+
+    # This is n from the table.
+    sample_size = np.array([
+        4, 10, 22, 80, 200, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 4
+        9.4162, 4.9212, 3.9582, 3.5449, 3.3166,
+        3.1727, 3.0742, 3.0028, 2.9489, 2.9068,
+        # n = 10
+        3.6167, 2.8193, 2.5709, 2.4481, 2.3748,
+        2.3265, 2.2923, 2.2671, 2.2477, 2.2324,
+        # n = 22
+        2.5979, 2.2631, 2.1429, 2.0791, 2.0393,
+        2.0120, 1.9921, 1.9770, 1.9652, 1.9558,
+        # n = 80
+        2.0282, 1.9056, 1.8562, 1.8281, 1.8097,
+        1.7964, 1.7864, 1.7785, 1.7721, 1.7668,
+        # n = 200
+        1.8657, 1.7973, 1.7686, 1.7520, 1.7409,
+        1.7328, 1.7266, 1.7216, 1.7176, 1.7142,
+        # n = 1000
+        1.7359, 1.7086, 1.6967, 1.6897, 1.6850,
+        1.6815, 1.6788, 1.6767, 1.6749, 1.6734,
+        # n = infinity
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+        1.6449, 1.6449, 1.6449, 1.6449, 1.6449,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF8(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.95
+
+    confidence = 0.99
+
+    # This is n from the table.
+    sample_size = np.array([
+        2, 9, 17, 40, 150, 500, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 2
+        182.7201, 23.1159, 11.9855, 8.7010, 7.1975,
+        6.3481, 5.8059, 5.4311, 5.1573, 4.9489,
+        # n = 9
+        4.5810, 3.4807, 3.1443, 2.9784, 2.8793,
+        2.8136, 2.7670, 2.7324, 2.7057, 2.6846,
+        # n = 17
+        3.3641, 2.8501, 2.6716, 2.5784, 2.5207,
+        2.4814, 2.4529, 2.4314, 2.4147, 2.4013,
+        # n = 40
+        2.6836, 2.4425, 2.3498, 2.2987, 2.2658,
+        2.2427, 2.2254, 2.2120, 2.2013, 2.1926,
+        # n = 150
+        2.2712, 2.1740, 2.1336, 2.1103, 2.0948,
+        2.0835, 2.0749, 2.0681, 2.0625, 2.0578,
+        # n = 500
+        2.1175, 2.0697, 2.0492, 2.0372, 2.0291,
+        2.0231, 2.0185, 2.0149, 2.0118, 2.0093,
+        # n = infinity
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+        1.9600, 1.9600, 1.9600, 1.9600, 1.9600,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
+
+
+class TestIsoF9(BaseTestIso.TestIsoTableF):
+
+    coverage = 0.99
+
+    confidence = 0.99
+
+    # This is n from the table.
+    sample_size = np.array([
+        3, 7, 15, 28, 70, 200, 1000, np.inf,
+    ])
+
+    # This is m from the table.
+    number_of_samples = np.arange(1, 11)
+
+    factor_k5 = np.array([
+        # n = 3
+        28.5857, 10.6204, 7.6599, 6.4888, 5.8628,
+        5.4728, 5.2065, 5.0131, 4.8663, 4.7512,
+        # n = 7
+        7.1908, 5.0656, 4.4559, 4.1605, 3.9847,
+        3.8678, 3.7844, 3.7220, 3.6736, 3.6350,
+        # n = 15
+        4.6212, 3.8478, 3.5825, 3.4441, 3.3581,
+        3.2992, 3.2564, 3.2238, 3.1983, 3.1777,
+        # n = 28
+        3.8042, 3.3792, 3.2209, 3.1350, 3.0801,
+        3.0418, 3.0135, 2.9916, 2.9742, 2.9600,
+        # n = 70
+        3.2284, 3.0179, 2.9334, 2.8857, 2.8544,
+        2.8319, 2.8150, 2.8016, 2.7908, 2.7818,
+        # n = 200
+        2.9215, 2.8144, 2.7695, 2.7434, 2.7260,
+        2.7133, 2.7036, 2.6958, 2.6894, 2.6841,
+        # n = 1000
+        2.7184, 2.6756, 2.6570, 2.6461, 2.6387,
+        2.6332, 2.6290, 2.6257, 2.6229, 2.6205,
+        # n = infinity
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+        2.5759, 2.5759, 2.5759, 2.5759, 2.5759,
+    ])
+
+    factor_k5 = factor_k5.reshape(sample_size.size, number_of_samples.size)
diff --git a/tests/test_twoside_normal_factor_mhe.py b/tests/test_twoside_normal_factor_mhe.py
new file mode 100644
index 0000000..e488e39
--- /dev/null
+++ b/tests/test_twoside_normal_factor_mhe.py
@@ -0,0 +1,317 @@
+# Author: Copyright (c) 2021 Jed Ludlow
+# License: MIT License
+
+"""
+Test normal_factor against tables standard values from Meeker, Hahn, and Escobar.
+
+Meeker, William Q.; Hahn, Gerald J.; Escobar, Luis A.. Statistical
+Intervals: A Guide for Practitioners and Researchers (Wiley Series in
+Probability and Statistics). Wiley. Kindle Edition.
+
+Tables J.5a and J.5b provide tolerance factors for various combinations
+of sample size, coverage, and confidence.
+
+"""
+
+import numpy as np
+import toleranceinterval.twoside as ts
+import unittest
+
+
+class BaseTestMHE:
+
+    class TestMeekerHahnEscobarJ5(unittest.TestCase):
+
+        def test_tolerance_factor(self):
+            for row_idx, row in enumerate(self.factor_g):
+                for col_idx, g in enumerate(row):
+                    k = ts.normal_factor(
+                        self.sample_size[row_idx],
+                        self.coverage[col_idx],
+                        self.confidence[col_idx],
+                        method='exact')
+                    self.assertAlmostEqual(k, g, places=3)
+
+
+class TestMeekerHahnEscobarJ5a(BaseTestMHE.TestMeekerHahnEscobarJ5):
+
+    # Table J.5a from Meeker, William Q.; Hahn, Gerald J.; Escobar, Luis A..
+    # Statistical Intervals: A Guide for Practitioners and Researchers
+    # (Wiley Series in Probability and Statistics) (p. 535). Wiley.
+
+    coverage = np.array([
+        0.5, 0.5, 0.5, 0.5, 0.5,
+        0.7, 0.7, 0.7, 0.7, 0.7,
+        0.8, 0.8, 0.8, 0.8, 0.8,
+    ])
+
+    confidence = np.array([
+        0.5, 0.8, 0.9, 0.95, 0.99,
+        0.5, 0.8, 0.9, 0.95, 0.99,
+        0.5, 0.8, 0.9, 0.95, 0.99,
+    ])
+
+    sample_size = np.array([
+        2, 3, 4, 5, 6, 7, 8, 9, 10,
+        11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+        21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
+        35, 40, 50, 60, 120, 240, 480, np.inf,
+    ])
+
+    factor_g = np.array([
+        # n = 2
+        1.243, 3.369, 6.808, 13.652, 68.316, 1.865, 5.023, 10.142,
+        20.331, 101.732, 2.275, 6.110, 12.333, 24.722, 123.699,
+        # n = 3
+        0.942, 1.700, 2.492, 3.585, 8.122, 1.430, 2.562, 3.747,
+        5.382, 12.181, 1.755, 3.134, 4.577, 6.572, 14.867,
+        # n = 4
+        0.852, 1.335, 1.766, 2.288, 4.028, 1.300, 2.026, 2.673,
+        3.456, 6.073, 1.600, 2.486, 3.276, 4.233, 7.431,
+        # n = 5
+        0.808, 1.173, 1.473, 1.812, 2.824, 1.236, 1.788, 2.239,
+        2.750, 4.274, 1.523, 2.198, 2.750, 3.375, 5.240,
+        # n = 6
+        0.782, 1.081, 1.314, 1.566, 2.270, 1.198, 1.651, 2.003,
+        2.384, 3.446, 1.477, 2.034, 2.464, 2.930, 4.231,
+        # n = 7
+        0.764, 1.021, 1.213, 1.415, 1.954, 1.172, 1.562, 1.853,
+        2.159, 2.973, 1.446, 1.925, 2.282, 2.657, 3.656,
+        # n = 8
+        0.752, 0.979, 1.143, 1.313, 1.750, 1.153, 1.499, 1.749,
+        2.006, 2.668, 1.423, 1.849, 2.156, 2.472, 3.284,
+        # n = 9
+        0.742, 0.947, 1.092, 1.239, 1.608, 1.139, 1.451, 1.672,
+        1.896, 2.455, 1.407, 1.791, 2.062, 2.337, 3.024,
+        # n = 10
+        0.735, 0.923, 1.053, 1.183, 1.503, 1.128, 1.414, 1.613,
+        1.811, 2.297, 1.393, 1.746, 1.990, 2.234, 2.831,
+        # n = 11
+        0.729, 0.903, 1.021, 1.139, 1.422, 1.119, 1.385, 1.566,
+        1.744, 2.175, 1.382, 1.710, 1.932, 2.152, 2.682,
+        # n = 12
+        0.724, 0.886, 0.996, 1.103, 1.357, 1.112, 1.360, 1.527,
+        1.690, 2.078, 1.374, 1.680, 1.885, 2.086, 2.563,
+        # n = 13
+        0.720, 0.873, 0.974, 1.073, 1.305, 1.106, 1.339, 1.495,
+        1.645, 1.999, 1.366, 1.654, 1.846, 2.031, 2.466,
+        # n = 14
+        0.717, 0.861, 0.956, 1.048, 1.261, 1.100, 1.322, 1.467,
+        1.608, 1.933, 1.360, 1.633, 1.812, 1.985, 2.386,
+        # n = 15
+        0.714, 0.851, 0.941, 1.027, 1.224, 1.096, 1.307, 1.444,
+        1.575, 1.877, 1.355, 1.614, 1.783, 1.945, 2.317,
+        # n = 16
+        0.711, 0.842, 0.927, 1.008, 1.193, 1.092, 1.293, 1.423,
+        1.547, 1.829, 1.350, 1.598, 1.758, 1.911, 2.259,
+        # n = 17
+        0.709, 0.835, 0.915, 0.992, 1.165, 1.089, 1.282, 1.405,
+        1.522, 1.788, 1.346, 1.584, 1.736, 1.881, 2.207,
+        # n = 18
+        0.707, 0.828, 0.905, 0.978, 1.141, 1.086, 1.271, 1.389,
+        1.501, 1.751, 1.342, 1.571, 1.717, 1.854, 2.163,
+        # n = 19
+        0.705, 0.822, 0.895, 0.965, 1.120, 1.083, 1.262, 1.375,
+        1.481, 1.719, 1.339, 1.559, 1.699, 1.830, 2.123,
+        # n = 20
+        0.704, 0.816, 0.887, 0.953, 1.101, 1.081, 1.253, 1.362,
+        1.464, 1.690, 1.336, 1.549, 1.683, 1.809, 2.087,
+        # n = 21
+        0.702, 0.811, 0.879, 0.943, 1.084, 1.079, 1.246, 1.350,
+        1.448, 1.664, 1.333, 1.540, 1.669, 1.789, 2.055,
+        # n = 22
+        0.701, 0.806, 0.872, 0.934, 1.068, 1.077, 1.239, 1.340,
+        1.434, 1.640, 1.331, 1.531, 1.656, 1.772, 2.026,
+        # n = 23
+        0.700, 0.802, 0.866, 0.925, 1.054, 1.075, 1.232, 1.330,
+        1.420, 1.619, 1.329, 1.523, 1.644, 1.755, 2.000,
+        # n = 24
+        0.699, 0.798, 0.860, 0.917, 1.042, 1.073, 1.226, 1.321,
+        1.408, 1.599, 1.327, 1.516, 1.633, 1.741, 1.976,
+        # n = 25
+        0.698, 0.795, 0.855, 0.910, 1.030, 1.072, 1.221, 1.313,
+        1.397, 1.581, 1.325, 1.509, 1.623, 1.727, 1.954,
+        # n = 26
+        0.697, 0.791, 0.850, 0.903, 1.019, 1.070, 1.216, 1.305,
+        1.387, 1.565, 1.323, 1.503, 1.613, 1.714, 1.934,
+        # n = 27
+        0.696, 0.788, 0.845, 0.897, 1.009, 1.069, 1.211, 1.298,
+        1.378, 1.550, 1.322, 1.497, 1.604, 1.703, 1.915,
+        # n = 28
+        0.695, 0.785, 0.841, 0.891, 1.000, 1.068, 1.207, 1.291,
+        1.369, 1.535, 1.320, 1.492, 1.596, 1.692, 1.898,
+        # n = 29
+        0.694, 0.783, 0.837, 0.886, 0.991, 1.067, 1.202, 1.285,
+        1.360, 1.522, 1.319, 1.487, 1.589, 1.682, 1.882,
+        # n = 30
+        0.694, 0.780, 0.833, 0.881, 0.983, 1.066, 1.199, 1.279,
+        1.353, 1.510, 1.318, 1.482, 1.581, 1.672, 1.866,
+        # n = 35
+        0.691, 0.770, 0.817, 0.859, 0.950, 1.061, 1.182, 1.255,
+        1.320, 1.459, 1.312, 1.462, 1.551, 1.632, 1.803,
+        # n = 40
+        0.689, 0.761, 0.805, 0.843, 0.925, 1.058, 1.170, 1.236,
+        1.296, 1.420, 1.308, 1.446, 1.528, 1.602, 1.756,
+        # n = 50
+        0.686, 0.750, 0.787, 0.820, 0.889, 1.054, 1.152, 1.209,
+        1.260, 1.365, 1.303, 1.424, 1.495, 1.558, 1.688,
+        # n = 60
+        0.684, 0.741, 0.775, 0.804, 0.864, 1.051, 1.139, 1.190,
+        1.235, 1.327, 1.299, 1.408, 1.471, 1.527, 1.641,
+        # n = 120
+        0.679, 0.718, 0.740, 0.759, 0.797, 1.044, 1.104, 1.137,
+        1.166, 1.225, 1.290, 1.365, 1.406, 1.442, 1.514,
+        # n = 240
+        0.677, 0.704, 0.719, 0.731, 0.756, 1.040, 1.082, 1.104,
+        1.124, 1.162, 1.286, 1.337, 1.365, 1.390, 1.437,
+        # n = 480
+        0.676, 0.694, 0.705, 0.713, 0.730, 1.038, 1.067, 1.083,
+        1.096, 1.122, 1.284, 1.320, 1.339, 1.355, 1.387,
+        # n = infinity
+        0.674, 0.674, 0.674, 0.674, 0.674, 1.036, 1.036, 1.036,
+        1.036, 1.036, 1.282, 1.282, 1.282, 1.282, 1.282,
+    ])
+
+    factor_g = factor_g.reshape(sample_size.size, coverage.size)
+
+
+class TestMeekerHahnEscobarJ5b(BaseTestMHE.TestMeekerHahnEscobarJ5):
+
+    # Table J.5b from Meeker, William Q.; Hahn, Gerald J.; Escobar, Luis A..
+    # Statistical Intervals: A Guide for Practitioners and Researchers
+    # (Wiley Series in Probability and Statistics) (p. 535). Wiley.
+
+    coverage = np.array([
+        0.90, 0.90, 0.90, 0.90, 0.90,
+        0.95, 0.95, 0.95, 0.95, 0.95,
+        0.99, 0.99, 0.99, 0.99, 0.99,
+    ])
+
+    confidence = np.array([
+        0.5, 0.8, 0.9, 0.95, 0.99,
+        0.5, 0.8, 0.9, 0.95, 0.99,
+        0.5, 0.8, 0.9, 0.95, 0.99,
+    ])
+
+    sample_size = np.array([
+        2, 3, 4, 5, 6, 7, 8, 9, 10,
+        11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
+        21, 22, 23, 24, 25, 26, 27, 28, 29, 30,
+        35, 40, 50, 60, 120, 240, 480, np.inf,
+    ])
+
+    factor_g = np.array([
+        # n = 2
+        2.869, 7.688, 15.512, 31.092, 155.569, 3.376, 9.032, 18.221,
+        36.519, 182.720, 4.348, 11.613, 23.423, 46.944, 234.877,
+        # n = 3
+        2.229, 3.967, 5.788, 8.306, 18.782, 2.634, 4.679, 6.823, 9.789,
+        22.131, 3.415, 6.051, 8.819, 12.647, 28.586,
+        # n = 4
+        2.039, 3.159, 4.157, 5.368, 9.416, 2.416, 3.736, 4.913, 6.341,
+        11.118, 3.144, 4.850, 6.372, 8.221, 14.405,
+        # n = 5
+        1.945, 2.801, 3.499, 4.291, 6.655, 2.308, 3.318, 4.142, 5.077,
+        7.870, 3.010, 4.318, 5.387, 6.598, 10.220,
+        # n = 6
+        1.888, 2.595, 3.141, 3.733, 5.383, 2.243, 3.078, 3.723, 4.422,
+        6.373, 2.930, 4.013, 4.850, 5.758, 8.292,
+        # n = 7
+        1.850, 2.460, 2.913, 3.390, 4.658, 2.199, 2.920, 3.456, 4.020,
+        5.520, 2.876, 3.813, 4.508, 5.241, 7.191,
+        # n = 8
+        1.823, 2.364, 2.754, 3.156, 4.189, 2.167, 2.808, 3.270, 3.746,
+        4.968, 2.836, 3.670, 4.271, 4.889, 6.479,
+        # n = 9
+        1.802, 2.292, 2.637, 2.986, 3.860, 2.143, 2.723, 3.132, 3.546,
+        4.581, 2.806, 3.562, 4.094, 4.633, 5.980,
+        # n = 10
+        1.785, 2.235, 2.546, 2.856, 3.617, 2.124, 2.657, 3.026, 3.393,
+        4.294, 2.783, 3.478, 3.958, 4.437, 5.610,
+        # n = 11
+        1.772, 2.189, 2.473, 2.754, 3.429, 2.109, 2.604, 2.941, 3.273,
+        4.073, 2.764, 3.410, 3.849, 4.282, 5.324,
+        # n = 12
+        1.761, 2.152, 2.414, 2.670, 3.279, 2.096, 2.560, 2.871, 3.175,
+        3.896, 2.748, 3.353, 3.759, 4.156, 5.096,
+        # n = 13
+        1.752, 2.120, 2.364, 2.601, 3.156, 2.085, 2.522, 2.812, 3.093,
+        3.751, 2.735, 3.306, 3.684, 4.051, 4.909,
+        # n = 14
+        1.744, 2.093, 2.322, 2.542, 3.054, 2.076, 2.490, 2.762, 3.024,
+        3.631, 2.723, 3.265, 3.620, 3.962, 4.753,
+        # n = 15
+        1.737, 2.069, 2.285, 2.492, 2.967, 2.068, 2.463, 2.720, 2.965,
+        3.529, 2.714, 3.229, 3.565, 3.885, 4.621,
+        # n = 16
+        1.731, 2.049, 2.254, 2.449, 2.893, 2.061, 2.439, 2.682, 2.913,
+        3.441, 2.705, 3.198, 3.517, 3.819, 4.507,
+        # n = 17
+        1.726, 2.030, 2.226, 2.410, 2.828, 2.055, 2.417, 2.649, 2.868,
+        3.364, 2.698, 3.171, 3.474, 3.761, 4.408,
+        # n = 18
+        1.721, 2.014, 2.201, 2.376, 2.771, 2.050, 2.398, 2.620, 2.828,
+        3.297, 2.691, 3.146, 3.436, 3.709, 4.321,
+        # n = 19
+        1.717, 2.000, 2.178, 2.346, 2.720, 2.045, 2.381, 2.593, 2.793,
+        3.237, 2.685, 3.124, 3.402, 3.663, 4.244,
+        # n = 20
+        1.714, 1.987, 2.158, 2.319, 2.675, 2.041, 2.365, 2.570, 2.760,
+        3.184, 2.680, 3.104, 3.372, 3.621, 4.175,
+        # n = 21
+        1.710, 1.975, 2.140, 2.294, 2.635, 2.037, 2.351, 2.548, 2.731,
+        3.136, 2.675, 3.086, 3.344, 3.583, 4.113,
+        # n = 22
+        1.707, 1.964, 2.123, 2.272, 2.598, 2.034, 2.338, 2.528, 2.705,
+        3.092, 2.670, 3.070, 3.318, 3.549, 4.056,
+        # n = 23
+        1.705, 1.954, 2.108, 2.251, 2.564, 2.030, 2.327, 2.510, 2.681,
+        3.053, 2.666, 3.054, 3.295, 3.518, 4.005,
+        # n = 24
+        1.702, 1.944, 2.094, 2.232, 2.534, 2.027, 2.316, 2.494, 2.658,
+        3.017, 2.662, 3.040, 3.274, 3.489, 3.958,
+        # n = 25
+        1.700, 1.936, 2.081, 2.215, 2.506, 2.025, 2.306, 2.479, 2.638,
+        2.984, 2.659, 3.027, 3.254, 3.462, 3.915,
+        # n = 26
+        1.698, 1.928, 2.069, 2.199, 2.480, 2.022, 2.296, 2.464, 2.619,
+        2.953, 2.656, 3.015, 3.235, 3.437, 3.875,
+        # n = 27
+        1.696, 1.921, 2.058, 2.184, 2.456, 2.020, 2.288, 2.451, 2.601,
+        2.925, 2.653, 3.004, 3.218, 3.415, 3.838,
+        # n = 28
+        1.694, 1.914, 2.048, 2.170, 2.434, 2.018, 2.279, 2.439, 2.585,
+        2.898, 2.650, 2.993, 3.202, 3.393, 3.804,
+        # n = 29
+        1.692, 1.907, 2.038, 2.157, 2.413, 2.016, 2.272, 2.427, 2.569,
+        2.874, 2.648, 2.983, 3.187, 3.373, 3.772,
+        # n = 30
+        1.691, 1.901, 2.029, 2.145, 2.394, 2.014, 2.265, 2.417, 2.555,
+        2.851, 2.645, 2.974, 3.173, 3.355, 3.742,
+        # n = 35
+        1.684, 1.876, 1.991, 2.094, 2.314, 2.006, 2.234, 2.371, 2.495,
+        2.756, 2.636, 2.935, 3.114, 3.276, 3.618,
+        # n = 40
+        1.679, 1.856, 1.961, 2.055, 2.253, 2.001, 2.211, 2.336, 2.448,
+        2.684, 2.628, 2.905, 3.069, 3.216, 3.524,
+        # n = 50
+        1.672, 1.827, 1.918, 1.999, 2.166, 1.992, 2.177, 2.285, 2.382,
+        2.580, 2.618, 2.861, 3.003, 3.129, 3.390,
+        # n = 60
+        1.668, 1.807, 1.888, 1.960, 2.106, 1.987, 2.154, 2.250, 2.335,
+        2.509, 2.611, 2.830, 2.956, 3.068, 3.297,
+        # n = 120
+        1.656, 1.752, 1.805, 1.851, 1.943, 1.974, 2.087, 2.151, 2.206,
+        2.315, 2.594, 2.743, 2.826, 2.899, 3.043,
+        # n = 240
+        1.651, 1.716, 1.753, 1.783, 1.844, 1.967, 2.045, 2.088, 2.125,
+        2.197, 2.585, 2.688, 2.744, 2.793, 2.887,
+        # n = 480
+        1.648, 1.694, 1.718, 1.739, 1.780, 1.963, 2.018, 2.048, 2.073,
+        2.121, 2.580, 2.652, 2.691, 2.724, 2.787,
+        # n = infinity
+        1.645, 1.645, 1.645, 1.645, 1.645, 1.960, 1.960, 1.960, 1.960,
+        1.960, 2.576, 2.576, 2.576, 2.576, 2.576,
+    ])
+
+    factor_g = factor_g.reshape(sample_size.size, coverage.size)