Integrating Functions in R

[omario]

I am using the R programming language. I am trying to evaluate the (multivariable) integral (through numerical integration) of the following function ("fitness") I defined :

fitness <- function(mu1, mu2, sigma1, sigma2, sigma12) {
    
    rho = sigma12/(sigma1*sigma2)
    a = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.4395244 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.4395244 - mu1)/sigma1) * ((3.7150650 - mu2)/sigma2)  + ((3.7150650 - mu2)/sigma2)^2
    b = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.7698225 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.7698225 - mu1)/sigma1) * ((2.4609162 - mu2)/sigma2)  + ((2.4609162 - mu2)/sigma2)^2
    c =  -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((2.5587083 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((2.5587083 - mu1)/sigma1) * ((0.7349388 - mu2)/sigma2)  + ((0.7349388 - mu2)/sigma2)^2
    d = -1/(2 * (1 - ((sigma12 / (rho))^2))) *  ((1.0705084 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.0705084 - mu1)/sigma1) * ((1.3131471 - mu2)/sigma2)  + ((1.3131471 - mu2)/sigma2)^2
    e = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((1.1292877 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.1292877 - mu1)/sigma1) * ((1.5543380 - mu2)/sigma2)  + ((1.5543380 - mu2)/sigma2)^2
    
    
    answer = answer = 1/ ((2*pi*sigma1*sigma2*sqrt((1 - rho^2)) / ((sigma1 * sigma2) * exp( a * b * c * d * e))))
    
}

I tried to integrate this function using the "integrate" function in R:

integrate(fitness, lower = c(-5, -5, -5, -5, -5), upper = c(5, 5, 5, 5, 5))

But I got the following error:

Error in integrate(fitness, lower = c(-5, -5, -5, -5, -5), upper = c(5,  : 
  length(lower) == 1 is not TRUE

Looking at the documentation for the "integrate" function (integrate function - RDocumentation), it seems like the "integrate" function is only designed to evaluate the integral single variable function.

What I tried so far: Looking at similar questions (e.g. numerical integration of multivariable functions), I saw that sometimes evaluating multivariable integrals requires you to "vectorize" your function:

Nd_vectorized_for <- function(s) {
  result <- numeric(length(s))
  for (i in 1:length(s)) {
    result[i] <- fitness(s[i])
    }
  result  ## or `return(result)`
  }

I tried again to evaluate the integral of this function:

Nd_vectorized_for( c(-5, -5, -5, -5, -5) :c(5, 5, 5, 5, 5) )

But it still returns an error:

Error in fitness(s[i]) : argument "sigma12" is missing, with no default

Can someone please show me how to fix this error?

Thanks!

[nirgrahamuk]

I'm afraid the stack overflow you found is a red herring for you as it is helping someone evaluate the integral at a vector of inputs, but the integration happening is only over 1-dimension. and not 5 like yours .

I think you need to use package curbature, but theres also a package called calculus that is a convenient wrapper for it.

example:

library(cubature)
library(calculus)

fitness <- function (a,b,c,d,e){
  a*5+4*b+c*3+d*2+e +1
}

integral(f=fitness,
         bounds = list(a=c(-1,1),
                       b=c(-1,1),
                       c=c(-1,1),
                       d=c(-1,1),
                       e=c(-1,1)))

I think however there are other issues with your function, as with division by zero potentially happening, and sqrt() of negative numbers, your function would often produce NaN, and integration result would likely be NaN.

Your function has an odd double assignment of answer , and I think its not returning a value. its good practice generally to test your function manually before attempting to use it in other functions, mistakes happen, good to catch them early.

[omario]

Thank you for your reply! I tried your code on another function:

fitness <- function(mu1, mu2, sigma1, sigma2, sigma12) {

rho = sigma12/(sigma1*sigma2)
a = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.4395244 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.4395244 - mu1)/sigma1) * ((3.7150650 - mu2)/sigma2)  + ((3.7150650 - mu2)/sigma2)^2
b = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.7698225 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.7698225 - mu1)/sigma1) * ((2.4609162 - mu2)/sigma2)  + ((2.4609162 - mu2)/sigma2)^2
c =  -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((2.5587083 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((2.5587083 - mu1)/sigma1) * ((0.7349388 - mu2)/sigma2)  + ((0.7349388 - mu2)/sigma2)^2
d = -1/(2 * (1 - ((sigma12 / (rho))^2))) *  ((1.0705084 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.0705084 - mu1)/sigma1) * ((1.3131471 - mu2)/sigma2)  + ((1.3131471 - mu2)/sigma2)^2
e = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((1.1292877 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.1292877 - mu1)/sigma1) * ((1.5543380 - mu2)/sigma2)  + ((1.5543380 - mu2)/sigma2)^2


answer = answer = 1/ ((2*pi*sigma1*sigma2*sqrt((1 - rho^2)) / ((sigma1 * sigma2) * exp( a * b * c * d * e))))

}

integral(f=fitness,
         bounds = list(mu1=c(-1,1),
                       mu2=c(-1,1),
                       sigma1=c(-1,1),
                       sigma2=c(-1,1),
                       sigma12=c(-1,1)))

This ran, but it gave a strange answer:

$value
[1] NaN

$error
[1] NaN

$cuba
$cuba$integral
[1] NaN

$cuba$error
[1] NaN

$cuba$neval
[1] 93

$cuba$returnCode
[1] 0

Does this look correct?

Thank you!

[nirgrahamuk]

I anticipated that, please reread my post.
Also test that your function gives results when you call it. I don't think you've done that.