Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to compute integral of partial derivative in R?

Tags:

r

integral

First time trying out the Deriv package. I was trying to compute a very simple integral of partial derivative as a start:

image

Here's my attempt:

junk <- function (m) {
  -m*eval(deriv(~((exp(0.1*(100-94.5)^2))/(exp(0.1*(100-94.5)^2)+exp(0.1* 
(m-95)^2)+exp(0.1*(m-96)^2))),"m"))
}
integrate(junk, lower = 100, upper = 100.5)

which gives me a wrong answer of -24.30757, rather than 12.383. Can anyone shed some light on how to fix this please? Thanks!

like image 641
wearpants Avatar asked Nov 25 '25 00:11

wearpants


1 Answers

The problem here is, that deriv() does not return an expression. For uni variate derivation use D()

junk <- function (m) {
  -m*eval(D(expression((exp(0.1*(100-94.5)^2))/(exp(0.1*(100-94.5)^2)+exp(0.1*(m-95)^2)+exp(0.1*(m-96)^2))),"m")
}
integrate(junk, lower = 100, upper = 100.5)

12.38271 with absolute error < 1.4e-13

For comparison:

return of deriv() - call

deriv(expression((exp(0.1*(100-94.5)^2))/(exp(0.1*(100-94.5)^2)+exp(0.1* 
                                                                        (m-95)^2)+exp(0.1*(m-96)^2))),"m"))
 expression({
     .expr4 <- exp(0.1 * (100 - 94.5)^2)
     .expr5 <- m - 95
     .expr8 <- exp(0.1 * .expr5^2)
     .expr10 <- m - 96
     .expr13 <- exp(0.1 * .expr10^2)
     .expr14 <- .expr4 + .expr8 + .expr13
     .value <- .expr4/.expr14
     .grad <- array(0, c(length(.value), 1L), list(NULL, c("m")))
     .grad[, "m"] <- -(.expr4 * (.expr8 * (0.1 * (2 * .expr5)) + 
         .expr13 * (0.1 * (2 * .expr10)))/.expr14^2)
     attr(.value, "gradient") <- .grad
     .value })

return of D() - expression

D(expression((exp(0.1*(100-94.5)^2))/(exp(0.1*(100-94.5)^2)+exp(0.1*(m-95)^2)+exp(0.1*(m-96)^2))),"m")
 -((exp(0.1 * (100 - 94.5)^2)) * (exp(0.1 * (m - 95)^2) * (0.1 * 
     (2 * (m - 95))) + exp(0.1 * (m - 96)^2) * (0.1 * (2 * (m - 
     96))))/(exp(0.1 * (100 - 94.5)^2) + exp(0.1 * (m - 95)^2) + 
     exp(0.1 * (m - 96)^2))^2)

so with your syntax and the eval() you need another expression which is provided by D(). Or you use the call object, either way is possible.

like image 144
mischva11 Avatar answered Nov 27 '25 14:11

mischva11



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!