[R] Numerical Derivatives in R

[Paul Gilbert]

(This code looks vaguely familiar.)

Is anyone interested in participating in an effort to make a self 
contained package for numerical derivatives? I would be happy to extract 
the Richardson extrapolation code for first and second derivatives from 
my package in the devel area of CRAN, but I'm a bit too busy to spend 
much time on it myself right now. (Also, one thing missing is good 
documentation, and I find it helps to have more than one person look at 
the documentation.)

Paul Gilbert

Tolga Uzuner wrote:
 > Actually, I did implement this using richardson extrapolation, but am 
having trouble vectorising it. For some reason, it fails within integrate...
 >
 > Anyone willing to look over the below and let me know what I am doing 
wrong, helps much appreciated. You can cut paste the below into the 
console..
 >
 > XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 >
 > richardson.grad <- function(func, x10, d=0.01, eps=1e-4, r=6, show=F){
 > sapply(x10,function(x){
 >
 >   v <- 2               # reduction factor.
 >   n <- length(x)       # Integer, number of variables.
 >   a.mtr <- matrix(1, r, n)
 >   b.mtr <- matrix(1, (r - 1), n)
 >
 >   h <- abs(d*x)+eps*(x==0.0)
 >
 >   for(k in 1:r)  { # successively reduce h                    for(i 
in 1:n)  {
 >          x1.vct <- x2.vct <- x
 >          x1.vct[i]  <- x[i] + h[i]
 >          x2.vct[i]  <- x[i] - h[i]
 >          if(k == 1) a.mtr[k,i] <- (func(x1.vct) - func(x2.vct))/(2*h[i])
 >          else{
 >            if(abs(a.mtr[(k-1),i])>1e-20)
 >                  # some functions are unstable near 0.0 
                 a.mtr[k,i] <- (func(x1.vct)-func(x2.vct))/(2*h[i])
 >             else  a.mtr[k, i] <- 0
 >           }
 >       }
 >      h <- h/v     # Reduced h by 1/v.
 >     }
 >    if(show) {
 >
 >         cat("\n","first order approximations", "\n") 
print(a.mtr, 12)
 >     }
 >   for(m in 1:(r - 1)) {
 >      for(i in 1:(r - m)) b.mtr[i,]<- 
(a.mtr[(i+1),]*(4^m)-a.mtr[i,])/(4^m-1)
 >      if(show & m!=(r-1) )  {
 >         cat("\n","Richarson improvement group No. ", m, "\n") 
        print(a.mtr[1:(r-m),], 12)
 >       }
 >    }
 > a.mtr[length(a.mtr)]})
 > }
 >
 > ## try it out
 >
 > richardson.grad(function(x){x3},2)
 >
 > #works fine... should return 12.
 >
 > # now try integrating something simple
 >
 > integrate(function(i) richardson.grad(function(x) x2,i),0,1)
 >
 > #also works fine, but instead try this:
 >
 > CDFLHP <-function(x,D,B)
 > pnorm((sqrt(1-B2)*qnorm(x)-D)/B)
 >
 >  integrate(function(i) richardson.grad(function(x) 
CDFLHP(x,-2,0.1),i),0,1)
 >
 > # fails, for some annoying reason, even tho richardson.grad is 
vectorised...
 >
 > XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
 > This is the error message:
 > Error in if (abs(a.mtr[(k - 1), i]) > 1e-20) a.mtr[k, i] <- 
(func(x1.vct) -  :
 >         missing value where TRUE/FALSE needed
 > In addition: Warning message:
 > NaNs produced in: qnorm(p, mean, sd, lower.tail, log.p)
 >
 >
 >
 >
 > Peter Dalgaard wrote:
 >
 >> "Gray Calhoun" <gcalhoun at ucsd.edu> writes:
 >>
 >>
 >>
 >>> Tolga,
 >>>
 >>> Look at numericDeriv.
 >>>
 >>>
 >>>> arbfun <- function(x) x2
 >>>> x <- 3
 >>>> numericDeriv(quote(arbfun(x)), "x")
 >>>>
 >>> [1] 9
 >>> attr(,"gradient")
 >>>     [,1]
 >>> [1,]    6
 >>>
 >> However, numericDeriv is not particularly intelligent. It is
 >> effectively doing what Tolga was trying not to. A more refined
 >> function could be a good idea, e.g. implementing higher order
 >> approximations, a tunable stepsize, box constraints...
 >>
 >>
 >>
 >>> --Gray
 >>>
 >>> On 3/12/06, Tolga Uzuner <tolga at coubros.com> wrote:
 >>>
 >>>> Hi,
 >>>>
 >>>> Suppose I have an arbitrary function:
 >>>>
 >>>> arbfun<-function(x) {...}
 >>>>
 >>>> Is there a robust implementation of a numerical derivative routine 
in R
 >>>> which I can use to take it's derivative ? Something a bit more than
 >>>> simple division by delta of the difference of evaluating the 
function at
 >>>> x and x+delta...
 >>>>
 >>>> Perhaps there is a way to do this using D or deriv but I could not
 >>>> figure it out. Trying:
 >>>>
 >>>> eval(deriv(function(x) arbfun(x),"x"),1)
 >>>>
 >>>> does not seem to work.
 >>>>
 >>>> Thanks,
 >>>> Tolga
 >>>>
 >>>> ______________________________________________
 >>>> R-help at stat.math.ethz.ch mailing list
 >>>> https://stat.ethz.ch/mailman/listinfo/r-help
 >>>> PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html
 >>>>
 >>>>
 >>> --
 >>> Gray Calhoun
 >>>
 >>> Economics Department
 >>> UC San Diego
 >>>
 >>> ______________________________________________
 >>> R-help at stat.math.ethz.ch mailing list
 >>> https://stat.ethz.ch/mailman/listinfo/r-help
 >>> PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html
 >>>
 >>>
 >>
 >>
 >
 > ______________________________________________
 > R-help at stat.math.ethz.ch mailing list
 > https://stat.ethz.ch/mailman/listinfo/r-help
 > PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html