Sylvester’s four-point problem in a square, Geometric Probability

Problem 45.1

Let $P_1, P_2, P_3, P_4$ be four points chosen independently and uniformly at random in the unit square. What is the probability that they form the vertices of a convex quadrilateral?

Solution. Exactly as in Chapter 43, the disjoint-events dichotomy gives $\mathbb{P}(\text{convex}) = 1 - 4 \cdot \mathbb{P}(P_4 \in \triangle P_1 P_2 P_3),$ and the probability on the right equals $\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)] / \operatorname{Area}(\text{square})$ .

By the expected-area computation in Chapter 44, $\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)] = \frac{11}{144},$ and the unit square has area $1$ . Therefore $\mathbb{P}(\text{convex}) = 1 - 4 \cdot \tfrac{11}{144} = 1 - \tfrac{44}{144} = \tfrac{100}{144},$ which reduces to $\boxed{\mathbb{P}(\text{convex}) = \frac{25}{36} \approx 0.694.} \qedhere$

Two points of comparison, showing how the convex-position probability depends on the shape of the region $K$ :

$K$	$\mathbb{P}(\text{convex})$	numerically
triangle	$\tfrac23$	$0.6667\ldots$
square	$\tfrac{25}{36}$	$0.6944\ldots$
disc	$1 - \tfrac{35}{12 \pi^2}$	$0.7045\ldots$

The monotonic increase reflects intuitively that a “rounder” shape makes it harder for one of the four random points to be trapped inside the triangle of the other three. Sylvester’s original 1864 problem was to find the shape minimising $\mathbb{P}(\text{convex})$ ; Blaschke proved in 1917 that the triangle is the extremal case (and so the value $\tfrac23$ of Chapter 43 is the minimum over all convex regions).

import numpy as np
N = 10**6
u, v = np.random.rand(2, 4, N)
def inside_tri(pu, pv, au, av, bu, bv, cu, cv):
    det = (bv-cv)*(au-cu) + (cu-bu)*(av-cv)
    w1 = ((bv-cv)*(pu-cu) + (cu-bu)*(pv-cv)) / det
    w2 = ((cv-av)*(pu-cu) + (au-cu)*(pv-cv)) / det
    w3 = 1 - w1 - w2
    return (w1 > 0) & (w2 > 0) & (w3 > 0)
any_in = np.zeros(N, dtype=bool)
for i in range(4):
    j, k, l = [m for m in range(4) if m != i]
    any_in |= inside_tri(u[i], v[i], u[j], v[j], u[k], v[k], u[l], v[l])
print(f"sim: {(1 - any_in).mean():.5f}   exact: {25/36:.5f}")
# sim: 0.69465   exact: 0.69444