[R] Numerical Derivatives in R

[Tolga Uzuner]

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){x^3},2)

#works fine... should return 12.

# now try integrating something simple

integrate(function(i) richardson.grad(function(x) x^2,i),0,1)

#also works fine, but instead try this:

CDFLHP <-function(x,D,B)
pnorm((sqrt(1-B^2)*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) x^2
>>>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
>>
>>    
>>
>
>  
>