如何计算python中矩母函数的导数？

Question

到目前为止，这是我的代码，我认为我可以使用 scipy，但它没有给出二阶导数 moment(0, 2) 的正确答案。我的猜测是我没有正确应用 scipy.misc.derivative 并且我应该使用 sympy 的 diffs_exp 但我也无法让它工作..

from scipy import misc
import numpy as np

def mgf(s):
    mu = 2
    sigma = 0.5
    mgf = np.exp(mu*s + ((sigma**2)*(s**2))/2)
    return mgf

def moment(s, i): 
    mo = misc.derivative(mgf, s, dx=0.000000001, n=i)
    return mo

moment(s, i) 在 i=1 时正确计算，但在 i>1 时不正确。 moment(0,2) 应该等于 sigma^2 或 .25 但函数 returns 0.0 当前

函数只有在s=0时才会求值，更重要的是微分正确。

Answer 1

所以我在误读了早期答案（已删除）中的问题后，在源代码中乱搞。

Scipy misc.derivative默认计算二阶导数为

lim h->0 (f(x+h)-2*f(x)+f(x-h))/h^2

这里出现问题是因为 np.exp() 的输出是 float64，其精度有限，即尾数为 52 位，指数为 11 位。当我们减少 dx 时，项之间的差异出现在高位数字中，由于精度有限，这些数字不存在。总而言之，这消失为零。

供参考，在上述函数中，值为对于 f(x+h)、f(x)、f(x-h) 分别为 0.999999998、1、1.000000002，其中 x=0 且 h =1e-9。解决方案可以使用具有更高精度的任一函数，但这将涉及更改 scipy 源代码并且不是一件小事。（Python pow 函数不是任意精度）

其他（实用）选项是使用较小的 dx 值。 dx=1e-2 似乎实际上给出了足够接近的答案，即 4.2501848996 与实际二阶导数 4.25 相比。

Answer 2

为有限差分方案选择合适的步长是一项棘手的工作。由于舍入误差（如您所见），步幅太小您注定失败。步骤太大，方案太粗糙（正如您所发现的那样）。 scipy.misc.derivative 的默认步骤不是很有用，顺便说一句。有一些关于如何选择合理步骤的文献。例如，Numerical Recipes有一个简单方案的简单介绍。

在这种特殊情况下，找到合理的步骤相当容易：

In [41]: from scipy.misc import derivative

In [42]: def f(x):
   ....:     arg = 2.*x + (0.5*x)**2 / 2.
   ....:     return np.exp(arg)
   ....: 

In [53]: derivative(f, 0., dx=1e-5, n=2)
Out[53]: 4.2499981312005266

另一种方法是使用一个包来进行更智能的步长选择（internet/literature 搜索的一个关键字是 Romberg 外推法）。例如，numdifftools:

In [57]: import numdifftools as nd

In [59]: fdd = nd.Derivative(f, n=2)

In [60]: fdd(0)
Out[60]: array([ 4.25])

Answer 3

以下是如何使用 sympy 进行象征性的操作，并以数值方式评估特定 mu、sigma 和 s

的结果

In [1]: from sympy import *

In [2]: mu, sigma, s = symbols("mu sigma s")

In [3]: expr = exp(mu*s+(sigma*s)**2/2)

In [4]: f = lambdify((mu, sigma, s), expr.diff(s, 2))

In [5]: f(2, 0.5, 0)
Out[5]: 4.25

如何计算python中矩母函数的导数？

How to calculate derivative of moment generating function in python?

python

statistics

sympy

scipy