Repository: cjekel/tolerance_interval_py (language: python)

Issue:
code exact solution for two sided normal distribution

code exact solution for two sided normal distribution



Discussion / hints:
@cjekel, I just released a Python module that computes the exact tolerance factors for two-sided normal tolerance intervals:

https://github.com/jedludlow/norm-tol-int

It could use a home inside of a larger package for computing tolerance intervals. Are you interested in perhaps folding this into your package?
Your python module looks nice!

It's fine if it replaces the current `twoside.normal` function, and the tests you are using are better than mine.

The only thing that will take a bit of work is  some tweaks to match the api of function inputs and outputs.
- Taking in a 1d or 2d array of samples (as oppose to n)
- Returning the left and right bound for each sample set (as oppose to just the factor)

If you're willing to match this style, I'd happily accept a PR and maintain the code! Thanks for reaching out.

Yes, I could adapt the code to conform to your standard `(x, p, g)` interface. That shouldn't be too difficult.

A couple of thoughts to consider:

* There are various approximate methods for the two-sided normal tolerance interval. It might be useful from a comparative perspective to have them represented here instead of just outright replacing the current `normal` with the exact solution. You could simply append something to each function name, e.g., `normal_exact` or `normal_guenther` or `normal_howe`. Or you could go a step further and use subpackages to namespace things, e.g., `toleranceinterval.twoside.normal.exact(x, p, g)` or `toleranceinterval.twoside.normal.howe(x, p, g)`, but that might be overkill. In the former case, `toleranceinterval.twoside.normal` could simply be an alias to `normal_exact`. Thoughts?
* It's actually really useful to have some way to access `tolerance_factor` directly since in some cases you may only have `x_bar`, `s`, and `n` provided. It is also helpful in the planning stages of a study, too. Furthermore, almost all published results for comparison in unit tests would take the form of tolerance factors. Would you be open to exposing `tolerance_factor_exact` or, say, `tolerance_factor_guenther` directly on the API for the various two-sided normal interval methods?
I should add, I'd be willing to address the above ideas in the PR I develop. Just trying to get a feel for your thoughts before getting in too deep.
> There are various approximate methods for the two-sided normal tolerance interval. It might be useful from a comparative perspective to have them represented here instead of just outright replacing the current normal with the exact solution. You could simply append something to each function name, e.g., normal_exact or normal_guenther or normal_howe. Or you could go a step further and use subpackages to namespace things, e.g., toleranceinterval.twoside.normal.exact(x, p, g) or toleranceinterval.twoside.normal.howe(x, p, g), but that might be overkill. In the former case, toleranceinterval.twoside.normal could simply be an alias to normal_exact. Thoughts?

how about this api
```python
toleranceinterval.twoside.normal(x, p, g, method='exact')
toleranceinterval.twoside.normal(x, p, g, method='howe')
toleranceinterval.twoside.normal(x, p, g, method='guenther')
```

> It's actually really useful to have some way to access tolerance_factor directly since in some cases you may only have x_bar, s, and n provided. It is also helpful in the planning stages of a study, too. Furthermore, almost all published results for comparison in unit tests would take the form of tolerance factors. Would you be open to exposing tolerance_factor_exact or, say, tolerance_factor_guenther directly on the API for the various two-sided normal interval methods?

Yes, let's make a factors module somewhere. Not sure about this as it's getting long. I want it to be intuitive to explore the api.  `toleranceinterval.twoside.factor` which can contain `normal_exact, normal_guenther,...`?

--- 

