From: Paul Gilbert <pgilbert_at_bank-banque-canada.ca>

Date: Tue 14 Mar 2006 - 04:08:32 EST

* > 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")
*

* > }
*

> 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@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@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@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 Received on Tue Mar 14 04:40:15 2006

Date: Tue 14 Mar 2006 - 04:08:32 EST

(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
**

in 1:n) {

> x1.vct <- x2.vct <- x

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.> }

print(a.mtr, 12)

> for(m in 1:(r - 1)) {

(a.mtr[(i+1),]*(4^m)-a.mtr[i,])/(4^m-1)

> if(show & m!=(r-1) ) {

print(a.mtr[1:(r-m),], 12)

> }

> a.mtr[length(a.mtr)]})

CDFLHP(x,-2,0.1),i),0,1)

> # fails, for some annoying reason, even tho richardson.grad is

vectorised...

(func(x1.vct) - :

> missing value where TRUE/FALSE needed

in R

>>>> which I can use to take it's derivative ? Something a bit more than

function at

>>>> x and x+delta...

http://www.R-project.org/posting-guide.html

>>> Gray Calhoun

http://www.R-project.org/posting-guide.html

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 Received on Tue Mar 14 04:40:15 2006

*
This archive was generated by hypermail 2.1.8
: Tue 14 Mar 2006 - 09:09:13 EST
*