scipy quad 仅使用 1 个细分并给出错误的结果
scipy quad uses only 1 subdivision and gives wrong result
我想使用 quad 来获得高斯分布的均值。
我的第一次尝试和第二次尝试得到了不同的结果。 quad 的第二次尝试仅使用 1 个细分。
mu =1
sigma =2
import scipy as sp
import scipy.integrate as si
import scipy.stats as ss
f = lambda x: x * ss.norm(loc=mu, scale=sigma).pdf(x)
a = si.quad(f, -999., 1001., full_output=True)
print a[0]
#print sum(a[2]["rlist"][:a[2]["last"]])
print a[2]["last"]
b = si.quad(f, -1001., 1001., full_output=True)
print b[0]
#print sum(b[2]["rlist"][:b[2]["last"]])
print b[2]["last"]
print sorted(a[2]["alist"][:a[2]["last"]])
print sorted(b[2]["alist"][:b[2]["last"]])
这是输出:
1.0
16
0.0
1
[-999.0, -499.0, -249.0, -124.0, -61.5, -30.25, -14.625, -6.8125, 1.0, 8.8125, 16.625, 32.25, 63.5, 126.0, 251.0, 501.0]
[-1001.0]
我有没有搞错?
因为积分的极限远远超出了高斯分布的尾部,所以您误以为 quad
认为函数是相同的 0:
In [104]: f(-1000)
Out[104]: -0.0
In [105]: f(-500)
Out[105]: -0.0
In [106]: f(-80)
Out[106]: -0.0
In [107]: f(-50)
Out[107]: -6.2929842629835128e-141
您可以通过多种方式解决此问题,其中之一是将参数 points=[mu]
添加到对 quad
的调用中:
In [110]: b = si.quad(f, -1001., 1001., full_output=True, points=[mu])
b
In [111]: b[0]
Out[111]: 1.0000000000000002
我想使用 quad 来获得高斯分布的均值。 我的第一次尝试和第二次尝试得到了不同的结果。 quad 的第二次尝试仅使用 1 个细分。
mu =1
sigma =2
import scipy as sp
import scipy.integrate as si
import scipy.stats as ss
f = lambda x: x * ss.norm(loc=mu, scale=sigma).pdf(x)
a = si.quad(f, -999., 1001., full_output=True)
print a[0]
#print sum(a[2]["rlist"][:a[2]["last"]])
print a[2]["last"]
b = si.quad(f, -1001., 1001., full_output=True)
print b[0]
#print sum(b[2]["rlist"][:b[2]["last"]])
print b[2]["last"]
print sorted(a[2]["alist"][:a[2]["last"]])
print sorted(b[2]["alist"][:b[2]["last"]])
这是输出:
1.0
16
0.0
1
[-999.0, -499.0, -249.0, -124.0, -61.5, -30.25, -14.625, -6.8125, 1.0, 8.8125, 16.625, 32.25, 63.5, 126.0, 251.0, 501.0]
[-1001.0]
我有没有搞错?
因为积分的极限远远超出了高斯分布的尾部,所以您误以为 quad
认为函数是相同的 0:
In [104]: f(-1000)
Out[104]: -0.0
In [105]: f(-500)
Out[105]: -0.0
In [106]: f(-80)
Out[106]: -0.0
In [107]: f(-50)
Out[107]: -6.2929842629835128e-141
您可以通过多种方式解决此问题,其中之一是将参数 points=[mu]
添加到对 quad
的调用中:
In [110]: b = si.quad(f, -1001., 1001., full_output=True, points=[mu])
b
In [111]: b[0]
Out[111]: 1.0000000000000002