python - Scipy 的四元组在 "pretty smooth"函数上失败

Question

数据生成函数相当困惑，所以我会尽力说得清楚。如果做得不够好，请发表评论。

我有一个包含多个变量的函数，并尝试对它进行积分，同时保持另一个变量不变。然而，对于次要变量的不同值，积分过程会产生完全不同的结果 - 尽管函数变化不大!

这是我的示例代码(不幸的是，不可重现):

fig, ax = plt.subplots()
tPrimeGrid = [0.16, 0.18, 0.22]
from scipy import integrate
for ttPrime in tPrimeGrid:
    # compute grid
    lowerBar = continuousTime.getPLowerBarAnalytical(ttPrime, Param.S, Param)
    upperBar = Param.r

    # integrate function
    h = integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, ttPrime, Param), full_output=True)
    print('tPrime {} value {} erorr {}'.format(ttPrime, h[0], h[1]))
    # plot function to be evaluated
    pGrid = np.linspace(lowerBar, upperBar, 100)
    plt.plot(pGrid, [myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid], label='t {} (h: {})'.format(ttPrime, h[0]))
plt.legend()

quad 的输出为

tPrime 0.16 value 6.63310536371 erorr 0.000345564616621
tPrime 0.18 value 5.93645658492 erorr 0.00045956820487
tPrime 0.22 value 0.359208009237 erorr 3.98801002485e-15

错误都在我的容忍范围内，但值跳跃了很多。这是图表:

所以这些函数具有非常相似的形状。用肉眼看，积分的差异确实没有意义。然后我“手动”计算它们:

(np.array([myFunc(pp, 0.6, ttPrime, Param) for pp in pGrid])).sum()*(pGrid[1]-pGrid[0])
0.35538925924926557

这意味着较小的值实际上是正确的。因此，我还将在此处添加不良输出之一的 quad() 的完整输出:

