[R] Partial Derivatives in R
Gabor Grothendieck
ggrothendieck at gmail.com
Tue May 12 13:23:12 CEST 2009
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.graves at 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 at 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 at 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.
>
More information about the R-help
mailing list