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.
