The PDF of the area of a random triangle in a square

Three points are chosen uniformly and independently at random on the unit rectangle $(0,1)\times(0,1)$.

What is the probability distribution (in the form of a cumulative distribution function and/or a probability density function supported on $[0,1/2]$) of the area of the triangle whose vertices are so chosen? What are the moments of this distribution?

The area might be expressed as the norm of the sum of three three-dimensional vector cross products, as

$$\frac12\|p_1\times p_2+p_2\times p_3+p_3\times p_1\|.$$

Closely related: Probability of lies a point in a random triangle. See also http://mathworld.wolfram.com/SquareTrianglePicking.html and http://www.diva-portal.org/smash/get/diva2:644463/FULLTEXT01.pdf

Let $X$ be the collection of random variables taking values from $[0,1]$.
Let $U$ be the sub-collection of uniform random variables.

For each $f \in X$, let $\rho_f \in [0,1] \to \mathbb{R}\cup \{\infty\}$ be its PDF whenever it is defined.
In particular, if $f \in U$, then $\rho_f(t) = 1$ for $t \in [0,1]$.

Let $(x_1,y_1)$, $(x_2,y_2)$, $(x_3,y_3)$ be 3 random points from $U^2$. Let $T$ be the triangle formed from them and $A$ be its area. Introduce a bunch of variables $x_\ell, x_m, x_u, y_\ell, y_m, y_u$ such that $$\begin{cases} x_\ell \le x_m \le x_u,& \{ x_\ell, x_m, x_u \} = \{ x_1, x_2, x_3 \}\\ y_\ell \le y_m \le y_u,& \{ y_\ell, y_m, y_u \} = \{ y_1, y_2, y_3 \}\\ \end{cases}$$

Let $w = x_u - x_\ell$ and $h = y_u - y_\ell$, they are the width and height of the bounding box of $T$.
It is easy to check $w, h \in X$ and

$$\rho_w(t) = \rho_h(t) = 6t(1-t)\quad\text{ for } t \in [0,1]$$

From this, we can deduce $wh \in X$ with

$$\begin{align} \rho_{wh}(\ell) &= \int_0^1 \int_0^1 \rho_w(s)\rho_h(t) \delta(st - \ell) ds dt = \int_\ell^1 \rho_w(s)\rho_h\left(\frac{\ell}{s}\right) \frac{ds}{s}\\ &= 36\ell \left[ -(\ell+1)\log \ell + 2(\ell-1) \right] \end{align} $$

Let $\displaystyle\;u = \frac1w ( x_m - x_\ell)$, $\displaystyle\;v = \frac1h ( y_m - y_\ell )$ and $\eta$ be the ratio $\displaystyle\;\frac{2A}{wh}$.

Aside from events of probability zero ($w = 0$ or $h = 0$), $u, v$ and $\eta$ are well defined. Furthermore, $u, v$ behave as if they belongs to $U$.

Given any realization of $( x_\ell, y_\ell)$, $( x_m, y_m )$, $( x_u, y_u )$, there are $(3!)^2 = 36$ ways to assign the coordinates to the 3 vertices $(x_1, y_1)$, $(x_2,y_2)$ and $(x_3,y_3)$.

• Case I - In $12$ of the ways, two of the vertices occupy an pair of opposite corners of the bounding box. One example is $$ \begin{cases} (x_1,y_1) &= (x_\ell,y_\ell),\\ (x_2,y_2) &= (x_m,y_m)\\ (x_3,y_3) &= (x_u, y_u) \end{cases} \quad\implies\quad \eta = \left|\begin{matrix} 0 & 0 & 1\\ u & v & 1\\ 1 & 1 & 1 \end{matrix}\right| = |u-v|$$ In these cases, $\eta$ behaves like a random variable in $X$ with $$\rho_{\eta,I}(t) = 2(1-t)$$

• Case II - For the remaining $24$ ways, only one vertices is occupying an corner of the bounding box and the remaining two vertices lying on edges of the bounding box. One example is $$ \begin{cases} (x_1,y_1) &= (x_\ell,y_\ell),\\ (x_2,y_2) &= (x_u, y_m),\\ (x_3,y_3) &= (x_m, y_u) \end{cases} \quad\implies\quad \eta = \left|\begin{matrix} 0 & 0 & 1\\ 1 & v & 1\\ u & 1 & 1 \end{matrix}\right| = 1-uv$$ In these cases, $\eta$ behaves like a random variable $\in X$ with $$\rho_{\eta,II}(t) = \rho_{uv}(1-t) = \int_{1-t}^1 \frac{ds}{s} = -\log(1-t)$$

Combine these two scenarios, we find aside from events of probability zero, we can treat $\eta$ as a random variable in $X$ with

$$\rho(t) = \frac13 \rho_{\eta,I}(t) + \frac23 \rho_{\eta,II}(t) = \frac23\left[(1-t) - \log(1-t)\right]$$

This implies we can treat $2A = \eta wh$ as a random variable $\in X$ with PDF

$$\rho_{2A}(t) = \int_{t}^1 \rho_{wh}(s)\rho_\eta\left(\frac{t}{s}\right)\frac{ds}{s}$$ With help of a CAS, this evaluates to an ugly mess: $$6(1-t)(1 - (5t+1)\log(1-t)) + 6t^2\log(t)(\log(t) - 5) + 12t(t+2)\left(\mathrm{Li}_2(t) - \frac{\pi^2}{6}\right) $$ where $$\mathrm{Li}_2(z) = - \int_0^z \frac{\log(1-t)}{t} dt = -\int_0^1 \frac{\log(1-zt)}{t} dt$$ is the Dilogarithm function.

As a double check, one can use above expression of $\rho_{2A}(t)$ to compute the expected area of $T$: $$\mathbb{E}[A] = \frac12 \int_0^1 \rho_{2A}(t) t dt$$ This gives us an expected area $\frac{11}{144}$ matching what is in other answer.