3

scipy.integrate.quad を使用して、非常に広い範囲 (0..10,000) で関数を統合しようとしています。関数はその範囲のほとんどでゼロですが、非常に狭い範囲 (1,602..1,618 など) でスパイクがあります。

統合すると、出力が正であると予想されますが、どういうわけかクワッドの推測アルゴリズムが混乱してゼロを出力していると思います。私が知りたいのは、これを克服する方法はありますか (たとえば、別のアルゴリズムや他のパラメーターを使用するなど)? 通常、どこでスパイクが発生するかはわかりません。そのため、積分範囲を分割して部分を合計することはできません (誰かがそれを行う方法について良いアイデアを持っていない限り)。

ありがとう！

出力例:

>>>scipy.integrate.quad(weighted_ftag_2, 0, 10000)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 0, 1602)
(0.0, 0.0)
>>>scipy.integrate.quad(weighted_ftag_2, 1602, 1618)
(3.2710994652983256, 3.6297354011338712e-014)
>>>scipy.integrate.quad(weighted_ftag_2, 1618, 10000)
(0.0, 0.0)
4

2 に答える 2

3

メソッドなど、他の統合方法を試してみることをお勧めしますintegrate.romberg()。

あるいは、関数が大きい点の位置を で取得し、weighted_ftag_2(x_samples).argmax()いくつかのヒューリスティックを使用して、関数の最大値付近で積分区間をカットすることができます (位置はです。サンプリングされた横座標のリスト ( )x_samples[….argmax()]をテイラーする必要があります)。x_samplesあなたの問題に対して：関数が最大になる領域にあるポイントを常に含む必要があります。

より一般的には、積分する関数に関する特定の情報は、その積分の適切な値を得るのに役立ちます。関数に適した方法 (Scipy が提供する多くの方法の 1 つ) と、統合間隔の合理的な分割 (たとえば、上記で提案された行に沿って) を組み合わせます。

于 2010-07-06T13:26:34.480 に答える
1

各整数範囲 [x, x+1) で関数 f() を評価し、たとえばromb()、EOL が示唆するように > 0 の場合は合計します。

from __future__ import division
import numpy as np
from scipy.integrate import romb

def romb_non0( f, a=0, b=10000, nromb=2**6+1, verbose=1 ):
    """ sum romb() over the [x, x+1) where f != 0 """
    sum_romb = 0
    for x in xrange( a, b ):
        y = f( np.arange( x, x+1, 1./nromb ))
        if y.any():
            r = romb( y, 1./nromb )
            sum_romb += r
            if verbose:
                print "info romb_non0: %d %.3g" % (x, r)  # , y
    return sum_romb

#...........................................................................
if __name__ == "__main__":
    np.set_printoptions( 2, threshold=100, suppress=True )  # .2f

    def f(x):
        return x if (10 <= x).all() and (x <= 12).all() \
            else np.zeros_like(x)

    romb_non0( f, verbose=1 )
于 2010-07-07T10:35:46.713 に答える