Tnx very much Dimitris, your code does what I need. I've just adapted it to my needs (e.g., I don't deal with scalar functions), and so solved my problem.
Given this, is there a way to use the deriv function in the base package, within this context (variable length vector of indipendent variables)? Best, Antonio, Fabio Di Narzo. On 9/25/05, Dimitris Rizopoulos <[EMAIL PROTECTED]> wrote: > maybe you can find the following function useful (any comments are > greatly appreciated): > > fd <- function(x, f, scalar = TRUE, ..., eps = > sqrt(.Machine$double.neg.eps)){ > f <- match.fun(f) > out <- if(scalar){ > if(length(f0 <- f(x, ...)) != length(x)) > stop("'f' must be vectorized") > x. <- x + eps * pmax(abs(x), 1) > c(f(x., ...) - f0) / (x. - x) > } else{ > n <- length(x) > res <- array(0, c(n, n)) > f0 <- f(x, ...) > ex <- pmax(abs(x), 1) > for(i in 1:n){ > x. <- x > x.[i] <- x[i] + eps * ex[i] > res[, i] <- c(f(x., ...) - f0) / (x.[i] - x[i]) > } > res > } > out > } > > > ## Examples > > x <- seq(-3.3, 3.3, 0.1) > all.equal(fd(x, pnorm, mean = 0.5), dnorm(x, mean = 0.5)) > > > # Approximate the Hessian matrix for a logistic regression > > # the score vector function > gn <- function(b, y, X){ > p <- as.vector(plogis(X %*% b)) > -colSums(X * (y - p)) > } > > # We simulate some data and fit the logistic regression > n <- 800 > x1 <- runif(n,-3, 3); x2 <- runif(n, -3, 3) > pr <- plogis(0.8 + 0.4 * x1 - 0.3 * x2) > y <- rbinom(n, 1, pr) > fm <- glm(y ~ x1 + x2, binomial) > > ## The Hessian using forward difference approximation > fd(fm$coef, gn, scalar = FALSE, y = y, X = cbind(1, x1, x2)) > > ## The true Hessian > solve(summary(fm)$cov.unscaled) > > > I hope it helps. > > Best, > Dimitris > > ---- > Dimitris Rizopoulos > Ph.D. Student > Biostatistical Centre > School of Public Health > Catholic University of Leuven > > Address: Kapucijnenvoer 35, Leuven, Belgium > Tel: +32/(0)16/336899 > Fax: +32/(0)16/337015 > Web: http://www.med.kuleuven.be/biostat/ > http://www.student.kuleuven.be/~m0390867/dimitris.htm > > > ----- Original Message ----- > From: "Antonio, Fabio Di Narzo" <[EMAIL PROTECTED]> > To: <R-help@stat.math.ethz.ch> > Sent: Sunday, September 25, 2005 11:37 AM > Subject: [R] getting variable length numerical gradient > > > > Hi all. > > I have a numerical function f(x), with x being a vector of generic > > size (say k=4), and I wanna take the numerically computed gradient, > > using deriv or numericDeriv (or something else). > > > > My difficulties here are that in deriv and numericDeric the function > > is passed as an expression, and one have to pass the list of > > variables > > involved as a char vector... So, it's a pure R programming question. > > > > > > Have a nice sunday, > > Antonio, Fabio Di Narzo. > > > > ______________________________________________ > > 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 > > > > > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm > > ______________________________________________ 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