(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