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: mlev...@purdue.edu Phone: +1-765-496-7571 Fax: +1-765-494-0558 URL: www.stat.purdue.edu/~mlevins ________________________________ From: Ivan Krylov <ikry...@disroot.org> Sent: Thursday, June 13, 2024 6:21 AM To: Levine, Michael <mlev...@purdue.edu> Cc: r-help@r-project.org <r-help@r-project.org> Subject: Re: [R] Integration of functions with a vector argument [You don't often get email from ikry...@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" <mlev...@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]]
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
