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" <antonio.fabio at gmail.com>
To: <R-help at 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.