Scipy free implementation of Normal distribution inverse CDF

normal_inverse_cdf.py

def inverse_normal_cdf(p, mean, std):
    """
    This is the inverse to a normal distribution's CDF.

    While much slower this means you do not need Scipy as a project requirement.
    :param p: list of p = (0, 1)
    :param mean:
    :param std:
    :return:
    """
    p = np.array(p)
    z_scores = np.zeros(p.shape)

    def i_inverse_cdf(p_i):
        if type(p_i) is float:
            assert 0 < p_i < 1, 'P = (0, 1)'

        def polevl(x, coef):
            accum = 0
            for c in coef:
                accum = x * accum + c
            return accum

        P0 = [
            -5.99633501014107895267E1,
            9.80010754185999661536E1,
            -5.66762857469070293439E1,
            1.39312609387279679503E1,
            -1.23916583867381258016E0,
        ]
        Q0 = [
            1,
            1.95448858338141759834E0,
            4.67627912898881538453E0,
            8.63602421390890590575E1,
            -2.25462687854119370527E2,
            2.00260212380060660359E2,
            -8.20372256168333339912E1,
            1.59056225126211695515E1,
            -1.18331621121330003142E0,
        ]
        P1 = [
            4.05544892305962419923E0,
            3.15251094599893866154E1,
            5.71628192246421288162E1,
            4.40805073893200834700E1,
            1.46849561928858024014E1,
            2.18663306850790267539E0,
            -1.40256079171354495875E-1,
            -3.50424626827848203418E-2,
            -8.57456785154685413611E-4,
        ]
        Q1 = [
            1,
            1.57799883256466749731E1,
            4.53907635128879210584E1,
            4.13172038254672030440E1,
            1.50425385692907503408E1,
            2.50464946208309415979E0,
            -1.42182922854787788574E-1,
            -3.80806407691578277194E-2,
            -9.33259480895457427372E-4,
        ]
        P2 = [
            3.23774891776946035970E0,
            6.91522889068984211695E0,
            3.93881025292474443415E0,
            1.33303460815807542389E0,
            2.01485389549179081538E-1,
            1.23716634817820021358E-2,
            3.01581553508235416007E-4,
            2.65806974686737550832E-6,
            6.23974539184983293730E-9,
        ]
        Q2 = [
            1,
            6.02427039364742014255E0,
            3.67983563856160859403E0,
            1.37702099489081330271E0,
            2.16236993594496635890E-1,
            1.34204006088543189037E-2,
            3.28014464682127739104E-4,
            2.89247864745380683936E-6,
            6.79019408009981274425E-9,
        ]

        negate = True
        y = p_i
        if y > 1.0 - 0.13533528323661269189:
            y = 1.0 - y
            negate = False

        if y > 0.13533528323661269189:
            y -= 0.5
            y2 = y * y
            x = y + y * (y2 * polevl(y2, P0) / polevl(y2, Q0))
            x = x * np.sqrt(2 * np.pi)
            return x

        x = np.sqrt(-2.0 * np.log(y))
        x0 = x - np.log(x) / x

        z = 1.0 / x
        if x < 8.0:
            x1 = z * polevl(z, P1) / polevl(z, Q1)
        else:
            x1 = z * polevl(z, P2) / polevl(z, Q2)
        x = x0 - x1
        if negate:
            x = -x
        return x

    for z_i in np.nditer(z_scores, op_flags=['readwrite']):
        z_i[...] = i_inverse_cdf(z_i)

    # Transform by with by mean and std
    return z_scores * std + mean