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

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 

Но для этого точно таких же функций я использую вольфрам и получить -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