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

Я пытаюсь численно взять вторую производную функции 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.


r poisson
person pretzelman    schedule 18.10.2015    source источник

Ответы (1)


arrow_upward
1
arrow_downward

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

  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
person josliber♦    schedule 18.10.2015