间隔-x - >无穷大的高斯 - 勒让德:有效地变换权重和节点的自适应算法

Question

好吧我知道之前有一个有限的例子已经被要求缩放[-1, 1]间隔[a, b] 高斯 - 勒让德正交在numpy中的不同间隔但是没有人发布如何概括它[-a, Infinity](如下所述,但不是(还)快).这也显示了如何使用多个实现调用复杂函数(无论如何在定量选项定价中).有基准quad代码,后面跟着leggauss代码示例链接到如何实现自适应算法.我已经解决了大多数链接的adaptive algorithm困难 - 它目前打印分割积分的总和以显示它正常工作.在这里,你会发现功能的范围从转换[-1, 1]到[0, 1]以[a, Infinity](感谢@AlexisClarembeau).要使用自适应算法,我不得不创建另一个函数来转换[-1, 1]到[a, b]被反馈到[a, Infinity]功能.

import numpy as np
from scipy.stats import norm, lognorm
from scipy.integrate import quad

a = 0
degrees = 50

flag=-1.0000
F = 1.2075
K = 0.1251
vol = 0.43
T2 = 0.0411
T1 = 0.0047

def integrand(x, flag, F, K, vol, T2, T1):
    d1 = (np.log(x / (x+K)) + 0.5 * (vol**2) * (T2-T1)) / (vol * np.sqrt(T2 - T1))
    d2 = d1 - vol*np.sqrt(T2 - T1)
    mu = np.log(F) - 0.5 *vol **2 * T1
    sigma = vol * np.sqrt(T1)
    return lognorm.pdf(x, mu, sigma) * (flag * x*norm.cdf(flag * d1) - flag * (x+K)*norm.cdf(flag * d2))

def transform_integral_0_1_to_Infinity(x, a): 
    return integrand(a+(x/(1-x)), flag, F, K, vol, T2, T1) *(1/(1-x)**2); 

def transform_integral_negative1_1_to_0_1(x, a): 
    return 0.5 * transform_integral_0_1_to_Infinity((x+1)/2, a)

def transform_integral_negative1_1_to_a_b(x, w, a, b):
    return np.sum(w*(0.5 * transform_integral_0_1_to_Infinity(((x+1)/2*(b-a)+a), a)))

def adaptive_integration(x, w, a=-1, b=1, lastsplit=False, precision=1e-10):
    #split the integral in half assuming [-1, 1] range
    midpoint = (a+b)/2
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint)
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b)
    return interval1+interval2 #just shows this is correct for splitting the interval

def integrate(x, w, a): 
    return np.sum(w*transform_integral_negative1_1_to_0_1(x, a))

x, w = np.polynomial.legendre.leggauss(degrees) 
quadresult = quad(integrand, a, np.Inf, args=(flag, F, K, vol, T2, T1), epsabs=1e-1000)[0]
GL = integrate(x, w, a)
print("Adaptive Sum Result:")
print(adaptive_integration(x, w))
print("GL result"); 
print(GL)
print("QUAD result")
print(quadresult)

