如何在Python中计算正态分布百分点函数

Question

我该如何在不使用Scipy的情况下完成与scipy.stats.norm.ppf等效的操作。我已经在Python的Math模块中找到了内置的erf函数，但我似乎无法重现它的功能。

PS：我不能使用Scipy，因为Heroku不允许你安装它，而使用其他的buildpacks会超出300Mb的最大slug大小限制。

- SARose

Answer 1

函数ppf是y = (1+erf(x/sqrt(2))/2的反函数。因此，我们需要解决这个方程，找到x，给定y在0和1之间。下面的代码使用二分法来完成这个任务。我导入了SciPy函数来说明结果是相同的。

from math import erf, sqrt
from scipy.stats import norm         # only for comparison
y = 0.123  

z = 2*y-1
a = 0
while erf(a) > z or erf(a+1) < z:    # looking for initial bracket of size 1
    if erf(a) > z:
        a -= 1
    else:
        a += 1
b = a+1                              # found a bracket, proceed to refine it
while b-a > 1e-15:                   # 1e-15 ought to be enough precision 
   c = (a+b)/2.0                     # bisection method
   if erf(c) > z:
       b = c
   else:
       a = c

print sqrt(2)*(a+b)/2.0              # this is the answer 
print norm.ppf(y)                    # SciPy for comparison

你需要完成以下工作：

• 进行初步边界检查（y必须在0和1之间）
• 如果需要其他平均值/方差，则进行缩放和移位；代码适用于标准正态分布（平均值为0，方差为1）。

- user6655984

---

非常赞赏您使用二分法（并提醒我使用它），我只是想看看这有多快，然后我会将其标记为正确的。 - SARose

为了实际使用，可以使用割线法或相关的方法，如伊利诺伊的弗尔西规则变种，以实现超线性收敛。 - Lutz Lehmann

@LutzL 是的，这不是一个快速的解决方案，但是erf函数远离0点时非常平坦，这会让我担心割线法或牛顿法（牛顿法是一种选择，因为erf的导数是显式的）。Regula falsi可以工作。无论如何，通用的根查找器无法与SciPy自己实现的ppf相竞争，因此这个答案更像是“如果我们没有SciPy源代码会怎样”。 - user6655984

Answer 2

使用erf实现norm.ppf并不简单，因为norm.ppf与erf的反函数有关。相反，这里提供了一个纯Python实现代码，它与scipy中的代码完全相同。您会发现函数ndtri返回与norm.ppf完全相同的值：

import math

s2pi = 2.50662827463100050242E0

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,
]

def ndtri(y0):
    if y0 <= 0 or y0 >= 1:
        raise ValueError("ndtri(x) needs 0 < x < 1")
    negate = True
    y = y0
    if y > 1.0 - 0.13533528323661269189:
        y = 1.0 - y
        negate = False

    if y > 0.13533528323661269189:
        y = y - 0.5
        y2 = y * y
        x = y + y * (y2 * polevl(y2, P0) / polevl(y2, Q0))
        x = x * s2pi
        return x

    x = math.sqrt(-2.0 * math.log(y))
    x0 = x - math.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

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

- K. A. Buhr

---

这些数字是什么意思？ - SARose

它们是预先计算的值，用于在最小计算量的情况下获得最准确的值。您可以在此处找到SciPy使用的原始代码（https://raw.githubusercontent.com/scipy/scipy/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/special/cephes/ndtri.c）；其中包含有关其中一些值含义的注释，但您不应该考虑它们的含义。 - K. A. Buhr

同时，如果您要将标准分数返回到原始输出，则必须进行转换。公式为 z * std + mean。 - SARose