Chức năng grad trong cả thư viện {pracma} và {numDeriv} của R cho kết quả sai lệch

Question

Tôi quan tâm đến đạo hàm bậc 1 của hàm tự định nghĩa pTgh_y(q,g,h) đối với q. Đối với trường hợp đặc biệt, pTgh_y (q, 0,0) = pnorm (q). Nói cách khác, pTgh_y(q,g,h) được giảm xuống CDF của chuẩn bình thường khi g = h = 0 (xem hình bên dưới).

Điều này có nghĩa rằng d pTgh_y (0,0,0)/dq phải bằng những điều sau đây

dnorm(0) 

0,3989423

grad(pnorm,0) 

0,3989423

Dưới đây là một số nỗ lực của tôi với hàm grad trong thư viện {pracma}.

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) 

Dưới đây là một số nỗ lực của tôi với chức năng grad trong {} numDeriv thư viện.

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') 

Lỗi trong Pmin (x, .Machine $ double.xmax): inp không hợp lệ ut type Lỗi trong hàm grad.default (hàm (x) {: không chấp nhận đối số phức tạp theo yêu cầu của phương thức 'phức tạp'.

Không có chức năng nào trong số này cung cấp kết quả chính xác.

chức năng My pTgh_y(q,g,h) được định nghĩa như sau

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

Chức năng của bạn là phẳng khoảng 0, do đó tính toán một dẫn xuất của 0 là đúng:

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 

Bạn cần phải tăng độ chính xác của chức năng của bạn 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) 
    } 
} 

Và bây giờ bạn nhận được kết quả mà bạn muốn có:

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