数值计算二阶导数误差大

Large error when numerically computing second derivative

我正在尝试在 R 中对函数 l(泊松分布在向量 xlambda=6 上的对数)进行数值计算,这是我的代码:

    x=c(2,3)
    t=6
    delta=1e-12
    h=1e-12

    L=function(x,t) dpois(x,t)
    l<-function(x,t) log(prod(L(x,t)))
    ld<-function(x,t) (l(x,t+delta)-l(x,t))/delta
    ldd<-function(x,t) (ld(x,t+h)-ld(x,t))/h
    ld(x,t)
    ldd(x,t)

我的输出是

> ld(x,t)
[1] -1.167066
> ldd(x,t)
[1] 888178420

但是对于这个完全相同的函数,我使用 wolfram 并得到 -7/6~~-1.16667 for the first derivative and -5/36~~-0.138889 作为二阶导数。在过去的两个小时里,我一直在试图弄清楚为什么我的函数有这么大的错误。

注意:这是一个 class 项目,所以我不能在 R 中使用导数函数。

我认为您看到的问题是由于四舍五入和其他数值问题造成的。我建议两种方法:

  1. 提高您的 deltah 值;由于您的计算机精度有限,因此在许多情况下 tt+1e-12 将完全相同。当然,这在数值计算导数时会造成很大的问题。
  2. 为了数值稳定性,将 log(prod(x)) 替换为 sum(log(x)) 通常是个好主意。

通过这两个调整,我们得到了更好的结果:

delta = 1e-5
h = 1e-4
L=function(x,t) dpois(x,t)
l<-function(x,t) sum(log(L(x,t)))
ld<-function(x,t) (l(x,t+delta)-l(x,t))/delta
ldd<-function(x,t) (ld(x,t+h)-ld(x,t))/h
ld(x,t)
# [1] -1.166667
ldd(x,t)
# [1] -0.1388853