python - Numerical 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 函数密切相关，后者在 scipy.special 中定义为一个 ufunc。具体来说，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 不同。