[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:

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
>    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"))
>>
>>>
>>>
>>
>>             x01      x02
>> [1,] -0.01891354 0.058918
>>
>>>
>>>
>>
>>            x01         x02
>> [1,] -0.1509395 -0.06258685
>>
>>>
>>>
>>
>>          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
>>> 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