# [R] Numerical Derivatives in R

Mon Mar 13 18:08:32 CET 2006

```(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
> 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")
>>>>
>>>  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