Scipy积分出错

-1

我试图集成一个函数。该功能是保证非负：Scipy积分出错

def function(x): 
    something = ... 
    something_else = ... 
    return exp(something)/sqrt(something_else)

现在我正在整合它：

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b)

我得到的结果不是我所期望的，所以要检查我这样做：

for x in range(0, 10000): 
    if integrand(0,x+1) < integrand(0,x): 
     raise ValueError("Weird!")

果然，我越来越'奇怪'的例外。怎么可能？

来源

2014-11-09 zmbq

你能举一个实际的例子来证明这个问题吗？ – BrenBarn 2014-11-09 05:38:04

这很难，因为实际的例子涉及10,000点C14年龄校准曲线...... – zmbq 2014-11-09 05:39:29

那么，如果没有一个小的测试用例，调试问题也很难。尝试将问题简化为一个小测试用例。 – BrenBarn 2014-11-09 05:52:45

没有一个工作的例子，很难调试你的代码，但这里有一些基于你所显示的建议。

quad返回一个元组，其中包含两个值：积分估计值和结果绝对误差估计值。您似乎没有使用代码中的第一个值。尝试调用quad后更改integrand到

def integrand(a, b): 
    return scipy.integrate.quad(function, a, b)[0]

注意[0]。

接下来，当您发现“奇怪”的情况时，您应该打印integrand(0, x)和integrand(0, x+1)的值。如果它们大致相同，那么正常的数值不精确可能会破坏您理论上预期的顺序。看看quad返回的第二个值。如果该值大于integrand(0, x)和integrand(0, x+1)之间的差值，则不能真正期望积分的估计值单调递增，因为预期的增加量小于计算中的误差。如果这是问题，您可以尝试将epsabs和/或epsrel设置为比1.49e-8的默认值更小的值，但这只会让您感到满意。最终你会遇到浮点计算的限制。

来源

2014-11-09 15:00:49

Scipy积分出错

回答

相关问题