如何将sympy导数与numpy数组表达式相结合？

Question

我希望我能把我的问题讲清楚。在

我用scipy.optimize.minimize来寻找一些实验数据的最大似然估计量。数据是2d（NxN），因此在我的注释中，对数似然函数（我需要最大化的函数）如下所示：

{1美元^

因此是一个2d数组，它依赖于我的特定模型，而n_ij是测量的二维数据。在代码中，我用这种方式编写了相同的函数：

def logll(params, *args):

    A, x0, y0, bkg = params
    pico, F = args

    x, y = np.arange(pico.shape[0]), np.arange(pico.shape[1])

    erfi = erf((x + 1 - x0) / F) - erf((x - x0) / F)
    erfj = erf((y + 1 - y0) / F) - erf((y - y0) / F)

    lambda_p = A * F**2 * np.pi * erfi[:, np.newaxis] * erfj / 4 + bkg

    return - np.sum(pico * np.log(lambda_p) - lambda_p)

这里需要注意的一点是，我的数据是2d数组的一部分，因此乘积erfi * erfj和{}是{}数组。我所知道的编写代码的唯一方法是如上所述，即将erfi和{}写成1d数组，并将它们乘以erfi[:, np.newaxis] * erfj以得到正确的形状数组。在

整个过程都很顺利，我只要打电话就能找到我的估计器

^{pr2}$

总之，这个Powell方法不够快，所以我想用一个像Newton-CG这样的方法，它应该更快。缺点是它需要jacobian和hessian，为此，我必须区分logll关于我的4个参数。我不想用手来做，所以我用sympy试试。我不得不从表达式中去掉[:, np.newaxis]部分，但其他的是导数没问题。在

jacobian ll_jac是包含{}相对于四个拟合参数的导数的一维数组。我把A、x0、y0、bkg、x和{}定义为带有{}的符号，然后用这种方式进行区分：

logll.diff(A)
logll.diff(x0)
logll.diff(y0)
logll.diff(bkg)

hessian ll_hess是一个2d数组，它包含logll关于四个参数的二阶导数，我通过

logll.diff(A, A)
logll.diff(A, x0)

等等。在这里，在做任何微分之前，我必须从logll表达式中去掉[:, np.newaxis]部分，以便sympy工作。在

所以，最后，有了雅各布和黑森，我试着

minimize(logll, x0=[A, x0, y0, bkg], args=(peak, F), 
         method='Newton-CG', jac=ll_jac, hess=ll_hess)

我得到了以下错误：

File "/usr/lib/python2.7/site-packages/scipy/optimize/_minimize.py", line 351, in minimize
**options)

File "/usr/lib/python2.7/site-packages/scipy/optimize/optimize.py", line 1365, in _minimize_newtoncg
update = alphak * pk

TypeError: unsupported operand type(s) for *: 'NoneType' and 'float'

我怀疑这与erfi和erfj是1d数组，bkg是2d数组有关。在

所以我的问题是：我需要如何写下这些表达式，以便numpy和{}不会发生冲突？有没有更好的方法来编写乘积erfi * erfj并且仍然得到我想要的2d数组？在