Dear Petr Savicky, this helped indeed. Thank you very much.
Gildas ----- Mail original ----- > De: "Petr Savicky" <savi...@cs.cas.cz> > À: r-help@r-project.org > Envoyé: Vendredi 23 Mars 2012 09:39:37 > Objet: Re: [R] Computing High Order Derivatives (Numerically) > On Fri, Mar 23, 2012 at 12:35:57AM +0100, Gildas Mazo wrote: > > Dear R users, > > > > Let f be a function over d variables x1,..,xd. I want to compute the > > k^th-order derivative with respect to x1,..,xk (k<=d). I have a by > > hand solution (see below) using an iterating code using D. However, > > I expect d to be high and f to be complicated. Then I want a vector > > x to be the input, instead of x1,..,xd. How to avoid the x1 <- x[1]; > > x2 <- x[2], etc steps in the code below? Moreover, D uses symbolic > > differentation and then eval evaluates the output to get a numerical > > result. But is there a way to compute the desired derivatives > > numerically directly (without using symbolic calculus at all)? > > Finally, what is the most efficient and fast way to get a numerical > > result for such derivatives? > > > > Thank you very much in advance, > > Gildas > > > > ### Code ### > > ### dif takes a function f, an order k, and a vector x as input. f > > must be a function of x1,..,xd with d >= k. The correspondance is > > done between xi and x[i]. The expression for f must be at the last > > row of the body function. > > dif <- function(f,k,x){ > > o <- list() > > n <- length(body(f)) > > o[[1]] <- body(f)[[n]] > > for (i in 1:k){ > > xi <- paste("x",i,sep="") > > o[[i+1]] <- D(o[[i]],name=xi) > > } > > x1 <- x[1] > > x2 <- x[2] > > x3 <- x[3] > > eval(o[[k+1]]) > > } > > > > ### Examples ### > > ## function to differentiate > > f <- function(x){ > > x1 <- x[1] > > x2 <- x[2] > > x3 <- x[3] > > 0.5*x1*x2*x3^2 > > } > > ## derivative w.r.t. x1, x2 and x3 at the point (1,2,3). > > dif(f,3,c(1,2,3)) > > > > ### My Questions ### > > ## how to avoid to write by hand xi <- x[i] ?? > > ## is there a way in R to compute such derivatives without using > > symbolic calculation but numerical compuation instead. > > Hi. > > For the first question, try the following > > dif <- function(f,k,x){ > o <- list() > n <- length(body(f)) > o[[1]] <- body(f)[[n]] > for (i in 1:k){ > xi <- paste("x",i,sep="") > o[[i+1]] <- D(o[[i]],name=xi) > assign(xi, x[i]) > } > eval(o[[k+1]]) > } > > For the second question, try the following. > > x <- c(1, 2, 3) > k <- length(x) > grid <- as.matrix(expand.grid(rep(list(c(0, 1)), times=k))) > signs <- 1 - 2*(rowSums(1 - grid) %% 2) > for (eps in 2^-(5:20)) { > xeps <- eps*grid + rep(x, each=nrow(grid)) > print(sum(signs*apply(xeps, 1, FUN=f))/eps^k) > } > > [1] 3.015625 > [1] 3.007812 > [1] 3.003906 > [1] 3.001953 > [1] 3.000977 > [1] 3.000488 > [1] 3.000244 > [1] 3.000122 > [1] 3 > [1] 3 > [1] 3 > [1] 3 > [1] 4 > [1] 0 > [1] 0 > [1] 0 > > If the above is computed in an exact arithmetic, then > with "eps" converging to zero, the result converges to > the required derivative. Since the numerical computations > are done with a rounding error, too small "eps" yields > a completely wrong result. The choice of a good "eps" > depends on the function and on "k". For a high "k", there > may even be no good "eps". See the considerations at > > http://en.wikipedia.org/wiki/Numerical_derivative > > where the choice of "eps" is discussed in the simplest > case of a univariate function. > > Hope this helps. > > Petr Savicky. > > ______________________________________________ > R-help@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. -- Gildas Mazo PhD student MISTIS team at INRIA Grenoble, France ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.