scipytegral.quad返回不正确的值 [英] scipy integrate.quad return an incorrect value

问题描述

我使用scipy integration.quad计算正态分布的cdf:

i use scipy integrate.quad to calc cdf of normal distribution:

def nor(delta, mu, x):
    return 1 / (math.sqrt(2 * math.pi) * delta) * np.exp(-np.square(x - mu) / (2 * np.square(delta)))


delta = 0.1
mu = 0
t = np.arange(4.0, 10.0, 1)
nor_int = lambda t: integrate.quad(lambda x: nor(delta, mu, x), -np.inf, t)
nor_int_vec = np.vectorize(nor_int)

s = nor_int_vec(t)
for i in zip(s[0],s[1]): 
    print i

其打印如下:

(1.0000000000000002, 1.2506543424265854e-08)
(1.9563704110140217e-11, 3.5403445591955275e-11)
(1.0000000000001916, 1.2616577562700088e-08)
(1.0842532749783998e-34, 1.9621183122960244e-34)
(4.234531567162006e-09, 7.753407284370446e-09)
(1.0000000000001334, 1.757986959115912e-10)

对于某些x，它返回的值近似为零，应该返回1. 有人可以告诉我怎么了吗?

for some x, it return a value approximate to zero, it should be return 1. can somebody tell me what is wrong?

推荐答案

与

Same reason as in why does quad return both zeros when integrating a simple Gaussian pdf at a very small variance? but seeing as I can't mark it as a duplicate, here goes:

您正在集成一个很大的(实际上是无限的)时间间隔内具有紧密定位(按比例增量)的函数.积分例程可以简单地忽略该函数实质上不同于0的区间部分，而是将其判断为0.需要一些指导.可以使用参数points来达到此目的(请参阅链接的问题)，但是由于在无限间隔内的quad不支持该参数，因此必须手动分割间隔，如下所示:

You are integrating a function with tight localization (at scale delta) over a very large (in fact infinite) interval. The integration routine can simply miss the part of the interval where the function is substantially different from 0, judging it to be 0 instead. Some guidance is required. The parameter points can be used to this effect (see the linked question) but since quad over an infinite interval does not support it, the interval has to be manually split, like so:

for t in range(4, 10):
    int1 = integrate.quad(lambda x: nor(delta, mu, x), -np.inf, mu - 10*delta)[0] 
    int2 = integrate.quad(lambda x: nor(delta, mu, x), mu - 10*delta, t)[0] 
    print(int1 + int2)

每次打印1或将近1.我选择mu-10*delta作为切入点，无论mu和delta是多少，大多数功能都位于它的右侧.

This prints 1 or nearly 1 every time. I picked mu-10*delta as a point to split on, figuring most of the function lies to the right of it, no matter what mu and delta are.

注意:

1. 使用np.sqrt等；通常没有理由在NumPy代码中放置math函数. NumPy版本可用并且已矢量化.
2. 将np.vectorize应用于quad并没有做任何事情，除了使代码更长且更难阅读.使用普通的Python循环或列表理解.请参见具有集成功能的NumPy矢量化

1. Use np.sqrt etc; there is usually no reason for put math functions in NumPy code. The NumPy versions are available and are vectorized.
2. Applying np.vectorize to quad is not doing anything besides making the code longer and slightly harder to read. Use a normal Python loop or list comprehension. See NumPy vectorization with integration

这篇关于scipytegral.quad返回不正确的值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！