integrate.quad(myFunc, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True)
Out[19]: 
(6.634157704675579,
 0.004721148834250418,
 {'alist': array([ 1.99994895,  1.78826738,  1.86060326,  1.79090489,  1.93030163,
          1.72120652,  1.96515082,  1.75605571,  1.98257541,  1.7734803 ,
          1.9912877 ,  1.7821926 ,  1.99564385,  1.78872682,  1.99782193,
          1.78654874,  1.99940018,  1.78763778,  1.99945548,  1.99891096,
          1.78845456,  1.99972774,  1.99918322,  1.78831843,  1.99988939,
          1.99931935,  1.7881823 ,  1.99999149,  1.99942145,  1.7882844 ,
          1.9998979 ,  1.9999447 ,  1.99940443,  1.78825036,  1.99996597,
          1.99986387,  1.99991492,  1.99995746,  1.99938742,  1.78827589,
          1.99993194,  1.99990641,  1.99998298,  1.99988089,  1.99995321,
          1.99939592,  1.78827164,  1.99994044,  1.99995108,  1.99940231]),
  'blist': array([ 1.99995108,  1.78827164,  1.93030163,  1.86060326,  1.96515082,
          1.75605571,  1.98257541,  1.7734803 ,  1.9912877 ,  1.7821926 ,
          1.99564385,  1.78654874,  1.99782193,  1.79090489,  1.99891096,
          1.78763778,  1.99940231,  1.7881823 ,  1.99972774,  1.99918322,
          1.78872682,  1.99986387,  1.99931935,  1.78845456,  1.9998979 ,
          1.99938742,  1.78825036,  2.        ,  1.99945548,  1.78831843,
          1.99990641,  1.99994895,  1.99942145,  1.78826738,  1.99998298,
          1.99988089,  1.99993194,  1.99996597,  1.99939592,  1.7882844 ,
          1.99994044,  1.99991492,  1.99999149,  1.99988939,  1.99995746,
          1.99940018,  1.78827589,  1.9999447 ,  1.99995321,  1.99940443]),
  'elist': array([  4.61134496e-03,   4.14135579e-06,   1.27804393e-15,
           1.55808958e-15,   5.54429458e-16,   0.00000000e+00,
           2.90617202e-16,   0.00000000e+00,   2.27720143e-15,
           0.00000000e+00,   3.91514838e-15,   0.00000000e+00,
           2.15987477e-14,   5.40749179e-17,   1.19682144e-12,
           0.00000000e+00,   1.03521880e-04,   0.00000000e+00,
           2.48275100e-08,   6.85735811e-13,   6.78454062e-18,
           8.65204618e-09,   1.64268148e-13,   3.39437556e-18,
           4.65487056e-08,   1.12644353e-13,   0.00000000e+00,
           6.19443638e-08,   4.08066769e-09,   8.48813317e-19,
           5.99147851e-07,   1.21478039e-06,   9.39741081e-10,
           0.00000000e+00,   6.32087778e-14,   2.19912539e-10,
           6.97448944e-09,   1.85114515e-13,   1.28558440e-13,
           2.12217047e-19,   8.89451362e-08,   8.45167083e-09,
           1.25621961e-14,   2.59393721e-08,   2.45738052e-15,
           1.19106919e-12,   1.06110581e-19,   4.91812027e-08,
           2.72831861e-15,   3.06819516e-12]),
  'iord': array([ 1, 17, 50, 49, 48, 47, 46, 45, 44, 43, 42, 29, 40, 39, 38, 23, 35,
         11, 34,  3,  7, 21, 30, 18, 12,  8,  6,  0,  0,  0,  0,  0,  0,  0,
          0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0], dtype=int32),
  'last': 50,
  'neval': 2079,
  'rlist': array([  1.04205921e-01,   6.52426361e-06,   1.15115963e-01,
           1.40340233e-01,   4.99385660e-02,   0.00000000e+00,
           2.33230318e-02,   0.00000000e+00,   1.12779305e-02,
           0.00000000e+00,   5.54633015e-03,   0.00000000e+00,
           2.75040018e-03,   4.87063561e-03,   1.36955727e-03,
           0.00000000e+00,   8.85474806e-03,   0.00000000e+00,
           1.73731989e+00,   3.41803845e-04,   6.11097092e-04,
           1.73678061e+00,   1.70814248e-04,   3.05738170e-04,
           2.00530044e-01,   8.53852175e-05,   0.00000000e+00,
           5.29681716e-06,   1.51991445e-01,   7.64543067e-05,
           2.17985628e-01,   2.00512965e-01,   7.26773425e-02,
           0.00000000e+00,   1.06564581e-05,   3.34545761e-01,
           5.59012296e-01,   5.32838374e-06,   1.06721256e-05,
           1.91148122e-05,   3.34512038e-01,   2.38773007e-01,
           5.32807434e-06,   1.85664300e-01,   2.66423055e-06,
           5.33597722e-06,   9.55759147e-06,   1.85647516e-01,
           1.33212494e-06,   8.93844378e-03])},
 'The maximum number of subdivisions (50) has been achieved.\n  If increasing the limit yields no improvement it is advised to analyze \n  the integrand in order to determine the difficulties.  If the position of a \n  local difficulty can be determined (singularity, discontinuity) one will \n  probably gain from splitting up the interval and calling the integrator \n  on the subranges.  Perhaps a special-purpose integrator should be used.')

这是我的相关问题:

1. 形状如此相似，怎么可能有些参数起作用而有些参数不起作用？我如何“预测” future 这些失败？
2. 完整的输出(如下)告诉我已经实现了最大分割数。它并没有告诉我该错误可能被错误估计。这一定是暗示的吧？ 5.9-0.0004 与 0.3 并不接近。
3. 由于增加限制没有帮助(见下文)，潜在的替代方案是什么？如何集成这个功能？

我尝试增加限制，但这只给了我以下(而不是更好)的输出:

integrate.quad(expectedUtility, lowerBar, upperBar, args=(0.6, 0.16, Param), full_output=True, limit=500)
Out[24]: 
(6.633112814769514,
 4.74743687572826e-06,
[...]
 'The occurrence of roundoff error is detected, which prevents \n  the requested tolerance from being achieved.  The error may be \n  underestimated.')

Answer 1

当 full_output 为 True 时，

quad 对于返回值有一个有点时髦的约定。如果quad没有检测到任何问题，它会返回y、abserr、infodict。如果quad检测到某种形式的故障，它会返回y、abserr、infodict、message。 infodict 不包含简单的 status 字段来指示失败，因此您必须检查第四个返回值是否存在。如果存在，则它是描述问题的字符串。 (在您的代码中，您可以添加类似 if len(h) > 3: ...handle the failure... 的内容。)在不良情况的完整输出中，您可以看到 message 是包含以下内容的字符串:

The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties.  If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges.  Perhaps a special-purpose integrator should be used.

这意味着数值积分失败。如果不了解有关被积函数 myFunc 的更多信息，就很难说出它失败的原因。