What are you going to do with `tolerance_factor` [def tolerance_factor(n, coverage, confidence, m=None, nu=None, d2=None](https://github.com/jedludlow/norm-tol-int/blob/d08e48a4c6b266f0bcdf5ed62b918226b87a9298/normtolint.py#L57) I think this is up to you.

I'd be fine if you made a small library on PyPi with your project, and then we add your project as a dependency to this package.

I'd also be fine to accept your code into this project (as I said, it looks well done!). 

I like your suggestion of including a keyword argument to set the method. I'll head that direction.

As for the tolerance factors, the normal distribution is kind of unique in that the width of the tolerance interval is just some multiplicative factor of the standard deviation. I that light, maybe it doesn't need a full module. What if we took a page from your suggestion and did this for the full interval API:

```
toleranceinterval.twoside.normal(x, p, g, method='exact')
toleranceinterval.twoside.normal(x, p, g, method='howe')
toleranceinterval.twoside.normal(x, p, g, method='guenther')
```

and then this for the factors API:

```
toleranceinterval.twoside.normal_factor(n, p, g, method='exact')
toleranceinterval.twoside.normal_factor(n, p, g, method='howe')
toleranceinterval.twoside.normal_factor(n, p, g, method='guenther')
```

> I'd be fine if you made a small library on PyPi with your project, and then we add your project as a dependency to this package.
> 
> I'd also be fine to accept your code into this project (as I said, it looks well done!).

At this point I would lean toward just integrating the code into your package. I'd also be happy to help maintain it if some issue arises. That feels like the cleanest solution. I may keep my original repo around just for historical purposes since it represents a very direct port of the MATLAB code in the way it is structured, but I'd edit my README to direct users here for the most up-to-date implementation. 
This all sounds good to me!

On Sat, Dec 4, 2021, 19:56 Jed Ludlow ***@***.***> wrote:

> I like your suggestion of including a keyword argument to set the method.
> I'll head that direction.
>
> As for the tolerance factors, the normal distribution is kind of unique in
> that the width of the tolerance interval is just some multiplicative factor
> of the standard deviation. I that light, maybe it doesn't need a full
> module. What if we took a page from your suggestion and did this for the
> full interval API:
>
> toleranceinterval.twoside.normal(x, p, g, method='exact')
> toleranceinterval.twoside.normal(x, p, g, method='howe')
> toleranceinterval.twoside.normal(x, p, g, method='guenther')
>
> and then this for the factors API:
>
> toleranceinterval.twoside.normal_factor(n, p, g, method='exact')
> toleranceinterval.twoside.normal_factor(n, p, g, method='howe')
> toleranceinterval.twoside.normal_factor(n, p, g, method='guenther')
>
> I'd be fine if you made a small library on PyPi with your project, and
> then we add your project as a dependency to this package.
>
> I'd also be fine to accept your code into this project (as I said, it
> looks well done!).
>
> At this point I would lean toward just integrating the code into your
> package. I'd also be happy to help maintain it if some issue arises. That
> feels like the cleanest solution. I may keep my original repo around just
> for historical purposes since it represents a very direct port of the
> MATLAB code in the way it is structured, but I'd edit my README to direct
> users here for the most up-to-date implementation.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/cjekel/tolerance_interval_py/issues/1#issuecomment-986161835>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ADJ5Z4MYH4K26NUDVW4MKZ3UPLPFHANCNFSM4JKINXPA>
> .
> Triage notifications on the go with GitHub Mobile for iOS
> <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675>
> or Android
> <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.
>
>

@cjekel, I just released a Python module that computes the exact tolerance factors for two-sided normal tolerance intervals:

https://github.com/jedludlow/norm-tol-int

It could use a home inside of a larger package for computing tolerance intervals. Are you interested in perhaps folding this into your package?
Your python module looks nice!

It's fine if it replaces the current `twoside.normal` function, and the tests you are using are better than mine.

The only thing that will take a bit of work is  some tweaks to match the api of function inputs and outputs.
- Taking in a 1d or 2d array of samples (as oppose to n)
- Returning the left and right bound for each sample set (as oppose to just the factor)

If you're willing to match this style, I'd happily accept a PR and maintain the code! Thanks for reaching out.

Yes, I could adapt the code to conform to your standard `(x, p, g)` interface. That shouldn't be too difficult.

A couple of thoughts to consider:

* There are various approximate methods for the two-sided normal tolerance interval. It might be useful from a comparative perspective to have them represented here instead of just outright replacing the current `normal` with the exact solution. You could simply append something to each function name, e.g., `normal_exact` or `normal_guenther` or `normal_howe`. Or you could go a step further and use subpackages to namespace things, e.g., `toleranceinterval.twoside.normal.exact(x, p, g)` or `toleranceinterval.twoside.normal.howe(x, p, g)`, but that might be overkill. In the former case, `toleranceinterval.twoside.normal` could simply be an alias to `normal_exact`. Thoughts?
* It's actually really useful to have some way to access `tolerance_factor` directly since in some cases you may only have `x_bar`, `s`, and `n` provided. It is also helpful in the planning stages of a study, too. Furthermore, almost all published results for comparison in unit tests would take the form of tolerance factors. Would you be open to exposing `tolerance_factor_exact` or, say, `tolerance_factor_guenther` directly on the API for the various two-sided normal interval methods?
I should add, I'd be willing to address the above ideas in the PR I develop. Just trying to get a feel for your thoughts before getting in too deep.
> There are various approximate methods for the two-sided normal tolerance interval. It might be useful from a comparative perspective to have them represented here instead of just outright replacing the current normal with the exact solution. You could simply append something to each function name, e.g., normal_exact or normal_guenther or normal_howe. Or you could go a step further and use subpackages to namespace things, e.g., toleranceinterval.twoside.normal.exact(x, p, g) or toleranceinterval.twoside.normal.howe(x, p, g), but that might be overkill. In the former case, toleranceinterval.twoside.normal could simply be an alias to normal_exact. Thoughts?

how about this api
```python
toleranceinterval.twoside.normal(x, p, g, method='exact')
toleranceinterval.twoside.normal(x, p, g, method='howe')
toleranceinterval.twoside.normal(x, p, g, method='guenther')
```

> It's actually really useful to have some way to access tolerance_factor directly since in some cases you may only have x_bar, s, and n provided. It is also helpful in the planning stages of a study, too. Furthermore, almost all published results for comparison in unit tests would take the form of tolerance factors. Would you be open to exposing tolerance_factor_exact or, say, tolerance_factor_guenther directly on the API for the various two-sided normal interval methods?

Yes, let's make a factors module somewhere. Not sure about this as it's getting long. I want it to be intuitive to explore the api.  `toleranceinterval.twoside.factor` which can contain `normal_exact, normal_guenther,...`?

--- 

What are you going to do with `tolerance_factor` [def tolerance_factor(n, coverage, confidence, m=None, nu=None, d2=None](https://github.com/jedludlow/norm-tol-int/blob/d08e48a4c6b266f0bcdf5ed62b918226b87a9298/normtolint.py#L57) I think this is up to you.

I'd be fine if you made a small library on PyPi with your project, and then we add your project as a dependency to this package.

I'd also be fine to accept your code into this project (as I said, it looks well done!). 

I like your suggestion of including a keyword argument to set the method. I'll head that direction.

As for the tolerance factors, the normal distribution is kind of unique in that the width of the tolerance interval is just some multiplicative factor of the standard deviation. I that light, maybe it doesn't need a full module. What if we took a page from your suggestion and did this for the full interval API:

```
toleranceinterval.twoside.normal(x, p, g, method='exact')
toleranceinterval.twoside.normal(x, p, g, method='howe')
toleranceinterval.twoside.normal(x, p, g, method='guenther')
```

and then this for the factors API:

```
toleranceinterval.twoside.normal_factor(n, p, g, method='exact')
toleranceinterval.twoside.normal_factor(n, p, g, method='howe')
toleranceinterval.twoside.normal_factor(n, p, g, method='guenther')
```

> I'd be fine if you made a small library on PyPi with your project, and then we add your project as a dependency to this package.
> 
> I'd also be fine to accept your code into this project (as I said, it looks well done!).

At this point I would lean toward just integrating the code into your package. I'd also be happy to help maintain it if some issue arises. That feels like the cleanest solution. I may keep my original repo around just for historical purposes since it represents a very direct port of the MATLAB code in the way it is structured, but I'd edit my README to direct users here for the most up-to-date implementation. 
This all sounds good to me!

On Sat, Dec 4, 2021, 19:56 Jed Ludlow ***@***.***> wrote:

> I like your suggestion of including a keyword argument to set the method.
> I'll head that direction.
>
> As for the tolerance factors, the normal distribution is kind of unique in
> that the width of the tolerance interval is just some multiplicative factor
> of the standard deviation. I that light, maybe it doesn't need a full
> module. What if we took a page from your suggestion and did this for the
> full interval API:
>
> toleranceinterval.twoside.normal(x, p, g, method='exact')
> toleranceinterval.twoside.normal(x, p, g, method='howe')
> toleranceinterval.twoside.normal(x, p, g, method='guenther')
>
> and then this for the factors API:
>
> toleranceinterval.twoside.normal_factor(n, p, g, method='exact')
> toleranceinterval.twoside.normal_factor(n, p, g, method='howe')
> toleranceinterval.twoside.normal_factor(n, p, g, method='guenther')
>
> I'd be fine if you made a small library on PyPi with your project, and
> then we add your project as a dependency to this package.
>
> I'd also be fine to accept your code into this project (as I said, it
> looks well done!).
>
> At this point I would lean toward just integrating the code into your
> package. I'd also be happy to help maintain it if some issue arises. That
> feels like the cleanest solution. I may keep my original repo around just
> for historical purposes since it represents a very direct port of the
> MATLAB code in the way it is structured, but I'd edit my README to direct
> users here for the most up-to-date implementation.
>
> —
> You are receiving this because you were mentioned.
> Reply to this email directly, view it on GitHub
> <https://github.com/cjekel/tolerance_interval_py/issues/1#issuecomment-986161835>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ADJ5Z4MYH4K26NUDVW4MKZ3UPLPFHANCNFSM4JKINXPA>
> .
> Triage notifications on the go with GitHub Mobile for iOS
> <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675>
> or Android
> <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.
>
>


Fix 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 @@
 
 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 @@
 
 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 @@ def test_nist_approx(self):
         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)