[R] Numerical Derivatives in R
Spencer Graves
spencer.graves at pdf.com
Tue Mar 14 00:37:16 CET 2006
Hi, Paul, Tolga, Ravi, et al.:
What about using splines to compute both derivatives and integrals?
These could be used, for example, to summarize complicated, expensive
Monte Carlo results in a relatively simple form could then be
manipulated analytically and quickly converted to other forms by a
variety of helper functions.
For example, I'm currently struggling with a multilevel modeling
application for which "lmer" won't help, because I can't transform it to
an additive normal situation. I need to integrate out the random
parameters to obtain the marginal likelihood, which I then maximize over
variations in hyperparameters. This problem could be simplified by
first building a multidimensional spline model over variations in the
variable of integration and the hyperparameters. The resulting spline
model could then be used for parameter estimation, computing confidence
intervals, etc. I'd like to extend this further for modeling
applications where known special cases and asymptotics might be blended
and extended using spline model adjustments or interpolants.
What do you think?
Best Wishes,
Spencer Graves
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
More information about the R-help
mailing list