[R] Partial Derivatives in R
spencerg
spencer.graves at prodsyse.com
Tue May 12 08:25:30 CEST 2009
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.
>>
>>
>>
>
>
More information about the R-help
mailing list