La funzione grad in entrambe le librerie {pracma} e {numDeriv} di R fornisce risultati errati

Question

Sono interessato alla derivata numerica del 1 ° ordine di una funzione auto-definita pTgh_y(q,g,h) rispetto a q. Per un caso speciale, pTgh_y (q, 0,0) = pnorm (q). In altre parole, pTgh_y(q,g,h) viene ridotto al CDF della norma standard quando g = h = 0 (vedere l'immagine sotto).

Ciò significa che d pTgh_y (0,0,0)/dq deve essere uguale al seguente

dnorm(0) 

0,3989423

grad(pnorm,0) 

0,3989423

Ecco alcuni dei miei tentativi con la funzione grad nella libreria {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) 

Ecco alcuni dei miei tentativi con la funzione grad nella {} numDeriv biblioteca.

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

Error in Pmin (x, .Machine $ double.xmax): inp non valida ut tipo Errore in grad.default (la funzione function (x) {: non accetta argomenti complessi come richiesto dal metodo 'complex'.

Nessuna di queste funzioni fornisce risultati corretti.

La mia funzione pTgh_y(q,g,h) è definito come segue

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

La vostra funzione è piatta circa 0, calcolando così un derivato di 0 è corretta:

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 

è necessario aumentare la precisione della funzione 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) 
    } 
} 

E ora ottieni il risultato che volevi avere:

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