[R] Integration of functions with a vector argument

[Ivan Krylov]

В Wed, 12 Jun 2024 23:42:18 +0000
"Levine, Michael" <mlevins using purdue.edu> пишет:

> f.int1 <- function(x,y) {Vectorize (function(y) f(c(x,y),H=H,j=j))}

Vectorize returns a callable function(y), so wrapping it in a
function(x,y) will not work. Since you'd like to integrate over y, we
can perform the same transformation manually using vapply:

# A version of f(...) that can be integrated over y
f.int1 <- function(y, x, H, j)
 vapply(y, function(y) f(c(x,y), H = H, j = j), numeric(1))

The same transformation can be performed using Vectorize as follows:

f.int1 <- Vectorize(
 function(y, x, H, j) f(c(x,y), H = H, j = j),
 vectorize.args = 'y'
)

>     mrg.d1 <- function(x){
>     integrate(f.int1,lower=-3,upper=3)$value
>     }

We can then integrate f.int1 over y, making sure to pass the remaining
scalar parameters to the function:

# mrg.d1(x, H, j) = ∫ f(y,x,H,j) dy from -3 to 3
mrg.d1 <- function(x, H, j)
 integrate(f.int1, lower=-3, upper=3, x = x, H = H, j = j)$value

> mrg.cdf1 <- function(x){
>     integrate(mrg.d1,lower=-Inf,upper=x)$value
>     }

Unfortunately, mrg.d1 is once again not vectorised and may give wrong
answers or raise exceptions with non-scalar x. We can now perform the
same transformation so that the function would work with integrate():

# Version of mrg.d1 integratable over x
mrg.d1.int <- function(x, H, j)
 vapply(x, mrg.d1, numeric(1), H = H, j = j)

Here, instead of having to combine arguments like in f.int1, we can
just use vapply() to call the original mrg.d1 for every scalar element
inside the vector x given to mrg.d1.int. Alternatively,

mrg.d1.int <- Vectorize(mrg.d1, 'x')

A vectorised function can then be integrated:

# mrg.cdf1(x, H, j) = ∫ mrg.d1(x', H, j) dx' from -∞ to x
mrg.cdf1 <- function(x, H, j)
 integrate(mrg.d1.int, lower=-Inf, upper=x, H = H, j = j)$value

This double nested loop (produced by Vectorize() and mapply() or
manually with vapply()) may be not very fast. If you rewrite your
function f to accept matrices containing values of (x, y) for many
evaluations of the function and to return vectors of the resulting
values, you'll be able to use the CRAN package 'cubature'
<https://CRAN.R-project.org/package=cubature> to integrate it over
multiple variables at once.

-- 
Best regards,
Ivan