大誤差的二階導數

Question

我試圖採取數值的函數l在R（泊松分佈在載體x和lambda=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 

但我用鎢，並得到-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