diff --git a/README.md b/README.md
index 77d8f4b..f2f031e 100644
--- a/README.md
+++ b/README.md
@@ -113,7 +113,7 @@ $ pip install pandas_ta
 
 Latest Version
 --------------
-Best choice! Version: *0.3.27b*
+Best choice! Version: *0.3.28b*
 * Includes all fixes and updates between **pypi** and what is covered in this README.
 ```sh
 $ pip install -U git+https://github.com/twopirllc/pandas-ta
diff --git a/pandas_ta/utils/__init__.py b/pandas_ta/utils/__init__.py
index 34b285f..eb26d1f 100644
--- a/pandas_ta/utils/__init__.py
+++ b/pandas_ta/utils/__init__.py
@@ -2,7 +2,8 @@
 from ._candles import *
 from ._core import *
 from ._math import *
-from ._signals import *
-from ._time import *
 from ._metrics import *
+from ._signals import *
+from ._stats import *
+from ._time import *
 from .data import *
diff --git a/pandas_ta/utils/_math.py b/pandas_ta/utils/_math.py
index 56e0dd7..41c0321 100644
--- a/pandas_ta/utils/_math.py
+++ b/pandas_ta/utils/_math.py
@@ -110,7 +110,7 @@ def geometric_mean(series: Series) -> float:
 
 
 def hpoly(array: npArray, x: Tuple[int, float]) ->  float:
-    """Horner Calculation for Polynomial Evaluation
+    """Horner Calculation for Polynomial Evaluation (hpoly)
 
     array: np.array of polynomial coefficients
         * Convert list or Series to np.array prior to calling the method for
@@ -194,6 +194,7 @@ def strided_window(array, length):
     """as_strided
     creates a view into the array given the exact strides and shape.
     * Recommended to avoid when possible.
+
     Source: https://numpy.org/devdocs/reference/generated/numpy.lib.stride_tricks.as_strided.html
     Pandas TA Issue: https://github.com/twopirllc/pandas-ta/issues/285
     """
diff --git a/pandas_ta/utils/_stats.py b/pandas_ta/utils/_stats.py
new file mode 100644
index 0000000..067a333
--- /dev/null
+++ b/pandas_ta/utils/_stats.py
@@ -0,0 +1,112 @@
+# -*- coding: utf-8 -*-
+from typing import Tuple
+
+from numpy import array as npArray
+from numpy import infty as npInfty
+from numpy import log as npLog
+from numpy import nan as npNaN
+from numpy import pi as npPi
+from numpy import sqrt as npSqrt
+
+from pandas_ta import Imports
+from ._math import hpoly
+
+
+
+def inv_norm(x0: Tuple[float, int]) -> Tuple[float, None]:
+    """Inverse Normal (inv_norm)
+    Calculates the 'x' in which the area under the Gaussian PDF is equal to x0.
+
+    If the user has package "statsmodels" installed, the method will call and
+    return norm().ppf(x0)
+
+
+    Source: https://github.com/scipy/scipy/blob/701ffcc8a6f04509d115aac5e5681c538b5265a2/scipy/special/cephes/ndtri.c
+    """
+
+    if Imports["statsmodels"]:
+        from scipy.stats import norm
+        return norm().ppf(x0)
+
+    # Polynomial Coefficients
+    p0 = npArray([
+        -5.99633501014107895267E1, 9.80010754185999661536E1,
+        -5.66762857469070293439E1, 1.39312609387279679503E1,
+        -1.23916583867381258016E0
+    ])
+    q0 = npArray([
+        1.00000000000000000000E0, 1.95448858338141759834E0,
+        4.67627912898881538453E0, 8.63602421390890590575E1,
+        -2.25462687854119370527E2, 2.00260212380060660359E2,
+        -8.20372256168333339912E1, 1.59056225126211695515E1,
+        -1.18331621121330003142E0
+    ])
+
+    p1 = npArray([
+        4.05544892305962419923E0, 3.15251094599893866154E1,
+        5.71628192246421288162E1, 4.40805073893200834700E1,
+        1.46849561928858024014E1, 2.18663306850790267539E0,
+        -1.40256079171354495875E-1, -3.50424626827848203418E-2,
+        -8.57456785154685413611E-4
+    ])
+    q1 = npArray([
+        1.00000000000000000000E0, 1.57799883256466749731E1,
+        4.53907635128879210584E1, 4.13172038254672030440E1,
+        1.50425385692907503408E1, 2.50464946208309415979E0,
+        -1.42182922854787788574E-1, -3.80806407691578277194E-2,
+        -9.33259480895457427372E-4
+    ])
+
+    p2 = npArray([
+        3.23774891776946035970E0, 6.91522889068984211695E0,
+        3.93881025292474443415E0, 1.33303460815807542389E0,
+        2.01485389549179081538E-1, 1.23716634817820021358E-2,
+        3.01581553508235416007E-4, 2.65806974686737550832E-6,
+        6.23974539184983293730E-9
+    ])
+    q2 = npArray([
+        1.00000000000000000000E0, 6.02427039364742014255E0,
+        3.67983563856160859403E0, 1.37702099489081330271E0,
+        2.16236993594496635890E-1, 1.34204006088543189037E-2,
+        3.28014464682127739104E-4, 2.89247864745380683936E-6,
+        6.79019408009981274425E-9
+    ])
+
+    if x0 == 0.0: return -npInfty
+    if x0 == 1.0: return npInfty
+    if x0 < 0.0 or x0 > 1.0: return npNaN
+
+    negate = True
+    x = x0
+
+    sqrt2pi = npSqrt(2 * npPi)
+    threshold = 0.13533528323661269189
+    if x > 1.0 - threshold:
+        x, negate = 1.0 - x, False
+
+    # 0 <= |x - 0.5| <= 3/8
+    if x > threshold:
+        x -= 0.5
+        x2 = x * x
+        y = x + x * (x2 * hpoly(p0, x2) / hpoly(q0, x2))
+        y *= sqrt2pi
+        return y
+
+    y = npSqrt(-2.0 * npLog(x))
+    y0 = y - npLog(y) / y
+
+    z = 1.0 / y
+    if y < 8.0:
+        # Approximation for interval z = sqrt(-2 log x ) between 2 and 8
+        #  i.e., x between exp(-2) = .135 and exp(-32) = 1.27e-14.
+        y1 = z * hpoly(p1, z) / hpoly(q1, z)
+    else:
+        # Approximation for interval z = sqrt(-2 log x ) between 8 and 64
+        # i.e., x between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
+        y1 = z * hpoly(p2, z) / hpoly(q2, z)
+
+    y = y0 - y1
+
+    if negate: y = -y
+
+    return y
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 13c01cf..0499efc 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,7 @@ setup(
         "pandas_ta.volatility",
         "pandas_ta.volume"
     ],
-    version=".".join(("0", "3", "27b")),
+    version=".".join(("0", "3", "28b")),
     description=long_description,
     long_description=long_description,
     author="Kevin Johnson",
diff --git a/tests/test_utils.py b/tests/test_utils.py
index ce76d24..804c7dc 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -223,6 +223,18 @@ class TestUtilities(TestCase):
         self.assertEqual(self.utils.hpoly([1, 0, 1], 1), 2)
         self.assertEqual(self.utils.hpoly([1, 1, 1], 1), 3)
 
+    def test_inv_norm(self):
+        np.testing.assert_equal(self.utils.inv_norm(-0.01), np.nan)
+        self.assertEqual(self.utils.inv_norm(0), -np.infty)
+        self.assertEqual(self.utils.inv_norm(1 - 0.96), -1.7506860712521692)
+        self.assertAlmostEqual(self.utils.inv_norm(1 - 0.8646), -1.101222112591979)
+        self.assertEqual(self.utils.inv_norm(0.5), 0)
+        self.assertAlmostEqual(self.utils.inv_norm(0.8646), 1.101222112591979)
+        self.assertEqual(self.utils.inv_norm(0.96), 1.7506860712521692)
+        self.assertEqual(self.utils.inv_norm(1), np.infty)
+        np.testing.assert_equal(self.utils.inv_norm(1.01), np.nan)
+
+
     def test_linear_regression(self):
         x = Series([1, 2, 3, 4, 5])
         y = Series([1.8, 2.1, 2.7, 3.2, 4])