[R] Numerical Derivatives in R

[Gabor Grothendieck]

Note that warning message.  Its trying to evaluate qnorm outside of
its allowable domain.

On 3/12/06, Tolga Uzuner <tolga at coubros.com> 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){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
> >>
> >>
> >>
> >
> >
> >
>
> ______________________________________________
> 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
>