We can visualize part of the following joint density function using as a 3D surface.
\[ f(x, y) = 6e^{-2x- 3y} \quad \quad (x, y > 0) \]
f <- function(x, y) {
6 * exp(-2 * x - 3 * y)
}
dx <- .1
dy <- .1
x <- seq(0, 3, by = dx)
y <- seq(0, 3, by = dy)
z <- outer(x, y, f)
persp(x, y, z, theta = 25, phi = 25)
Note that the function is defined or all \(x, y > 0\), but most of the density is in this region around the origin.
We represent the following probability graphically:
\[ P(X > 1, Y > .5) = \int_{y = .5}^{\infty} \int_{x = 1}^{\infty} f(x, y) dx, dy \]
z2 <- z
z2[1:10,] <- 0
z2[, 1:5] <- 0
persp(x, y, z2, theta = 25, phi = 25, zlim = c(0, max(z)))
On board.
We approximate the integral by discretizing the continuous volume into a series of rectangles of height \(f(x, y)\) and base area \(dxdy\) and summing the probabilities.
\[ P(X > 1, Y > .5) \approx \sum_{x,y\textrm{ in grid}} f(x, y) dx, dy \]
sum(z2 * dx * dy)
## [1] 0.03797101
We can visualize \(P(X > 1 \textrm{ or } Y > .5)\),
z3 <- z
z3[1:10, 1:5] <- 0
persp(x, y, z3, theta = 25, phi = 25, zlim = c(0, max(z)))
then approximate it computationally.
sum(z3 * dx * dy)
## [1] 0.4165207
Or we can visualize the complement, \(P(X < 1 \textrm{ or } Y < .5)\),
z4 <- z
z4[11:31, ] <- 0
z4[, 6:31] <- 0
persp(x, y, z4, theta = 25, phi = 25, zlim = c(0, max(z)))
and compute 1 minus its volume.
1 - sum(z4 * dx * dy)
## [1] 0.1421356
The fact that these methods aren’t in close agreement is likely due to the disagreement between the discrete and continuous densities around the origin.
We can visualize the process of integrating out \(x\) to get \(f_Y(y)\).
fy <- colSums(z)
plot(y, fy, pch = 16)
lines(y, fy, col = "tomato", lwd = 2)
And equivalently for \(f_X(x)\)
fx <- rowSums(z)
plot(y, fy, pch = 16)
lines(y, fy, col = "tomato", lwd = 2)
points(x, fx, pch = 16)
lines(x, fx, col = "steelblue", lwd = 2)