R的{pracma}和{numDeriv}庫中的grad函數給出了錯誤的結果

Question

我對自定義函數pTgh_y(q,g,h)關於q的一階數值導數感興趣。對於特殊情況，pTgh_y（q，0,0）= pnorm（q）。換句話說，當g = h = 0時（參見下圖），pTgh_y(q,g,h)減少到標準正態分佈的CDF。

這意味着，d pTgh_y（0,0,0）/ dQ的應該等於以下

dnorm(0) 

0.3989423

grad(pnorm,0) 

0.3989423

下面是我在{pracma}庫中使用grad函數的一些嘗試。

library(pracma) 
# load pTgh and all relevant functions 
grad(function(x){pTgh_y(x,0,0)},0) 

grad(function(x){pTgh_y(x,0,0)},0,heps=1e-10) 

下面是我的一些嘗試與在{} numDeriv庫中畢業的功能。

library(numDeriv) 
# load pTgh and all relevant functions 
grad(function(x){pTgh_y(x,0,0)},0,method='simple') 

0.3274016

grad(function(x){pTgh_y(x,0,0)},0,method='Richardson') 

-0.02505431

grad(function(x){pTgh_y(x,0,0)},0,method='complex') 

在PMIN錯誤（X，$。機double.xmax）：INP無效ut類型 grad.default中的錯誤（function（x）{： 函數不接受方法'complex'所要求的複雜參數。

這些函數都沒有給出正確的結果。

我pTgh_y(q,g,h)函數的定義如下

qTgh_y = function(p,g,h){        
    zp = qnorm(p) 
    if(g==0) q = zp          
    else  q = (exp(g*zp)-1)*exp(0.5*h*zp^2)/g  
    q[p==0] = -Inf 
    q[p==1] = Inf 
    return(q) 
} 

pTgh_y = function(q,g,h){      
    if (q==-Inf) return(0) 
    else if (q==Inf) return(1) 
    else { 
    p = uniroot(function(t){qTgh_y(t,g,h)-q},interval=c(0,1)) 
    return(p$root) 
    } 
} 

Answer 1

你的功能是平在0附近，從而計算的0衍生是正確的：

f = function(x){pTgh_y(x,0,0)} 
h = 0.00001; f(0+h); f(0-h) 
# [1] 0.5 
# [1] 0.5 

library(pracma) 
grad(f, 0, heps=1e-02); grad(f, 0, heps=1e-03); 
grad(f, 0, heps=1e-04); grad(f, 0, heps=1e-05) 
# [1] 0.3989059 
# [1] 0.399012 
# [1] 0.4688766 
# [1] 0 

你需要增加你的函數的精確度pTgh_y：

pTgh_y = function(q,g,h){      
    if (q==-Inf) return(0) 
    else if (q==Inf) return(1) 
    else { 
     p = uniroot(function(t){qTgh_y(t,g,h)-q},interval=c(0,1), 
        tol = .Machine$double.eps) 
     return(p$root) 
    } 
} 

現在你得到你想要的結果：

f = function(x){pTgh_y(x,0,0)} 
grad(f, 0) 
[1] 0.3989423