[R] Numerical Derivatives in R
Martin Maechler
maechler at stat.math.ethz.ch
Tue Mar 14 11:15:18 CET 2006
>>>>> "Spencer" == Spencer Graves <spencer.graves at pdf.com>
>>>>> on Mon, 13 Mar 2006 15:37:16 -0800 writes:
Spencer> Hi, Paul, Tolga, Ravi, et al.: What about using
Spencer> splines to compute both derivatives and integrals?
That seems good and pretty "obvious" to me, and for that reason there
have been the functions D1ss() and D2ss() in package 'sfsmisc'
"forever" (since 1992, longer than R exists)
install.package("sfsmisc")
library(sfsmisc)
?D1ss
## yes, even the help file says that I wrote these in 1992
They are not as nice as you suggest below,
(since at that time I knew a bit less than now :-)
but might still be a useful start.
Martin Maechler, ETH Zurich
Spencer> These could be used, for example, to summarize
Spencer> complicated, expensive Monte Carlo results in a
Spencer> relatively simple form could then be manipulated
Spencer> analytically and quickly converted to other forms
Spencer> by a variety of helper functions.
Spencer> For example, I'm currently struggling with a
Spencer> multilevel modeling application for which "lmer"
Spencer> won't help, because I can't transform it to an
Spencer> additive normal situation. I need to integrate out
Spencer> the random parameters to obtain the marginal
Spencer> likelihood, which I then maximize over variations
Spencer> in hyperparameters. This problem could be
Spencer> simplified by first building a multidimensional
Spencer> spline model over variations in the variable of
Spencer> integration and the hyperparameters. The resulting
Spencer> spline model could then be used for parameter
Spencer> estimation, computing confidence intervals, etc.
Spencer> I'd like to extend this further for modeling
Spencer> applications where known special cases and
Spencer> asymptotics might be blended and extended using
Spencer> spline model adjustments or interpolants.
Spencer> What do you think? Best Wishes, Spencer
Spencer> Graves
Spencer> Ravi Varadhan wrote:
>> 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
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> ______________________________________________
>>> 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
Spencer> ______________________________________________
Spencer> R-help at stat.math.ethz.ch mailing list
Spencer> https://stat.ethz.ch/mailman/listinfo/r-help PLEASE
Spencer> do read the posting guide!
Spencer> http://www.R-project.org/posting-guide.html
More information about the R-help
mailing list