在 IronPython 中计算正态概率

John D. Cook

5.00/5 (2投票s)

2009年3月10日

BSD

4分钟阅读

29471

用于计算正态概率的代码，以及关于在 IronPython 中什么可行、什么不可行的讨论。

引言

统计学中最常见的任务之一是计算正态（高斯）随机变量取两个边界之间的值的概率。这用于计算区间估计、构建假设检验以及为许多问题进行粗略估计。

Python 的 SciPy 数值库包含用于计算此类概率的函数。但是，SciPy 目前不适用于 IronPython。本文将解释为什么某些 Python 库适用于 IronPython，而另一些则不适用。此外，还包括用于计算正态概率的独立代码，该代码适用于任何版本的 Python。

背景

概率和统计学中最常见的分布是“正态”或高斯分布。人们通常称这种分布为“钟形曲线”。这种分布随处可见：身高、智商分数、测量误差等。

正态随机变量有两个参数。均值 μ 描述分布的中心在哪里。标准差 σ 描述分布的离散程度。如果 X 是均值为 μ 且标准差为 σ 的正态随机变量，则 X 取 a 和 b 之间值的概率由下面的积分给出。

$\frac{1}{\sqrt{2\pi}\, \sigma}\int_a^b \exp\left(-\frac{(t-\mu)^2}{2\sigma^2}\right)\, dt$

这个积分不能用更简单的函数来计算。多年前，人们会在表格书中查找积分值。当然，现在人们使用软件而不是表格，这就是本文的主题。

起初，它看起来像是一个包含四个变量的函数：a、b、μ 和 σ，但我们可以将问题简化为仅计算一个变量的函数。首先，我们可以假设 μ = 0 且 σ = 1。这是因为如果 X 是均值为 μ 且标准差为 σ 的正态随机变量，则 Z = (X-μ)/σ 是均值为 0 且标准差为 1 的正态随机变量。因此，概率 P(a < X < b) 等于概率 P((a-μ)/σ < Z < (b-μ)/σ。此外，我们可以假设 b 为 -∞。这是因为 P(a < Z < b) = P(Z < b) - P(Z < a)。因此，我们将问题简化为需要评估以下函数。

$\Phi(x) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^x e^{-t^2/2}\, dt$

（Iron）Python 中的数值计算

Python 不包含用于计算函数 Φ(z) 的函数。但是，SciPy 库包含大量像 Φ (Phi) 这样的数学函数。（“SciPy”代表“科学 Python”。）在 Python 中寻找高级数学支持时，首选 SciPy。该代码经过良好测试且高效（尽管并非总是得到很好的文档记录）。

以下示例演示了如何使用 SciPy 从 Python 命令行计算 phi(0.7)。

>>> from scipy.stats.distributions import norm
>>> norm.cdf(0.7)
0.75803634777692697

但是，有时您可能希望避免创建额外的库依赖性。例如，您可能想通过电子邮件向某人发送一个快速 Python 脚本，并且不希望假设他们已安装 SciPy。

有时，即使人们想使用 SciPy 也无法使用。IronPython 就是这种情况。

IronPython 是 Microsoft 为 .NET 提供的 Python 实现。 IronPython 允许您编写 Python 语法并访问任何 .NET 代码。但是，它并不总是让从 .NET 调用 Python 变得容易。从这个意义上讲，IronPython 是一扇单向门：它使 Python 代码可以轻松调用 .NET，但不一定使 .NET 代码可以轻松调用 Python。困难在于 Python 库，而不是 Python 语言。一些 Python 库是用 C 实现的，而另一些是用 Python 实现的。一个典型的 Python 程序员不需要意识到这种区别，但是一个 IronPython 程序员需要。 IronPython 可以轻松调用用纯 Python 编写的代码，但是用 C 实现的 Python 代码会导致 .NET 的托管代码模型出现问题。 Ironclad 项目解决了这个问题。此外，来自 Microsoft 的人员表示，随着 .NET 和 DLR（动态语言运行时，是 .NET 对 IronPython 和 IronRuby 等动态语言的补充）的未来版本，这种情况应该会得到改善。

Phi(x) 的独立 Python 代码

这是函数 Φ(x) 的独立 Python 实现。

import math

def phi(x):
    # constants
    a1 =  0.254829592
    a2 = -0.284496736
    a3 =  1.421413741
    a4 = -1.453152027
    a5 =  1.061405429
    p  =  0.3275911

    # Save the sign of x
    sign = 1
    if x < 0:
        sign = -1
    x = abs(x)/math.sqrt(2.0)

    # A&S formula 7.1.26
    t = 1.0/(1.0 + p*x)
    y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*math.exp(-x*x)

    return 0.5*(1.0 + sign*y)

该算法基于 Abramowitz 和 Stegun 的《数学函数手册》（被亲切地称为“A&S”）中的公式 7.1.26。 A&S 中给出的算法用于计算误差函数 erf(x)。误差函数与函数 Phi 通过简单转换相关。（SciPy 中相应的函数是 scipy.special.erf。）

因为此代码除了标准的 math 模块之外没有任何依赖项，所以它应该可以在 Python 运行的任何地方运行。它已经在 Python 2.5、Python 3.0 和 IronPython 上进行了测试。该代码应精确到小数点后七位。

现在我们有了 phi(x) 的实现，我们可以围绕它编写一个简单的包装器来计算正态概率。

# Compute the probability that a normal(mu, sigma) random variable
# takes on a value between a and b.

def normal_probability(a, b, mu = 0, sigma = 1):
    return phi((b-mu)/sigma) - phi((a-mu)/sigma)