[R] Numerical Derivatives in R
Ravi Varadhan
rvaradhan at jhmi.edu
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
> >
> > 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
>
> ______________________________________________
> 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
More information about the R-help
mailing list