[R] Integration of functions with a vector argument

[Levine, Michael]

Dear Ivan,

Sorry for my slow response. Your advice works well; however, as you predicted, the time it takes to obtain any results is completely unacceptable (for example, a sequence of integrations we discussed is done once per iteration and even one iteration takes about 2 hours to complete).

I have one more question for you - hope you don't mind. I have heard of several packages used for numerical integration in R - cubature that you mentioned, mvQuad, and pracma.  My impression is that you think that Cubature is the best in your opinion. Is that so? If yes, do you know of any detailed discussion of this package beyond the two vignettes available on CRAN?

Thank you very much again for your help!

Yours sincerely,
Michael

Michael Levine
Associate Professor, Statistics

Department of Statistics
Purdue University
250 North University Street
West Lafayette, IN 47907 USA

email: mlevins using purdue.edu
Phone: +1-765-496-7571
Fax:   +1-765-494-0558
URL:   www.stat.purdue.edu/~mlevins
________________________________
From: Ivan Krylov <ikrylov using disroot.org>
Sent: Thursday, June 13, 2024 6:21 AM
To: Levine, Michael <mlevins using purdue.edu>
Cc: r-help using r-project.org <r-help using r-project.org>
Subject: Re: [R] Integration of functions with a vector argument

[You don't often get email from ikrylov using disroot.org. Learn why this is important at https://aka.ms/LearnAboutSenderIdentification ]

---- External Email: Use caution with attachments, links, or sharing data ----


�� 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://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fcran.r-project.org%2Fpackage%3Dcubature&data=05%7C02%7Cmlevins%40purdue.edu%7C9753b609c1d34745311908dc8b9296b9%7C4130bd397c53419cb1e58758d6d63f21%7C0%7C0%7C638538708933465815%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C0%7C%7C%7C&sdata=uCwVWSrsLemT4j1nT9UiEt%2BQxH99Z8erBorg0LEYmWI%3D&reserved=0<https://cran.r-project.org/package=cubature>> to integrate it over
multiple variables at once.

--
Best regards,
Ivan

	[[alternative HTML version deleted]]