首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >加速scipy.integrate.quad()?

加速scipy.integrate.quad()?
EN

Stack Overflow用户
提问于 2013-06-20 02:42:08
回答 2查看 2.1K关注 0票数 1

我正在尝试加速下面的计算积分和的代码。为了获得更好的准确性,我需要增加L_max,但这也会使执行时间变得更长。下面的具体情况计算0.999999的概率曲线,大约需要65秒。我听说过cython及其加速代码的能力,但我不知道如何使用它,也不知道它在这种情况下能有什么帮助。有什么想法吗?

代码语言:javascript
复制
import math
from scipy import integrate
import numpy
from decimal import *
import time

start_time=time.time()
getcontext().prec=100
################################
def pt(t):
    first_term=math.exp(-lam*t)*((0.0001*I)**ni)*(t**(ni-1))*math.exp(-(0.0001*I)*t)/(math.factorial(ni-1))
    sum_term=0.0
    i=0
    while i<ni:
        sum_term=sum_term+((0.0001*I)**i)*(t**(i))*math.exp(-(0.0001*I)*t)/(math.factorial(i))
        i=i+1
    sum_term=lam*math.exp(-lam*t)*sum_term
    total=first_term+sum_term
    return total
#################################
def pLgt(t):
    return Decimal(((Decimal((0.0001*O)*t))**Decimal(L))*Decimal(math.exp(-(0.0001*O)*t)))/Decimal((math.factorial(L)))
######################################
def pL_t(t):
    return (pLgt(t))*Decimal(pt(t))
################################
lam=0.0001
beta=0.0001
ni=10
I=5969
O=48170
L_max=300.0
L=0.0
sum_term=0.0
sum_probability=0.0
while L<L_max:
    probability=(integrate.quad(lambda t: pL_t(t),0,800))[0]
    sum_probability=sum_probability+probability
    sum_term=sum_term+L*probability
    L=L+1.0
print time.time()-start_time
print sum_probability
print sum_term
print (sum_term-1)*0.46+6.5 
EN

回答 2

Stack Overflow用户

发布于 2013-06-20 03:10:40

用Decimal进行计算可能会大大减慢你的速度,而不会带来任何好处。小数计算比浮点数慢得多,大约慢100倍,正如Kozyarchuk在Stack Overflow here上所指出的那样。在Numpy数组中使用Decimal类型使您无法从Numpy获得速度优势。

同时,我不清楚scipy.integrate.quad的结果是否真的达到您想要的精度水平;如果您确实需要任意精度,您可能必须从头开始编写正交代码。

如果您确实需要使用Decimal数字,那么至少缓存这些数字的Decimal表示将为您提供一些速度优势。也就是说,使用

代码语言:javascript
复制
O=Decimal(48170)
L=Decimal(0.0)

告诉pLgt只使用O和L会更快。

票数 2
EN

Stack Overflow用户

发布于 2013-06-20 03:11:31

pL_t函数看起来像伽马分布的总和,在这种情况下,您应该能够将积分计算为部分不完整伽马函数的总和:http://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gdtr.html#scipy.special.gdtr

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/17198898

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档