[R] getting variable length numerical gradient

[Dimitris Rizopoulos]

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


Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm