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