仍然需要以更小的尺寸来提高速度和精度,因为我无法手动调整degrees范围-a以获得收敛.为了说明为什么这是一个问题,把这些值改为:a=-20,F=50,然后运行.degrees=1000如果没有智能应用,您可以增加并看到此Gauss-Legendre算法没有任何好处.我对速度的要求是每循环达到0.0004s,而最后一次算法我Cython化大约需要0.75s,这就是我尝试使用Gauss-Legendre的低度,高精度算法的原因.使用Cython和多线程,来自完全优化的Python实现的这个要求是每个循环大约0.007s(非向量化,循环,低效的例程,每个循环可以是0.1s degrees=20,即%timeit adaptive_integration(x,w).

我已经半实现的一个可能的解决方案是在第5/6页的http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals,adaptive integration而间隔a-b(在这种情况下,我编写了transform_integral_negative1_1_to_a_b函数),其中间隔被划分在2(@ 0.5)中,然后在这1/2个间隔上评估函数,并将两个0->0.5+的0.5->1和与整个范围的函数结果进行比较0->1.如果精度不在公差范围内,则范围进一步细分为0.25和0.75,再次对每个子区间评估函数,并与先前的1/2区间和@进行比较0.5.如果1侧在公差范围内(例如abs(0->0.5 - (0->0.25 + 0.25->0.5)) < precision),但另一侧不在公差范围内,则分裂在公差范围内停止,但在另一侧继续,直到precision达到.此时,将每个间隔切片的结果相加,以获得具有更高精度的完整积分.

有可能更快,更好的方法来解决这个问题.只要它快速准确,我就不在乎.以下是我参考的集成例程的最佳描述http://orion.math.iastate.edu/keinert/computation_notes/chapter5.pdf奖励是100分赏金+ 15分答案接受.感谢您协助快速准确地完成此代码!

编辑:

以下是我对adaptive_integration代码的更改- 如果有人能够快速完成这项工作,我可以接受答案并奖励赏金.第7页的这个Mathematica代码http://online.sfsu.edu/meredith/Numerical_Analysis/improper_integrals执行我尝试的例程.它有一个不能很好地收敛的例程,请参阅下面的变量.现在我的代码错误了:RecursionError: maximum recursion depth exceeded in comparison在一些输入上,或者如果degrees设置得太高,或者quad在它工作时没有接近结果,那么这里显然是错误的.

def adaptive_integration(x, w, a, b, integralA2B, remainingIterations, firstIteration, precision=1e-9):
    #split the integral in half assuming [-1, 1] range
    if remainingIterations == 0:
        print('Adaptive integration failed on the interval',a,'->',b)
    if np.isnan(integralA2B): return np.nan

    midpoint = (a+b)/2
    interval1 = transform_integral_negative1_1_to_a_b(x, w, a, midpoint)
    interval2 = transform_integral_negative1_1_to_a_b(x, w, midpoint, b)
    if np.abs(integralA2B - (interval1 + interval2))  < precision : 
        return(interval1 + interval2)       
    else:
        return adaptive_integration(x, w, a, midpoint, interval1, (remainingIterations-1), False) + adaptive_integration(x, w, midpoint, b, interval2, (remainingIterations-1), False) 

#This example doesn't converge to Quad

# non-converging interval inputs
a = 0 # AND a = -250
degrees = 10

flag= 1
F = 50
K = 0.1251
vol = 0.43
T2 = 0.0411
T1 = 0.0047

print(adaptive_integration(x, w, -1, 1, GL, 500, False))

与输出degrees=100(计算后GL用degrees=10000一个更好的初始估计,否则,该算法总是与它自己的精度显然同意和不调用从而未能每次自适应路径):

GL result:
60.065205169286379
Adaptive Sum Result:
RecursionError: maximum recursion depth exceeded in comparison
QUAD result:
68.72069173210338

Answer 1

我认为代码完成了这项工作:

import numpy as np 
import math

deg = 10
x, w = np.polynomial.legendre.leggauss(deg)

def function(x): 
    # the function to integrate
    return math.exp(-x)

def function2(x, a): 
    return function(a+x/(1-x))/((1-x)**2); 

def anotherOne(x, a): 
    return 0.5 * function2(x/2 + 1/2, a)

def integrate(deg, a): 
    sum = 0 
    x, w = np.polynomial.legendre.leggauss(deg)
    for i in range(deg): 
        print("sum({}) += {} * {} (eval in {})".format(sum, w[i], anotherOne(x[i], a), x[i]))
        sum += w[i]*anotherOne(x[i], a)
    return sum; 

print("result"); 
print(integrate(10, 1))

它结合了你的方程,从a到inf和方程的积分来改变积分的界限.

我希望它能解决你的问题(至少适用于exp(-x)):)

如果您想要内联计算,程序将总和:

它是以下组合:

和:

和: