时间戳

首页 > 解决方案 > scipy中的理论正态分布函数

python - scipy中的理论正态分布函数

问题描述

我需要为给定的垃圾箱边缘绘制正态累积分布：

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() 结果看起来更糟。我究竟做错了什么？

标签： pythonnumpyscipystatisticsprobability

解决方案

问题

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

在此处输入图像描述

你的结果与你的结果不一致的原因是你scipy的scipy积分比你做得更好。通过对您bin_edges有效定义的直方图“条形”区域求和，您可以有效地整合普通 PDF。在您的垃圾箱数量非常非常多（可能至少有数千个）之前，这不会产生相当准确的结果。您的规范化方法也已关闭，因为实际上您需要除以 PDF 的积分 from -infto inf，而不是 from 1.02to 28.7。

另一方面，Numpy 只是计算积分的封闭形式解的高精度数值近似。它使用的函数称为scipy.special.ndtr. 这是它在 Scipy 代码中的实现。

解决方案

-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.

推荐阅读