[R] Numerical Derivatives in R
Paul Gilbert
pgilbert at bank-banque-canada.ca
Mon Mar 13 18:08:32 CET 2006
(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
More information about the R-help
mailing list