如何使用scipy.stats expect函数？

Question

我正期待着

scipy.stats.norm.expect(loc = 55, scale = 1)

返回分布55的平均值，而不是返回3.9096876333292135e-108。我用错了吗？有人能给我解释一下这个函数是如何工作的吗？谢谢

Answer 1

我想这应该是个bug。您已经正确地理解了手册，但事实证明，如果"loc"大于35.6，那么情况就会变得异常。否则，它会像预期的那样工作。

Answer 2

问题是：

import scipy.stats
#expect(func, loc=0, scale=1, lb=None, ub=None, conditional=False, **kwds)

for i in range(26, 55):
    print scipy.stats.norm.expect(loc=i,scale=1),

输出：

26.0
27.0
28.0
29.0
30.0
31.0
32.0
33.0
34.0
35.0
3.36373206533e-10
1.24696133938e-13
1.7005527212e-17
8.53164410131e-22
1.57463854604e-26
1.06913991702e-31
2.670507049e-37
2.45390887262e-43
8.29523418554e-50
1.03158350625e-56
4.71938759035e-64
7.94277793838e-72
4.91773037097e-80
1.12011477121e-88
9.38568238273e-98
2.89317497775e-107
3.28087032536e-117
1.4932849185e-127
3.63870849201e-118

对于较大的scale (标准差)，行为的“截止”会进一步向外移动，正如您所预期的那样。

problem can be 'fixed' by setting the lower and upper bounds explicitly，如下所示：

import numpy as np

for i in np.arange(5,100,5):
    print i,scipy.stats.norm.expect(loc=55,lb=-i,ub=i,scale=1) 

输出：

5 0.0
10 0.0
15 0.0
20 2.01210143973e-267
25 1.05364770562e-196
30 7.87517644756e-137
35 8.61623210583e-88
40 1.40277331283e-49
45 3.46495136419e-22
50 1.42791169386e-05
55 27.1010577196
60 54.9999827474
65 55.0
70 55.0
75 55.0
80 55.0
85 55.0
90 55.0
95 55.0

但这里面肯定有个细微的bug。如果您查看source (参见Numpy init类rv_continous的方法)，您会发现默认限制是作为Numpy的'inf‘导入的。如果以+/- np.inf的形式显式运行限制，则会得到与所述操作相同的行为：

for i in np.arange(5,60,5):
    print i,scipy.stats.norm.expect(loc=i,lb=-np.inf,ub=np.inf,scale=1) 

输出：

5 5.0
10 10.0
15 15.0
20 20.0
25 25.0000000007
30 30.0
35 35.0
40 1.57463854604e-26
45 1.03158350625e-56
50 9.38568238273e-98
55 3.90968763333e-108

还要注意，在source中的expect定义中，集成警告是静默的：

    # Silence floating point warnings from integration.
    olderr = np.seterr(all='ignore')
    vals = integrate.quad(fun, lb, ub, **kwds)[0] / invfac

问题的根源很可能在于当给定为+/- np.inf时，integrate.quad如何处理限制。在Wikipedia中给出了实际的source for the integration is in Fortran，但粗略地描述了如何对无限区间进行数值积分(映射到有限范围，就像黎曼球的情况一样)。