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: [EMAIL PROTECTED] [mailto:r-help- > [EMAIL PROTECTED] On Behalf Of Paul Gilbert > Sent: Monday, March 13, 2006 12:09 PM > To: r-help@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" <[EMAIL PROTECTED]> 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 <[EMAIL PROTECTED]> 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@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@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@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@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@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