積分関数について

[masakicktashiro]

scipyのintegrate.nquadに値する関数を作ろうと思っているのですがメモリの問題と 精度の問題から思うような関数が作れていません。もし僕の考えたコードもあげるので もしよかったら改善したコードを作っていただきたいです。 ちなみに1個目はmeshgridを使った実装で2個目は乱数で実装です。

def integral(_func, start, end, n_split=500, limit=1000, args=None):
    #範囲が無限の時に1000のリミットをつけ、幅も細かく
    if np.isinf(start):
        start = limit if start > 0 else -limit
        n_split = 10000
    if np.isinf(end):
        end = limit if end > 0 else -limit
        n_split = 10000
    interval = (end - start) / n_split
    return sum([_func(i,*args)*interval for i in np.linspace(start, end, n_split)]) \
        if args is not None else um([_func(i)*interval for i in np.linspace(start, end, n_split)])

def limit_inf(num, limit=1000):
    if not np.isinf(num):
        return num
    else:
        return limit if num > 0 else -limit

def n_dimension_integral(_func, ranges, n_split=500, args=None):
    """
    ranges...([low,high],[low,high],...)
    """
    ranges = tuple(list(map(limit_inf ,_range)) for _range in ranges)
    interval = np.prod([(_range[1] - _range[0])/ n_split for _range in ranges], axis=0)
    ranges = tuple(list(map(lambda x: np.linspace(*x, n_split), ranges)))
    combinations = np.asarray([_range.ravel() for _range in np.meshgrid(*ranges)]).T
    return sum([_func(*combi,*args)*interval for combi in combinations]) \
        if args is not None else \
            sum([_func(*combi)*interval for combi in combinations])

def n_intergral(_func, ranges, num_random=10000, args=None):
    """
    ranges...([low,high],[low,high],...)
    """
    ranges = tuple(list(map(limit_inf ,_range)) for _range in ranges)
    interval = np.prod([(_range[1] - _range[0]) for _range in ranges], axis=0) / num_random
    combinations = np.asarray([
        np.random.uniform(*_range,num_random) for _range in ranges]).T
    print(interval,combinations.shape)
    return sum([_func(*combi,*args)*interval for combi in combinations]) \
        if args is not None else \
            sum([_func(*combi)*interval for combi in combinations])

[yasudadesu]

あこれissueじゃなくてPRで出していただけるとって感じですね、 皆さんでコードに対して意見だしあうならそちらの方がいいと思います.