Большая ошибка при численном вычислении второй производной

Question

Я пытаюсь численно взять вторую производную функции l (лог распределения Пуассона по вектору x а также lambda=6) в R и это мой код:

    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 для первой производной и -5/36~~-0.138889 для второй производной. Я пытался выяснить, почему моя функция имеет такую ​​большую ошибку за последние два часа.

Примечание: это для проекта класса, поэтому я не могу использовать производную функцию в R.

Answer 1

Я думаю, что проблемы, которые вы видите, связаны с округлением и другими числовыми проблемами. Я бы предложил два подхода:

1. Увеличьте свой delta а также h ценности; так как ваш компьютер имеет конечную точность, во многих случаях t а также t+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