3

我正在使用 mpmaths 不完整的 gamma 函数，z 为负，复杂的积分范围为 a，b。请参阅文档：文档

这不能用 scipy 的不完整 gamma 函数来完成。我将标记的 mpmath 放在 Mathematica 的替补席上，得到了相同的结果。不幸的是，mpmaths 函数一次只能评估一个值，这意味着我必须遍历包含 z、a、b 值的数组。

我相信矢量化很难，因为对于负 z，必须递归地评估积分部分（参见这篇文章和相应的 wiki：python 中的不完整 gamma 函数？）

有谁知道可以处理这个的模块？或者，如果它甚至可能的话，就没有这样的事情吗？

所以本质上：

from mpmath import gammainc
values = np.random.normal(0, 100, 1000) + 1j*np.random.normal(0, 100, 1000)
res_mp = np.array([gammainc(cc, 4, 10) for cc in values])

我实际上已经实现了一些我认为在某种程度上有效的东西：

import numpy as np

def new_quad_routine(func, c, a, b,x_list, w_list):
    c_1 = (b-a)/2.0
    c_2 = (b+a)/2.0

    eval_points = c_1*x_list+c_2             
    func_evals = func(eval_points, c)                              
    return c_1 * np.sum(func_evals * w_list[:,np.newaxis], axis=0)

def new_quad_gauss_7(func, c, a, b):   
    """
    integrates a complex function with bounds, a, b with an array of arguments c
    call for instance with the gamma function:
    new_quad_gauss_7(gamma_integrator, 4, 10)

    """
    x_gauss = np.array([-0.949107912342759, -0.741531185599394, -0.405845151377397, 0, 0.405845151377397, 0.741531185599394, 0.949107912342759])            
    w_gauss = np.array([0.129484966168870, 0.279705391489277, 0.381830050505119, 0.417959183673469, 0.381830050505119, 0.279705391489277,0.129484966168870])
    return new_quad_routine(func, c, a, b, x_gauss, w_gauss)

def gamma_integrator(t, c):
    return t[:, np.newaxis]**c*np.exp(-t[:,np.newaxis])

def gammainc_function(c, a, b):
    if type(c) is not np.ndarray:
        raise ValueError("Need a numpy array for function argument!")
    return new_quad_gauss_7(gamma_integrator, c-1, a, b)

这是基于本文中的数值积分器。

它非常快：

values=np.repeat(-2+3j-1, 100)
In [143]: %timeit new_quad_gauss_7(gamma_integrator, 4, 10, values)
107 µs ± 574 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [144]: %timeit [gammainc(cc, 4, 10) for cc in values]
224 ms ± 2.86 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

我不知道它有多好。所以这个问题稍微改变了我如何进行基准测试和/或改进它？此外，这似乎适用于固定范围。但是，我不知道如何实现两端无限边界的积分。

4

0 回答 0