2005/9/26, Dimitris Rizopoulos <[EMAIL PROTECTED]>: > AFAIK the deriv() is for symbolic derivatives, so I don't know if it > will work in your case.
numericDeriv surely computes numerical gradient. The problem here is its interface, which requires the vector of symbols involved in the expression to differenziate. > > Best, > Dimitris > > > ----- Original Message ----- > From: "Antonio, Fabio Di Narzo" <[EMAIL PROTECTED]> > To: "Dimitris Rizopoulos" <[EMAIL PROTECTED]> > Cc: <[email protected]> > Sent: Monday, September 26, 2005 10:58 AM > Subject: Re: getting variable length numerical gradient > > > 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: <[email protected]> > > 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. > > > > > > ______________________________________________ > > > [email protected] 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 > > > > > > > Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm > > ______________________________________________ [email protected] mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
