# [R] Partial Derivatives in R

Paul Heinrich Dietrich paul.heinrich.dietrich at gmail.com
Tue May 12 12:48:10 CEST 2009

```Thanks for showing me how to use genD correctly and that it matches deriv
here.  Like you say, there must be a problem with the manual way of doing
it, and I will look at it closely.  Thanks again.

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

--
View this message in context: http://www.nabble.com/Partial-Derivatives-in-R-tp23470413p23500226.html
Sent from the R help mailing list archive at Nabble.com.

```

More information about the R-help mailing list