Search code examples
rintegral

How to compute integral of partial derivative in R?

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!

Solution

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