scipy 中的理论正态分布函数

Question

我需要绘制给定 bin 边缘的正态累积分布：

bin_edges = np.array([1.02,  4.98,  8.93, 12.89, 16.84, 20.79, 24.75, 28.7])
mean = 15.425
standard_deviation = 6.159900567379315

我首先做了：

cdf = ((1 / (np.sqrt(2 * np.pi) * standard_deviation)) *
   np.exp(-0.5 * (1 / standard_deviation * (bin_edges - mean))**2))
cdf = cdf.cumsum()
cdf /= cdf[-1]

我发现的另一种方式：

cdf = scipy.stats.norm.cdf(bin_edges, loc=mean, scale=standard_deviation)

这两种方法的输出应该是相等的，但事实并非如此：

First: [0.0168047  0.07815162 0.22646339 0.46391741 0.71568769 0.89247475 
0.97468339 1.]
Second: [0.0096921  0.04493372 0.14591031 0.34010566 0.59087116 0.80832701
0.93495018 0.98444529]

对我来说 scipy cdf() 结果更糟。我做错了什么？

Answer 1

问题

您正在尝试通过计算每个 bin 边缘的以下积分值来计算每个 bin 边缘的 CDF：

您的结果与scipy的结果不一致的原因是scipy比您做得更好。通过对 bin_edges 有效定义的直方图的 "bars" 区域求和，您可以有效地整合普通 PDF。在您的 bin 计数高得多（可能至少有数千个）之前，这不会产生相当准确的结果。您的标准化方法也已关闭，因为您确实需要除以 PDF 从 -inf 到 inf 的积分，而不是从 1.02 到 28.7.

另一方面，Numpy 只是计算积分的封闭形式解的高精度数值近似。它使用的函数叫做scipy.special.ndtr. Here's it's implementation in the Scipy code。

解决方案

您可以从 -inf 到 x 进行实际的数值积分，而不是通过对条形面积求和来积分，以获得精度接近 scipy.stats.norm.cdf 的结果。下面是如何执行此操作的代码：

import scipy.integrate as snt

def pdf(x, mean, std):
    return ((1/((2*np.pi)**.5 * std)) * np.exp(-.5*((x - mean)/std)**2))

cdf = [snt.quad(pdf, -np.inf, x, args=(mean, std))[0] for x in bin_edges]

Scipy 的 ndtr 版本是用 C 语言编写的，但这里有一个接近 Python 的近似值以供比较：

import scipy.special as sps

def ndtr(x, mean, std):
    return .5 + .5*sps.erf((x - mean)/(std * 2**.5))

正在测试

import scipy.special as sps
import scipy.stats as sts
import scipy.integrate as snt

bin_edges = np.array([1.02,  4.98,  8.93, 12.89, 16.84, 20.79, 24.75, 28.7])
mean = 15.425
std = 6.159900567379315

with np.printoptions(linewidth=9999):
    print(np.array([snt.quad(pdf, -np.inf, x, args=(mean, std))[0] for x in bin_edges]))
    print(ndtr(bin_edges, mean, std))
    print(sts.norm.cdf(bin_edges, loc=mean, scale=std))

输出：

[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]
[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]
[0.00968036 0.04497664 0.14584988 0.34034101 0.59084202 0.80811081 0.93496465 0.98442171]

因此，当您准确积分时，您使用的方法的结果与 scipy.stats.norm.cdf 的结果具有高精度匹配。

scipy 中的理论正态分布函数

Theoretical normal distribution function in scipy

python

statistics

numpy

probability

scipy

问题

解决方案

正在测试