# [R] Numerical Derivatives in R

Mon Mar 13 21:24:56 CET 2006

```Hi Paul, Tolga, and others:

I had also written some codes to compute derivatives, jacobians, and
Hessians.  Please see the attached file for the code.  I will be happy to
help out with the development of a package and/or with the documentation
process.

Best,
Ravi.

> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch [mailto:r-help-
> bounces at stat.math.ethz.ch] On Behalf Of Paul Gilbert
> Sent: Monday, March 13, 2006 12:09 PM
> To: r-help at stat.math.ethz.ch
> Subject: Re: [R] Numerical Derivatives in R
>
> (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
>  >
>  >
>  > #works fine... should return 12.
>  >
>  > # now try integrating something simple
>  >
>  >
>  > #also works fine, but instead try this:
>  >
>  > CDFLHP <-function(x,D,B)
>  > pnorm((sqrt(1-B2)*qnorm(x)-D)/B)
>  >
> 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
>  >>>     [,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
> 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
> 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