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