Sylvester’s four-point problem in a triangle

Solution. Step 1: the disjoint-events dichotomy. Four points in the plane in general position either lie in convex position (their convex hull is a quadrilateral) or else exactly one of the four lies inside the triangle formed by the other three. Let $E_i$ denote the event “$P_i$ lies strictly inside $\triangle P_j P_k P_l$” (the triangle of the other three). The events $E_1, E_2, E_3, E_4$ are mutually exclusive, and $$\mathbb{P}(\text{convex}) + \mathbb{P}(E_1 \cup E_2 \cup E_3 \cup E_4) = 1.$$

Step 2: reduce to the expected-area lemma. By symmetry, $\mathbb{P}(E_i)$ is the same for each $i$. For $i = 4$, condition on $P_1, P_2, P_3$: $$\mathbb{P}(E_4) = \mathbb{E}\!\left[\frac{\operatorname{Area}(\triangle P_1 P_2 P_3)}{\operatorname{Area}(T)}\right] = \frac{\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)]}{\operatorname{Area}(T)}.$$

Step 3: the expected area is $\operatorname{Area}(T)/12$. We claim that for three uniform random points in a triangle $T$, $$\mathbb{E}[\operatorname{Area}(\triangle P_1 P_2 P_3)] = \frac{\operatorname{Area}(T)}{12}. \tag{$\star$}$$ By affine invariance of both sides under linear maps $T \to T$, we may take $T$ to be the standard triangle with vertices $(0,0), (1,0), (0,1)$, of area $\tfrac12$. Write $P_i = (X_i, Y_i)$. The signed area is $\tfrac12 D$ with $$D = (X_2 - X_1)(Y_3 - Y_1) - (X_3 - X_1)(Y_2 - Y_1),$$ so $\operatorname{Area}(\triangle P_1 P_2 P_3) = \tfrac12 |D|$. Direct integration (using the density $2$ on $T$ and repeated Fubini) gives $\mathbb{E}[|D|] = \tfrac1{12}$, hence $\mathbb{E}[\operatorname{Area}] = \tfrac1{24} = \tfrac1{12} \cdot \operatorname{Area}(T)$. The computation is tedious but elementary.

Step 4: assemble. From Steps 1–3, $$\boxed{\mathbb{P}(\text{convex}) = 1 - 4 \cdot \tfrac{1}{12} = 1 - \tfrac13 = \tfrac23.} \qedhere$$

import numpy as np
# 4 uniform points in standard triangle (fold-sampled); directly test convex
# position by checking no point is strictly inside the triangle of the other 3.
N = 10**6
u, v = np.random.rand(2, 4, N)
mask = u + v > 1
u, v = np.where(mask, 1-u, u), np.where(mask, 1-v, v)
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
    return (w1 > 0) & (w2 > 0) & (1 - w1 - w2 > 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: {2/3:.5f}")
# sim: 0.66736   exact: 0.66667