python - 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 有效定义的直方图的“条形”区域求和，您可以有效地集成普通 PDF。在您的 bin 计数高得多(可能至少有数千个)之前，这不会产生相当准确的结果。您的规范化方法也已关闭，因为您实际上需要除以从 -inf 到 inf 的 PDF 积分，而不是从 1.02到 28.7。

另一方面，Numpy 只是计算积分的封闭形式解的高精度数值近似。它使用的函数称为 scipy.special.ndtr .这是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]

ndtr 的 Scipy 版本是用 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 的结果高度匹配。