You could also use rSymPy to symbolically differentiate it. For example, based on the semi-automatic differentiation example on the rSymPy home page:
http://code.google.com/p/rsympy/#Semi-Automatic_Differentiation just replace the indented lines there with the indented lines here (I've also added # after each such line in case the email messes up the indentation.) The remaining lines are identical to that example. Note that x and bx hold symbolic variables so one must use [[ ... ]] rather than [ ... ] with them: > library(rSymPy) > # next file is sourced from devel version of rSymPy > source("http://rsympy.googlecode.com/svn/trunk/R/Sym.R";) > n <- 2 # > res <- x <- lapply(paste("x", 1:n, sep = ""), Sym) > sympy(paste("var(", shQuote(paste(x, collapse = " ")), ")")) [1] "(x1, x2)" > > b <- matrix(1:9, 3) # > bx <- list(NA, NA, NA) # > for(i in 1:3) bx[[i]] <- b[i, 1] + b[i, 2] * x[[1]] + b[i, 3] * x[[2]] # > result <- bx[[1]] / (bx[[1]] + bx[[2]] + bx[[3]]) # > > g <- sapply(x, deriv, expr = result, USE.NAMES = FALSE) > g2 <- sapply(g, sympy, USE.NAMES = FALSE) > g3 <- gsub("x([0-9]+)", "x[\\1]", g2) > gfun0 <- paste("grad <- function(x, n =", n, ") { c(", paste(g3, collapse = > ","), ") }") > gfun <- eval(parse(text = gfun0)) > gfun function(x, n = 2 ) { c( -15*(1 + 4*x[1] + 7*x[2])/(6 + 15*x[1] + 24*x[2])**2 + 4/(6 + 15*x[1] + 24*x[2]),-24*(1 + 4*x[1] + 7*x[2])/(6 + 15*x[1] + 24*x[2])**2 + 7/(6 + 15*x[1] + 24*x[2]) ) } etc. On Tue, May 12, 2009 at 2:25 AM, spencerg <spencer.gra...@prodsyse.com> wrote: > Hi, Paul: > > Your example is so complicated that I don't want to take the time to > check it. You apply "deriv" to an exponential divided by a sum of > exponentials, and I'm not convinced that your manual "Correct way" is > actually correct. It looks like you've followed the examples in the "deriv" > help page, so I would be more inclined to trust that, especially since it > matched the answer I got from genD, as follows. > > In your "genD" example, x01 and x02 should be x[1] and x[2]: > p1 <- function(x) {exp(b00.1+b01.1*x[1]+b02.1*x[2]) / > (exp(b00.1+b01.1*x[1]+b02.1*x[2])+ > exp(b00.2+b01.2*x[1]+b02.2*x[2])+ > exp(b00.3+b01.3*x[1]+b02.3*x[2])) - phat1} > test <- genD(p1, c(x01, x02)) > test$D > [,1] [,2] [,3] [,4] [,5] > [1,] -0.2012997 0.1296301 -0.03572875 0.07082898 -0.1261376 > > > The first two components of test$D here match your attr(eval(dp1.dx), > "gradient"). The next three are the lower triangular portion of the matrix > of second partials of the function "p1", per the "genD" documentation. > > The function numericGradient in the maxLik package could also be used > for this, I believe. However, I won't take the time here to test that. > > Hope this helps. Spencer Graves > > Paul Heinrich Dietrich wrote: >> >> Hi Spencer, >> Thanks for suggesting the genD function. In attempting it, I have >> rearranged my function from phat1 ~ ... to ... - 1, it apparently doesn't >> like the first one :) But when I run it, it tells me the partials are all >> zero. I'm trying out a simple MNL equation before I expand it to what I'm >> looking for. Here is what I tried (and I get different answers from a >> textbook solution, deriv(), and genD()): >> >> >>> >>> ### Variables for an observation >>> x01 <- rnorm(1,0,1) >>> x02 <- rnorm(1,0,1) >>> ### Parameters for an observation >>> b00.1 <- rnorm(1,0,1) >>> b00.2 <- rnorm(1,0,1) >>> b00.3 <- 0 >>> b01.1 <- rnorm(1,0,1) >>> b01.2 <- rnorm(1,0,1) >>> b01.3 <- 0 >>> b02.1 <- rnorm(1,0,1) >>> b02.2 <- rnorm(1,0,1) >>> b02.3 <- 0 >>> ### Predicted Probabilities for an observation >>> phat1 <- 0.6 >>> phat2 <- 0.3 >>> phat3 <- 0.1 >>> ### Correct way to calculate a partial derivative >>> partial.b01.1 <- phat1 * (b01.1 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b01.2 <- phat2 * (b01.2 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b01.3 <- phat3 * (b01.3 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b01.1; partial.b01.2; partial.b01.3 >>> >> >> [1] 0.04288663 >> [1] -0.1804876 >> [1] 0.1376010 >> >>> >>> partial.b02.1 <- phat1 * (b02.1 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b02.2 <- phat2 * (b02.2 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b02.3 <- phat3 * (b02.3 - (b01.1*phat1+b01.2*phat2+b01.3*phat3)) >>> partial.b02.1; partial.b02.2; partial.b02.3 >>> >> >> [1] 0.8633057 >> [1] 0.3171978 >> [1] 0.1376010 >> >>> >>> ### Derivatives for MNL >>> dp1.dx <- deriv(phat1 ~ exp(b00.1+b01.1*x01+b02.1*x02) / >>> >> >> + (exp(b00.1+b01.1*x01+b02.1*x02)+exp(b00.2+b01.2*x01+b02.2*x02)+ >> + exp(b00.3+b01.3*x01+b02.3*x02)), c("x01","x02")) >> >>> >>> dp2.dx <- deriv(phat2 ~ exp(b00.2+b01.2*x01+b02.2*x02) / >>> >> >> + (exp(b00.1+b01.1*x01+b02.1*x02)+exp(b00.2+b01.2*x01+b02.2*x02)+ >> + exp(b00.3+b01.3*x01+b02.3*x02)), c("x01","x02")) >> >>> >>> dp3.dx <- deriv(phat3 ~ exp(b00.3+b01.3*x01+b02.3*x02) / >>> >> >> + (exp(b00.1+b01.1*x01+b02.1*x02)+exp(b00.2+b01.2*x01+b02.2*x02)+ >> + exp(b00.3+b01.3*x01+b02.3*x02)), c("x01","x02")) >> >>> >>> attr(eval(dp1.dx), "gradient") >>> >> >> x01 x02 >> [1,] -0.01891354 0.058918 >> >>> >>> attr(eval(dp2.dx), "gradient") >>> >> >> x01 x02 >> [1,] -0.1509395 -0.06258685 >> >>> >>> attr(eval(dp3.dx), "gradient") >>> >> >> x01 x02 >> [1,] 0.169853 0.003668849 >> >>> >>> library(numDeriv) >>> dp1.dx <- function(x) {exp(b00.1+b01.1*x01+b02.1*x02) / >>> >> >> + (exp(b00.1+b01.1*x01+b02.1*x02)+exp(b00.2+b01.2*x01+b02.2*x02)+ >> + exp(b00.3+b01.3*x01+b02.3*x02)) - phat1} >> >>> >>> test <- genD(dp1.dx, c(phat1,b00.1,b01.1,b02.1,x01,x02)); test >>> >> >> $D >> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] >> [,14] >> [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] >> [,26] >> [1,] 0 0 0 0 0 0 0 0 0 0 0 >> 0 >> [,27] >> [1,] 0 >> >> $p >> [1] 6 >> >> $f0 >> [1] 0.05185856 >> >> $func >> function(x) {exp(b00.1+b01.1*x01+b02.1*x02) / >> (exp(b00.1+b01.1*x01+b02.1*x02)+exp(b00.2+b01.2*x01+b02.2*x02)+ >> exp(b00.3+b01.3*x01+b02.3*x02)) - phat1} >> >> $x >> [1] 0.600000000 1.401890082 -1.304531849 0.062833294 -0.143141379 >> [6] -0.005995477 >> >> $d >> [1] 1e-04 >> >> $method >> [1] "Richardson" >> >> $method.args >> $method.args$eps >> [1] 1e-04 >> >> $method.args$d >> [1] 1e-04 >> >> $method.args$zero.tol >> [1] 1.781029e-05 >> >> $method.args$r >> [1] 4 >> >> $method.args$v >> [1] 2 >> >> >> attr(,"class") >> [1] "Darray" >> >> >> >> >> >> spencerg wrote: >> >>> >>> Have you considered genD{numDeriv}? >>> If this does not answer your question, I suggest you try the >>> "RSiteSearch" package. The following will open a list of options in a web >>> browser, sorted by package most often found with your search term: >>> >>> library(RSiteSearch) >>> pd <- RSiteSearch.function('partial derivative') >>> pds <- RSiteSearch.function('partial derivatives') >>> attr(pd, 'hits') # 58 >>> attr(pds, 'hits')# 52 >>> summary(pd) >>> HTML(pd) >>> HTML(pds) >>> >>> The development version available via >>> 'install.packages("RSiteSearch", repos="http://R-Forge.R-project.org";)' also >>> supports the following: >>> pd. <- unionRSiteSearch(pd, pds) >>> attr(pd., 'hits')# 94 >>> HTML(pd.) >>> >>> >>> Hope this helps. Spencer Graves >>> >>> Paul Heinrich Dietrich wrote: >>> >>>> >>>> Quick question: >>>> >>>> Which function do you use to calculate partial derivatives from a model >>>> equation? >>>> >>>> I've looked at deriv(), but think it gives derivatives, not partial >>>> derivatives. Of course my equation isn't this simple, but as an >>>> example, >>>> I'm looking for something that let's you control whether it's a partial >>>> or >>>> not, such as: >>>> >>>> somefunction(y~a+bx, with respect to x, partial=TRUE) >>>> >>>> Is there anything like this in R? >>>> >>>> >>> >>> ______________________________________________ >>> 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. >>> >>> >>> >> >> > > ______________________________________________ > 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. > ______________________________________________ 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.
