My textbook, Introduction to Probability by Blitzstein and Hwang, presents the following example:
Example 7.2.2 (Expected distance between two Uniforms). Let $X$ and $Y$ be i.i.d. Unif$(0, 1)$ r.v.s. Find $E(|X - Y|)$.
Solution:
Since the joint PDF is 1 on the unit square $\{ (x, y) : x, y \in [0, 1] \}$, 2D LOTUS gives
$$\begin{align} E(|X - Y|) &= \int_0^1 \int_0^1 |x - y| \ dx dy \\ &= \int_0^1 \int_y^1 (x - y) \ dxdy + \int_0^1 \int_0^y (y - x) \ dxdy \\ &= 2 \int_0^1 \int_y^1 (x - y) \ dxdy = 1/3 \end{align}$$
First we broke up the integral into two parts so we could eliminate the absolute value; then we used symmetry.
The textbook defines the 2D LOTUS as follows:
Theorem 7.2.1 (2D LOTUS). Let $g$ be a function from $\mathbb{R}^2$ to $\mathbb{R}$. If $X$ and $Y$ are discrete, then
$$E(g(X, Y)) = \sum_x \sum_y g(x, y) P(X = x, Y = y).$$
If $X$ and $Y$ are continuous with joint PDF $f_{X, Y}$, then
$$E(g(X, Y)) = \int_{-\infty}^\infty \int_{-\infty}^\infty g(x, y) f_{X, Y}(x, y) \ dx dy.$$
LOTUS means Law of the Unconscious Statistician.
My multiple integral knowledge is rusty, so I would appreciate it if people could please take the time to explain what's going on for each step here:
$$\begin{align} E(|X - Y|) &= \int_0^1 \int_0^1 |x - y| \ dx dy \\ &= \int_0^1 \int_y^1 (x - y) \ dxdy + \int_0^1 \int_0^y (y - x) \ dxdy \\ &= 2 \int_0^1 \int_y^1 (x - y) \ dxdy = 1/3 \end{align}$$
Thank you.
EDIT:
The PDF of a continuous uniform random variable is $$f(x) = {\begin{cases}{\frac {1}{b-a}}&\mathrm {for} \ a\leq x\leq b,\\[8pt]0&\mathrm {for} \ x<a\ \mathrm {or} \ x>b\end{cases}}$$
Since $X$ and $Y$ are independent, my understanding is that the product of their individual PDFs is equal to the joint PDF:
$$p_{X, Y}(x, y) = p_X(x) p_Y(y) = (1)(1) = 1$$
By the LOTUS, we have $$E(\vert X - Y \vert) = \int_{0}^{1}\int_{0}^{1} \vert x - y \vert \, \mathrm{d}x \,\mathrm{d}y.$$ The $\vert x - y\vert$ term inside the integral takes one of the two forms $|x-y| = x-y$ for $x - y \geq 0$ and $|x -y| = y - x$ for $x - y < 0.$ So we split the region of interest into two regions (see attached image). 1
In the region defined by $x - y \geq 0$ (i.e. the region of the square below the line $y = x$), $y$ can take values between $0$ and $1$, and for each fixed $y$, $x$ runs from $y$ to $1$. Hence the integral of the function over this region can be written as $$\int_{0}^{1}\int_{y}^{1} (x-y)\,\mathrm{d}x\,\mathrm{d}y.$$ Similarly, for the region defined by $x - y < 0$, for each fixed $y$ between $0$ and $1$, $x$ runs from $0$ to $y$. So the integral of the function over this region is $$\int_{0}^{1}\int_{0}^{y} (y-x)\,\mathrm{d}x\,\mathrm{d}y.$$ This explains the step going from the first line to the second: $$\int_{0}^{1}\int_{0}^{1} \vert x - y \vert \, \mathrm{d}x \,\mathrm{d}y = \int_{0}^{1}\int_{y}^{1} (x-y)\,\mathrm{d}x\,\mathrm{d}y + \int_{0}^{1}\int_{0}^{y} (y-x)\,\mathrm{d}x\,\mathrm{d}y.$$ In fact, the region in the second integral can be specified in a different manner by noting that $x$ runs from $0$ to $1$ and for each fixed $x$, $y$ runs from $x$ to $1$. This allows us to reverse the order of integration to obtain $$\int_{0}^{1}\int_{0}^{y} (y-x)\,\mathrm{d}x\,\mathrm{d}y = \int_{0}^{1}\int_{x}^{1}(y-x)\,\mathrm{d}y\,\mathrm{d}x.$$ The value of the definite integral is independent of the dummy variables $x$ and $y$, so we may replace them as we please. In particular, we can switch their places to obtain $$\int_{0}^{1}\int_{0}^{y} (y-x)\,\mathrm{d}x\,\mathrm{d}y = \int_{0}^{1}\int_{y}^{1}(x-y)\,\mathrm{d}x\,\mathrm{d}y.$$ Note that this is an identical copy of the other integral in the second line, so this explains how the $2$ is obtained in the step from the second line to the third, and it remains to evaluate the result.