计算正态分布标准差的标准C或Python库

Question

假设我们有一个正态分布 n(x)，它的平均值是0，并且在区间从-a到a的总面积是P。

那么，计算这个分布的标准差最简单的方法是什么呢？也许在Python或C语言中有一些现成的库可以用来完成这个任务？

Answer 1

SciPy 有一个叫做 stats 的子包。

Answer 2

如果X是一个正态分布，平均值为0，标准差为sigma，那么必须满足以下条件：

P = Prob[ -a <= X <= a ] = Prob[ -a/sigma <= N <= a/sigma ]
  = 2 Prob[ 0 <= N <= a/sigma ]
  = 2 ( Prob[ N <= a/sigma ] - 1/2 )

这里的N也是一个正态分布，平均值为0，标准差为1。因此：

P/2 + 1/2 = Prob[ N <= a/sigma ] = Phi(a/sigma)

其中Phi是一个正态分布的累积分布函数（cdf），它的平均值是0，标准差是1。现在我们需要的是正态分布的反累积分布函数（或者叫“百分位点函数”），在Python中可以用scipy.stats.norm.ppf()来实现。下面是一个示例代码：

from scipy.stats import norm
P = 0.3456
a = 3.0

a_sigma = float(norm.ppf(P/2 + 0.5))   # a/sigma
sigma = a/a_sigma   # Here is the standard deviation

举个例子，我们知道N(0,1)这个变量落在区间[-1, 1]内的概率大约是0.682（在这个图中是深蓝色的区域）。如果你设定P = 0.682，a = 1.0，那么你会得到sigma大约等于1.0，这确实就是标准差。

Answer 3

一个均值为零的高斯分布的标准差，满足条件 Pr(-a < X < a) = P，可以用下面的公式表示：

a/(sqrt(2)*inverseErf(P))

这个公式就是你需要的，其中 inverseErf 是误差函数的反函数，通常我们称之为 erf。

如果你使用 C 语言，Gnu 科学库（GSL）是一个不错的资源。不过，它只提供了 erf 函数，没有 inverseErf，所以你需要自己去计算反函数（简单的二分查找就可以解决这个问题）。另外，这里有一个很好的方法来近似计算 erf 和 inverseErf：

http://homepages.physik.uni-muenchen.de/~Winitzki/erf-approx.pdf

如果你使用 Python，inverseErf 可以在 SciPy 库中找到，叫做 erfinv，所以你可以用下面的代码来计算标准差：

a/(math.sqrt(2)*erfinv(P))

附言：Stackoverflow 的链接显示有点问题，我无法在上面链接到 GSL： http://www.gnu.org/software/gsl。而且当我把上面的 PDF 链接做成正确的链接时，它也显示错误。