3

ネストされたscipy.integrate.quad呼び出しを使用して、2次元被積分関数を統合しています。被積分関数はnumpy関数で構成されているため、入力をループして入力ごとに1回呼び出すよりも、入力の配列を渡す方がはるかに効率的です。numpyの配列により、約2桁高速になります。

ただし、被積分関数を1つの次元のみに統合したいが、他の次元に入力の配列がある場合は、問題が発生します。「scipy」クアッドパックパッケージでは、それが何でもできないようです。 numpyは配列された入力を処理します。他の誰かがこれを見たことがありますか-そして/またはそれを修正する方法を見つけました-または私はそれを誤解していますか？クワッドから得られるエラーは次のとおりです。 

Traceback (most recent call last):
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 312, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
quadpack.error: Supplied function does not return a valid float.

私がやろうとしていることの漫画版を下に置きました-私が実際にやっていることはより複雑な被積分関数を持っていますが、これはジャイストです。

肉は上にあります-下は私のポイントを示すためにベンチマークを行っています。 

import numpy as np
import time

from scipy.integrate import quad


def Integrand(x, y):
    '''
        Integrand
    '''
    return np.sin(x)*np.sin( y )

def Integrate_x(y):
    '''
        Integrate over x given (y)
    '''
    return quad(Integrand, 0, np.pi/2, args=(y))[0]



def fnIntegrate_x(ystart, yend, nsteps, ArrayInput = False):
    '''

    '''

    yarray = np.arange(ystart,yend, (yend - ystart)/float(nsteps))
    I = np.zeros(nsteps)
    if ArrayInput :
        I = Integrate_x(yarray)
    else :
        for i,y in enumerate(yarray) :

            I[i] = Integrate_x(y)

    return y, I




NCALLS_SET = 1000
NSETS = 10

SETS_t = np.zeros(NSETS)

for i in np.arange(NSETS) :

    XInputs = np.random.rand(NCALLS_SET, 2)

    t0 = time.time()
    for x in XInputs :
        Integrand(x[0], x[1])
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET

print "Benchmarking Integrand - Single Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
Benchmarking Integrand - Single Values:
NCALLS_SET:  1000
NSETS:  10
TimePerCall(s):  1.23999834061e-05 4.06987868647e-06
'''








NCALLS_SET = 1000
NSETS = 10

SETS_t = np.zeros(NSETS)

for i in np.arange(NSETS) :

    XInputs = np.random.rand(NCALLS_SET, 2)

    t0 = time.time()
    Integrand(XInputs[:,0], XInputs[:,1])
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET

print "Benchmarking Integrand - Array Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
Benchmarking Integrand - Array Values:
NCALLS_SET:  1000
NSETS:  10
TimePerCall(s):  2.00009346008e-07 1.26497018465e-07
'''












NCALLS_SET = 1000
NSETS = 100

SETS_t = np.zeros(NSETS)

for i in np.arange(NSETS) :


    t0 = time.time()
    fnIntegrate_x(0, 1, NCALLS_SET, False)
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET

print "Benchmarking fnIntegrate_x - Single Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        
print "TotalTime: ",np.sum(SETS_t) * NCALLS_SET
'''
NCALLS_SET:  1000
NSETS:  100
TimePerCall(s):  0.000165750000477 8.61204306241e-07
TotalTime:  16.5750000477
'''








NCALLS_SET = 1000
NSETS = 100

SETS_t = np.zeros(NSETS)

for i in np.arange(NSETS) :


    t0 = time.time()
    fnIntegrate_x(0, 1, NCALLS_SET, True)
    t1 = time.time()
    SETS_t[i] = (t1 - t0)/NCALLS_SET

print "Benchmarking fnIntegrate_x - Array Values:"
print "NCALLS_SET: ", NCALLS_SET
print "NSETS: ", NSETS    
print "TimePerCall(s): ", np.mean(SETS_t) , np.std(SETS_t)/ np.sqrt(SETS_t.size)        

'''
****  Doesn't  work!!!! *****
Traceback (most recent call last):
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 159, in <module>
    fnIntegrate_x(0, 1, NCALLS_SET, True)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 35, in fnIntegrate_x
    I = Integrate_x(yarray)
  File "C:\Users\JP\Documents\Python\TestingQuad\TestingQuad_v2.py", line 23, in Integrate_x
    return quad(Integrand, 0, np.pi/2, args=(y))[0]
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 247, in quad
    retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)
  File "C:\Python27\lib\site-packages\scipy\integrate\quadpack.py", line 312, in _quad
    return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
quadpack.error: Supplied function does not return a valid float.

'''
4

4 に答える 4

2

ここで自分の質問に否定的に答えているのではないかと心配しています。私はそれが可能だとは思わない。クワッドは、何か他のもので書かれたライブラリのある種のポートのようです - そのように、物事がどのように行われるかを定義するのは内部のライブラリです - そのため、ライブラリ自体を再設計しない限り、私が望んでいたことを行うことはおそらく不可能です.

複数の D 統合でタイミングの問題を抱えている他の人にとっては、専用の統合ライブラリを使用するのが最善の方法であることがわかりました。「cuba」にはかなり効率的なマルチ D 統合ルーチンがいくつかあるようです。

http://www.feynarts.de/cuba/

これらのルーチンは c で書かれているので、最終的には SWIG を使用してそれらと対話することになりました - そして最終的に効率のために被積分関数を c で書き直しました - これにより負荷が高速化されました....

于 2013-03-03T18:02:54.440 に答える
2

numpy.vectorize 関数を介して可能です。私は長い間この問題を抱えていましたが、このベクトル化機能にたどり着きました。

次のように使用できます。

vectorized_function = numpy.vectorize(your_function)

output = vectorized_function(your_array_input)
于 2015-06-27T19:48:32.123 に答える
0

すべての次元で -np.inf から np.inf までの確率密度関数を統合する際に、この問題が発生していました。

*args を取り込むラッパー関数を作成し、args を numpy 配列に変換し、ラッパー関数を統合することで修正しました。

numpy の vectorize を使用すると、すべての値が等しい部分空間のみが統合されると思います。

次に例を示します。

from scipy.integrate import nquad
from scipy.stats import multivariate_normal

mean = [0., 0.]
cov = np.array([[1., 0.],
                [0., 1.]])

bivariate_normal = multivariate_normal(mean=mean, cov=cov)

def pdf(*args):
    x = np.array(args)
    return bivariate_normal.pdf(x)

integration_range = [[-18, 18], [-18, 18]]

nquad(pdf, integration_range)

Output: (1.000000000000001, 1.3429066352690133e-08)
于 2021-02-08T01:14:33.183 に答える