数值Python - 我如何使这成为一个ufunc？

Question

NumPy的新手可能没有正确搜索,所以如果这是一个常见问题,我会采取肿块...

我正在研究一个问题,我需要计算相对较大数字的log(n!) - 即.为了大计算阶乘,所以我写了以下函数:

def log_fact(n):
    x = 0
    for i in range(1,n+1):
        x += log(i)
    return x

现在问题是我想将此作为传递给curve_fit的函数的一部分:

def logfactfunc(x, a, b, c):
    return a*log_fact(x) + b*x + c

from scipy.optimize import curve_fit

curve_fit(logfactfunc, x, y)

但是,这会产生以下错误:

File "./fit2.py", line 16, in log_fact
    for i in range(1,n+1):
TypeError: only length-1 arrays can be converted to Python scalars

一点点搜索建议numpy.frompyfunc()将其转换为ufunc

curve_fit(np.frompyfunc(logfactfunc, 1, 1), data[k].step, data[k].sieve)

TypeError: <ufunc 'logfactfunc (vectorized)'> is not a Python function

试过这个:

def logfactfunc(x, a, b, c):
    return a*np.frompyfunc(log_fact, 1, 1)(x) + b*x + c

File "./fit2.py", line 30, in logfactfunc
    return a*np.frompyfunc(log_fact, 1, 1)(x) + b*x + c
TypeError: unsupported operand type(s) for +: 'numpy.ndarray' and 'numpy.float64

关于如何在curve_fit()函数中使用log_fact()函数的任何想法？

谢谢!

Answer 1

您的log_fact函数与函数密切相关,gammaln函数定义为ufunc scipy.special.具体来说,log_fact(n) == scipy.special.gammaln(n+1).对于即使适度的值n,这明显更快:

In [15]: %timeit log_fact(19)
10000 loops, best of 3: 24.4 us per loop
In [16]: %timeit scipy.special.gammaln(20)
1000000 loops, best of 3: 1.13 us per loop

而且,运行时间gammaln是独立的n,不像log_fact.