[R] Partial Derivatives in R

[Paul Heinrich Dietrich]

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

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