How do I find the constant of a continuous joint probability distribution function in R?

Question

I'm wondering how to set up the calculation for a double integral to solve for the value of c for the problem below.

Consider the joint probability density function $f_{XY} (x, y) = c(x+y)$ over the range $0 < x < 3$ and $x < y < x + 2$.

I know how to do this for a joint discrete distribution, e.g.

Consider the joint probability mass function $f_{XY} (x, y) = c(x + y)$ over the nine points with $x = 1, 2, 3$ and $y = 1, 2, 3$.

For this, I'd create two variables so that a function $f_{XY}$ would return all possible values of $x + y$:

x = c(rep(1,3), rep(2,3), rep(3,3))
y = rep(c(1:3),3) 
f_XY = x+y
f_XY
c = 1/sum(f_XY)

I'm not sure how to create the variables for x and y when they're continuous, so I'd really appreciate any help!

Answer 1

The integral2 function of the pracma package is a possibility:

library(pracma)

f <- function(x ,y) x + y

integral2(f, xmin = 0, xmax = 3, ymin = function(x) x, ymax = function(x) x+2)

---

If you are not allowed to use a package, you can nest the integrate function:

inner <- function(x) integrate(function(y) x+y, lower = x, upper = x + 2)$value
integrate(Vectorize(inner), lower = 0, upper = 3)

Both methods give 24, an approximate value of the double integral.

---

Finally, you can use the SimplicialCubature package after noticing that the region of integration can be split into two triangles (= simplicies). Moreover, the integrand is a polynomial function, and then SimplicialCubature offers the possibility to get the exact value of the integral.

library(SimplicialCubature)

S1 <- cbind(c(0,0), c(3,3), c(0,2)) # first triangle
S2 <- cbind(c(0,2), c(3,3), c(3,5)) # second triangle
S <- array(c(S1, S2), dim = c(2, 3, 2))

P <- definePoly(coef = c(1,1), k = cbind(c(1,0), c(0,1)))
printPoly(P) # x + y

integrateSimplexPolynomial(P, S)
# 24