The grad function in the {pracma} and {numDeriv} R libraries gives erroneous results

Question

I am interested in the first-order numerical derivative of the self-determined function pTgh_y(q,g,h) respect to q. For a special case, pTgh_y (q, 0,0) = pnorm (q). In other words, pTgh_y(q,g,h) reduces to the CDF of the standard normal when g = h = 0 (see Figure below).

This means that d pTgh_y (0,0,0) / dq should correspond to the following

 dnorm(0) 

0.3989423

 grad(pnorm,0) 

0.3989423

Here are some of my attempts with the grad function in the {pracma} library.

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

0

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

0

Here are some of my attempts with the grad function in the {numDeriv} library.

 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): invalid input type Error in method grad.default (function (x) {: function does not accept a complex argument, as required by the "complex" method.

None of these functions give the correct result.

The function My pTgh_y(q,g,h) is defined as follows

 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)/gq[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

Your function is flat around 0, so calculating the derivative of 0 is true:

 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 

You need to increase the accuracy of your pTgh_y function:

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

And now you get the result you want:

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