When $(X, Y)$ is not uniform, probability is still found by integrating, now with a density over a region of the plane. The limits come from sketching the region, exactly as in double integrals over general regions. This generalizes geometric probability and the joint pmf.
Key idea
A joint density satisfies $f(x, y) \ge 0$ and $\iint f\,dA = 1$, and for a region $R$, $$P((X, Y) \in R) = \iint_{R} f(x, y)\,dA.$$ Integrate only where the event holds and $f > 0$.
Why it works
As in one dimension, $f(x, y)$ is probability per unit area: a small rectangle near $(x, y)$ with sides $\Delta x$ and $\Delta y$ has probability about $f(x, y)\,\Delta x\,\Delta y$. Adding over a region gives a double Riemann sum, whose limit is the double integral. Outside the support the density is 0, so those parts of the region add nothing. The work is in the limits: sketch the support and the event, then for each $x$ read off the range of $y$ that lies in both.
Worked examples
Worked example: The constant
$f(x, y) = c(x + y)$ on the unit square. Find $c$.
- $$\int_{0}^{1}\int_{0}^{1} (x + y)\,dy\,dx = \int_{0}^{1} \left(x + \frac{1}{2}\right)dx = 1$$
- So $c = 1$.
Answer: $c = 1$
Worked example: A triangle region
With $f(x, y) = x + y$ on the unit square, find $P(X + Y < 1)$.
- For each $x$ in $[0, 1]$, $y$ runs from 0 to $1 - x$.
- Inner: $$\int_{0}^{1 - x} (x + y)\,dy = x(1 - x) + \frac{(1 - x)^{2}}{2} = \frac{1 - x^{2}}{2}$$
- Outer: $\int_{0}^{1} \frac{1 - x^{2}}{2}\,dx = \frac{1}{2} \cdot \frac{2}{3} = \frac{1}{3}$.
Answer: $\frac{1}{3}$
Worked example: A density on a triangle
$f(x, y) = 2$ for $0 \le x \le y \le 1$. Find $P(Y > 2X)$.
- The event $y > 2x$ inside the triangle: $x$ from 0 to $\frac{1}{2}$, and $y$ from $2x$ to 1 (then $y > x$ automatically).
- $$\int_{0}^{1/2}\int_{2x}^{1} 2\,dy\,dx = \int_{0}^{1/2} 2(1 - 2x)\,dx = 1 - \frac{1}{2} = \frac{1}{2}$$
Answer: $\frac{1}{2}$
Supports can be unbounded. If $f(x, y) = \lambda\mu e^{-\lambda x - \mu y}$ for $x, y > 0$ (two independent exponential times), then $$P(X < Y) = \int_{0}^{\infty}\int_{x}^{\infty} \lambda\mu e^{-\lambda x - \mu y}\,dy\,dx = \int_{0}^{\infty} \lambda e^{-(\lambda + \mu)x}\,dx = \frac{\lambda}{\lambda + \mu}:$$ the clock with the larger rate tends to ring first.
Common mistake
Integrating over the whole support instead of the part where the event happens, or using limits that run outside the support, where $f = 0$. Check: sketch both regions; the limits must describe their overlap, and the answer must lie between 0 and 1.