векторизация неполной гамма-функции с отрицательными аргументами - PullRequest
Новая
       95

векторизация неполной гамма-функции с отрицательными аргументами

3 голосов
/ 23 января 2020

Я использую неполную гамма-функцию mpmaths с отрицательным z и комплексные границы интегрирования a, b. См. Документацию: документы

Этого нельзя сделать с помощью неполной гамма-функции Сципи. Я сделал отметку mpmath против Mathematica и получил те же результаты. К сожалению, функция mpmaths может оценивать только одно значение за раз, что означает, что мне нужно l oop по массиву, который содержит значения z, a, b.

Я считаю, что векторизация трудна, потому что для отрицательного z нужно рекурсивно оценивать интеграл по частям (см. Этот пост и соответствующую вики: неполная гамма-функция в python? )

Кто-нибудь знает модуль, который может справиться с этим? Или если такой вещи нет, если это вообще возможно?

По сути: 

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)

Я не знаю, насколько это хорошо это хотя. Таким образом, вопрос немного меняется, как я могу оценить и / или улучшить это? Кроме того, это, кажется, работает для фиксированных границ. Однако я не знаю, как можно реализовать интеграл с бесконечными границами на обоих концах.

...