В 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://CRAN.R-project.org/package=cubature> to integrate it over
multiple variables at once.
--
Best regards,
Ivan
______________________________________________
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.
