What's the probability that three points determine an acute triangle?

Three distinct points are chosen at random from the unit square. The goal is to find the probability that they form an acute triangle. I started working on this because I want to know how to approach a problem of this sort, where the sample space seems to be something like $[0,1]^2$. However I've had no experience dealing with pdf's on $[0,1]^2$, $\mathbb{R}^2$, etc.

So far, given any two points, I can find the probability that the third point will work via a simple geometry/area argument (see below). But then what to do with that? I suspect integration, but of what and over what domain?

Suppose $(a,b)$ and $(c,d)$ are two points in the square. Construct a line containing $(a,b)$ and a line containing $(c,d)$ so that both are perpendicular to the segment joining the points. The probability that the third point will work is the area of the region between the two lines and inside the square, but outside the circle whose diameter is the segment joining the points. (If the third point is inside the circle then an obtuse triangle is formed, and if the third point is on the circle a right triangle is formed.)

Solution 1:

After working this out and writing it down, I found that it was already done in 1969: Eric Langford, The probability that a random triangle is obtuse, Biometrika (1969) 56 (3): 689-690.

For this kind of problem, it's important to make good use of the symmetries and choose convenient coordinates. A significant simplification is achieved by noting that a triangle has at most one obtuse angle; thus the three events of one of the angles being obtuse are mutually exclusive, and thus the probability of the triangle being obtuse is just three times the probability $p$ of any given one of the angles being obtuse. Thus you don't have to deal with both the circle and the strip between the perpendiculars, since the probability of the third point being in the circle is $p$ and the probability of the third point being in the strip is $1-2p$ and you can obtain the desired probability $1-3p$ from either of them. In fact, you don't even need the strip, since the probability of the third point being beyond one of the perpendiculars is also $p$.

So there are three different approaches; we can either integrate either the area beyond one of the perpendiculars or the area of the strip or the area of the circle over all positions of two points. The advantage of using the perpendiculars or the strip is that you have to deal only with straight lines and not with circles; the disadvantage is that the integration region is complicated due to the corners. The advantage of the second approach is that you avoid the corners and get a relatively simple integration region; the disadvantage is that you get somewhat complicated integrands for the circle segments that extend beyond the square.

I'll work out the circle approach here. The easy part is integrating over the area of the circle without regard to whether it's inside the square. This is

$$ \def\lint{\int\limits} \begin{eqnarray} && \lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\lint_0^1\mathrm dy_1\lint_0^1\mathrm dy_2\,\pi\left(\left(\frac{x_2-x_1}2\right)^2+\left(\frac{y_2-y_1}2\right)^2\right) \\ &=& 2\cdot\frac\pi4\lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\,(x_2-x_1)^2 \\ &=& \frac\pi2\left(\frac13-2\cdot\frac12\cdot\frac12+\frac13\right) \\ &=& \frac\pi{12}\;. \end{eqnarray} $$

From this we have to subtract the integral over the circle segments that extend beyond the square. By symmetry, this is four times the integral over the segments that extend beyond the bottom edge of the square. These segments depend on $x_1$ and $x_2$ only through $x:=|x_2-x_1|$, so we can get rid of one integral by the transformation

$$ \lint_0^1\mathrm dx_1\lint_0^1\mathrm dx_2\to2\lint_0^1\mathrm dx_1\lint_{x_1}^1\mathrm dx_2\to2\lint_0^1\mathrm dx\lint_x^1\mathrm dx_2\to2\lint_0^1\mathrm dx\,(1-x)\;. $$

It will also be convenient to transform from $y_1$ and $y_2$ to $s:=y_1+y_2$ and $a:=y_2-y_1$. We can restrict to $a\gt0$ by symmetry, and the symmetry factor of $2$ cancels the Jacobi determinant $2$ of the transformation:

$$\lint_0^1\mathrm dy_1\lint_0^1\mathrm dy_2\to\lint_0^1\mathrm da\lint_a^{2-a}\mathrm ds\;.$$

The centre of the circle is at $y=s/2$, and its squared radius is $(x/2)^2+(a/2)^2$, so the angle of the segment extending beyond the lower edge is

$$\alpha=2\arccos\frac s{\sqrt{x^2+a^2}}\;.$$

The area of a circular segment with angle $\alpha$ and radius $r$ is $(\alpha-\sin\alpha)r^2/2$. Thus, with $\sin2\arccos u=2u\sqrt{1-u^2}$, the area is

$$A(x,a,s)=\frac14\left(\left(x^2+a^2\right)\arccos\frac s{\sqrt{x^2+a^2}}-s\sqrt{x^2+a^2-s^2}\right)\;.$$

The factor $1/4$ cancels with the symmetry factor $4$ for the four edges of the square. We need to integrate this over all values for which there is a circle segment extending beyond the lower edge, that is, for $s^2\lt x^2+a^2$. We also have $s\lt2-a$. These two bounds coincide for $x^2+a^2=(2-a)^2$, that is $a=1-x^2/4$, so we need two integrals for $a$ above and below that value:

$$\lint_0^1\mathrm dx\,(1-x)\lint_0^{1-x^2/4}\mathrm da\lint_a^{\sqrt{x^2+a^2}}\mathrm ds\,A(x,a,s)+\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\lint_a^{2-a}\mathrm ds\,A(x,a,s)\;.$$

The integrals over the second term in $A(x,a,s)$ are relatively easy to evaluate:

$$\lint_0^1\mathrm dx\,(1-x)\lint_0^{1-x^2/4}\mathrm da\lint_a^{\sqrt{x^2+a^2}}\mathrm ds\,s\sqrt{x^2+a^2-s^2}=\frac{37}{2520}$$

(computation) and

$$\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\lint_a^{2-a}\mathrm ds\,s\sqrt{x^2+a^2-s^2}=\frac1{840}$$

(computation), and the sum is $\displaystyle\frac{37}{2520}+\frac1{840}=\frac1{63}$.

The first term is a bit harder to integrate; Wolfram|Alpha fails to do it in one go. If $F(s)$ is an antiderivative for this term with respect to $s$ (where I'm suppressing $a$ and $x$ to avoid notational clutter), then the integral we need is

$$\lint_0^1\mathrm dx\,(1-x)\left(\lint_0^{1-x^2/4}\mathrm da\,F\left(\sqrt{x^2+a^2}\right)+\lint_{1-x^2/4}^1\mathrm da\,F(2-a)-\lint_0^1\mathrm da\,F(a)\right)\;.$$

With

$$F(s)=\left(x^2+a^2\right)\left(s\arccos\frac s{\sqrt{x^2+a^2}}-\sqrt{x^2+a^2-s^2}\right)$$

(computation) we have

$$F\left(\sqrt{x^2+a^2}\right)=0\;,$$

$$F(2-a)=\left(x^2+a^2\right)\left((2-a)\arccos\frac{2-a}{\sqrt{x^2+a^2}}-\sqrt{x^2+4a-4}\right)\;,$$

$$F(a)=\left(x^2+a^2\right)\left(a\arccos\frac a{\sqrt{x^2+a^2}}-x\right)\;.$$

Wolfram|Alpha evaluates the integral over $F(a)$ in one go if you give it extra time:

$$\lint_0^1\mathrm dx\,(1-x)\lint_0^1\mathrm da\,F(a)=\frac{-52+21\pi-48\log2}{720}$$

(computation). The hardest part is the integral over $F(2-a)$. The antiderivative with respect to $a$ looks rather daunting at first, until you realize that it's zero at the lower limit $a=1-x^2/4$. Thus we only need it for $a=1$. The result is still complicated but can be simplified to

$$\frac1{12}\Bigg(\left(18x^2+5\right)\arccos\frac1{\sqrt{x^2+1}}-\frac1{15}x\left(6x^4+245x^2+75\right) \\ +x^3\left(8\log\left(x^2+1\right)-3x\arctan x\right)\Bigg).$$

A final integration then yields

$$\lint_0^1\mathrm dx\,(1-x)\lint_{1-x^2/4}^1\mathrm da\,F(2-a)=-\frac{517}{3150}+\frac{\pi-\log 2}{15}$$

(computation).

Putting it all together, we obtain

$$ \begin{eqnarray} p &=& \frac\pi{12}-2\left(-\frac{517}{3150}+\frac{\pi-\log 2}{15}-\frac{-52+21\pi-48\log2}{720}-\frac1{63}\right) \\ &=& \frac{97}{450}+\frac\pi{120}\;. \end{eqnarray} $$

(computation). Thus, the probability $1-3p$ of the triangle being acute is

$$\boxed{\displaystyle\frac{53}{150}-\frac\pi{40}\approx0.2748}$$

in agreement with numerical simulations.

While we're at it, we can also calculate the probability for points uniformly distributed on the boundary of the square. We can assume the first point to be on one side, say the left one, and average the probability of the angle at the third point being obtuse as the second and third points range over all sides.

If the second point is on the same side as the first, the angle is obtuse only if the third point is also on the same side and between the other two; the probability for that is $1/3$ by symmetry.

If the second point is on the opposite side from the first, the angle is obtuse if the third point is also on one of these sides and between the other two, for a contribution of $2\cdot1/3$. It can also be obtuse if the third point is on one of the remaining sides. The width $w$ of the interval on the bottom edge for which this is the case is determined by the circle equation:

$$ \left(\frac12\right)^2+\left(\frac{y_1-y_2}2\right)^2=\left(\frac w2\right)^2+\left(\frac{y_1+y_2}2\right)^2\;,\\ w^2=1-4y_1y_2 \;. $$

The corresponding integrals are

$$\lint_0^{1/4}\mathrm dy_1\lint_0^1\mathrm dy_2\,\sqrt{1-4y_1y_2}=\frac49-\frac{\log2}3$$

(computation) and

$$\lint_{1/4}^1\mathrm dy_1\lint_0^{1/(4y_1)}\mathrm dy_2\,\sqrt{1-4y_1y_2}=\frac{\log2}3$$

(computation) for a total contribution of $8/9$ (multiplied by $2$ for the two remaining edges).

If the second point is on a side adjacent to the left one, say, the bottom one, the angle is obtuse if the third point is on one of these two sides and between the other points, which gives two contributions of $1/2$. The angle can also be obtuse if the third point is on one of the remaining sides, say, the right one. Again, the width of the interval for which this is the case is determined by the circle equation:

$$ \left(\frac w2\right)^2+\left(1-\frac x2\right)^2=\left(\frac x2\right)^2+\left(\frac y2\right)^2\;,\\ w^2=y^2+4x-4 \;. $$

The corresponding integral is

$$\lint_0^1\mathrm dy\lint_{1-y^2/4}^1\mathrm dx\,\sqrt{y^2+4x-4}=\frac1{24}$$

(computation).

Putting it all together, we have

$$ \begin{eqnarray} p &=& \frac14\cdot\frac14\left(\frac13+\frac23+\frac89+2\left(2\cdot\frac12+2\cdot\frac1{24}\right)\right) \\ &=& \frac{73}{288}\;, \end{eqnarray} $$

and the probability $1-3p$ of the triangle being acute is

$$ \frac{23}{96}\approx0.2396 $$

in agreement with numerical simulations.

Here's a table with the probabilities of three points forming an acute triangle for various distributions of the points in the plane:

$$ \begin{array}{c|c|c} \text{circle}&\frac14&0.25\\ \hline \text{disk}&\frac4{\pi^2}-\frac18&0.2803\\ \hline \partial(\text{square})&\frac{23}{96}&0.2396\\ \hline \text{square}&\frac{53}{150}-\frac\pi{40}&0.2748\\ \hline \mathbb R^6\text{ symmetry}&\frac14&0.25 \end{array} $$

References for these values are at MathWorld under Disk Triangle Picking, Square Triangle Picking and Gaussian Triangle Picking and in an article by Stephen Portnoy, A Lewis Carroll Pillow Problem: Probability of an Obtuse Triangle, Statistical Science, Volume 9, Number 2 (1994), 279-284 that contains a good discussion of what it might mean to choose triangles randomly and ends with the brilliant sentence "I would suggest that the value and enjoyment of mathematics is greatest when it is used to provide imperfect and incomplete solutions to ill-formulated problems".

The entry "$\mathbb R^6\text{ symmetry}$" refers to a part of Portnoy's article where he shows that any distribution invariant under rotations in $\mathbb R^6$ (treating all point coordinates as a single vector in $\mathbb R^6$) must lead to the same result and derives that result for a Gaussian distribution. Here's a simple argument for the result $1/4$ that I haven't found anywhere: Parametrize the coordinates $x_1,y_1,x_2,y_2$ of the first two points with $\phi,a_1,a_2,h$, where $\phi$ is the angle through which they have to be rotated to let the vector from one to the other point along the positive $x$ direction, $a_1$ and $a_2$ are the $x$ coordinates of the rotated points and $h$ is their common $y$ coordinate. The Jacobian of this transformation is $|a_1-a_2|$. Thus, to find the probability that one of the angles at the first two points is obtuse, we can average over all values of $\phi$, $a_1$, $a_2$, $h$ and the rotated coordinates $a_3,h_3$ of the third point, weighting each configuration by the distance between the first two points. Since $\phi$, $h$ and $h_3$ are independent of $a_1$, $a_2$ and $a_3$ and the obtuseness of the angles doesn't depend on them, we only have to average over $a_1$, $a_2$ and $a_3$. By symmetry, all six permutations of $a_1$, $a_2$ and $a_3$ are equally likely. They form three pairs of mirror images. The weights of the two pairs that have one of the first two points in the middle add up to to the weight of the pair that has the third point in the middle. Thus the probability for one of two angles to be obtuse is $1/2$, and hence the probability for one of three angles to be obtuse is $3/4$.

Solution 2:

While I wouldn't want to quibble with @joriki's excellent answer, I can't help but think that there is a simpler solution. Remembering a bit from probability, here is an indication, as yet still incomplete, of how one might proceed.

I would start by formulating the problem as that of being given six i.i.d. uniform random variates as three pairs $\{(X_i,Y_i)\}_{i=1}^3\in[0,1]^{2\times3}$ considered as vertices $V_i$ of a triangle, each with angle $\Theta_i$ (allowing indices to "wrap" or "cycle" modulo $3$), and being asked for the probability that each $\Delta_{ij}=\operatorname{sign}(\cos\Theta_i)$ is positive. However, from the inner product formula $$\cos\Theta_i =(x_{i+1}-x_i)(x_{i-1}-x_i) +(y_{i+1}-y_i)(y_{i-1}-y_i)$$ we can see that the problem can be reformulated using intermediate variables $U=(x_1-x_3)(x_2-x_3)$ and $V$ independent but identically distributed to $U$, by noting that $$\cos\Theta_i =U+V\,.$$ Before computing $Pr(U+V>0)$, let us compute $Pr(U>u)$, from which the distribution of $U+V$ will follow. To compute this for $u\in[-1,1]$, let us call the variates $(x,y,z)$ rather than $(x_1,x_2,x_3)$ and start, for example, with $u=0$, for which $$ Pr(U>0)=\int_0^1 \left(z^2+(1-z)^2\right)dz=\int_0^1(2w+1)\,dz=\frac23 $$ This can be seen by noting that, for $z\in(0,1)$ fixed (and $w=z^2-z$), the product $(x-z)(y-z)$ is positive inside the squares going from $(0,0)$ to $(z,z)$ and from $(z,z)$ to $(1,1)$ (the grey & green regions below), and negative on the rest of $(x,y)\in[0,1]^2$ (in white & red). For the general integral $Pr(U > u)$ with $u\in[-1,1]$, instead of the bounding lines $x,y=z$, we have the bounding curve $(x-z)(y-z)=u$, which on $[0,1]^2$ has up to two components. Here is a graph showing what the curves (and regions) look like when $z=\frac25$ for $U>.08$ (in green), $U<-.08$ (red), and for $U\in(0,.08)$ (grey).

$$ Pr(U>u) =\int_{[0,1]^3} \chi\left((x-z)(y-z) > u\right) \,dx\,dy\,dz =\int_0^1 A(u,z)\,dz$$ for the area function $$ A(u,z) = \left\{\matrix{ 0 &\quad& u &\ge& v &\quad& u\ge\tfrac14~\vee~ u\in\left(0,\tfrac14\right) \wedge z\in\left[0,\sqrt u\right]\cup\left[1-\sqrt u,1\right] \\\\ A^+(u,z) &\quad& u &\in& \left(0,~ v\right) &\quad& u\in\left(0,\tfrac14\right) \wedge z\in\left(\sqrt u,\,1-\sqrt u\right) \\\\ 2w+1 &\quad& u & = & 0 &\quad& u=0 \wedge z\in\left[0,\,1\right] \\\\ A^-(u,z) &\quad& u &\in& \left(w,~ 0\right) &\quad& u\in\left(-\tfrac14,\,0\right) \wedge z\in\left(\tfrac12-\sqrt{u+\tfrac14},\,\tfrac12+\sqrt{u+\tfrac14}\right) \\\\ 0 &\quad& u &\le& w &\quad& u\le-\tfrac14 \vee \left|z-\frac12\right|>\sqrt{u+\frac14} }\right.$$ where $v=\min(z,~1-z)^2\in\left[0,~\frac14\right]$, $w=z(z-1)\in\left[-\frac14,~0\right]$, and $A^\pm(u,z)$ have yet to be determined but are fairly clear from the sketch above. Using symmetry about $y=x$ and, as reference values, the total area $1$ and $Pr(U>0)=\frac23$ as above (which has squares and no curves for all $z$), we can write $$ \eqalign{ A^+(u,z) &= z^2 + (1-z)^2 - 2 \int_{-z}^{-\min(z,\sqrt{u})}\left( \frac{u}{x}-x \right)\,dx - 2 \int_{\min(1-z,\sqrt{u})}^{1-z}\left( x-\frac{u}{x} \right)\,dx \\ &= z^2 + (1-z)^2 + \left[x^2-2u\ln|x|\right]_{-z,~1-z}^{-\min(z,\sqrt{u}),~\min(1-z,\sqrt{u})} \\ &= \min(z^2,u) - u\ln\frac{\min(z^2,u)}{z^2} + \min\left((1-z)^2,u\right) - u\ln\frac{\min\left((1-z)^2,u\right)}{(1-z)^2} \\ A^-(u,z) &= 1 -\int_{-z}^{-\min\left(z,\,\frac{z-1}{u}\right)}\left(1-\frac{u}{x}\right)\,dx -\int_{\min\left(1-z,~-\frac{z}{u}\right)}^{1-z}\left(1+\frac{u}{x}\right)\,dx \\ &= 1 - \Bigl[x\Bigr] _{ -z,~ \min\left(1-z,~-\frac{z}{u}\right)} ^{1-z,~-\min\left(z,\,\frac{z-1}{u}\right)} + u \Bigl[\ln|x|\Bigr] _{-z,~1-z} ^{-\min\left(z,\,\frac{z-1}{u}\right) ,~\min\left(1-z,~-\frac{z}{u}\right)} \\ &= 2z + \min\left( z,\frac{ z-1}{u }\right) - \min\left(1-z,\frac{-z }{u }\right) + u\ln\min\left(1 ,\frac{ z-1}{u z }\right) + u\ln\min\left(1 ,\frac{-z }{u(1-z)}\right) } $$ where I have used multiple endpoints in the bracketed antiderivative differences to save space. This is indeed so far quite similar to Langford's 1969 solution, the biggest differences being that he defines the CDF, $Pr(U\le u)$, and expresses the boundaries directly in terms of $u$ rather than implicitly in terms of $z$ (and $v$ & $v$), which I have now also done in the rightmost column above. From this, we can see that we need only concern ourselves with the two middle cases, $u\in\left(0,\tfrac14\right)$ and $u\in\left(-\tfrac14,0\right)$. In the former, $$ \eqalign{ Pr(U>u) &=\int_{\sqrt u}^{1-\sqrt u}A^+(u,z)\,dz\\ %&=\int_{\sqrt u}^{1-\sqrt u}\left[ % \min(z^2,u) % - u\ln\frac{\min(z^2,u)}{z^2} % + \min\left((1-z)^2,u\right) % - u\ln\frac{\min\left((1-z)^2,u\right)}{(1-z)^2} %\right]dz\\ &=\int_{\sqrt u}^{1-\sqrt u}\left[ u - u\ln\frac{u}{z^2} + u - u\ln\frac{u}{(1-z)^2} \right]dz\\ &=2u\int_{\sqrt u}^{1-\sqrt u}\Bigl[ 1-\ln u + \ln z + \ln(1-z) \Bigr]dz\\ &=2u\Bigl(1-\ln u\Bigr)\Bigl[z\Bigr]_{\sqrt u}^{1-\sqrt u} +\int_{\sqrt u}^{1-\sqrt u}\Bigl[\ln z-\ln z\Bigr]dz\\ &=2u\Bigl(1-\ln u\Bigr)\Bigl(1-2\sqrt u\Bigr) } $$ since $\int_a^b\ln(1-z)dz=-\int_{1-a}^{1-b}\ln z\,dz$ using the substitution $t=1-z$. In the latter case, setting $a=\tfrac12-\sqrt{u+\tfrac14}$, $b=\tfrac12+\sqrt{u+\tfrac14}$ where $\tfrac14 < u < 0 < z < 1$, note that $(z-\tfrac12)^2 < u+\tfrac14$ $\implies$ $z(z-1)\le u$ $\implies$ $z\ge\frac{u}{z-1}$, $\frac{z-1}{u}\ge z$ and $1-z\ge\frac{-u}{z}$ so that $$ \eqalign{ Pr(U>u) &=\int_a^b A^-(u,z)\,dz\\ &=\int_a^b \left[2z+z-\frac{-z}{u}+u\ln1+u\ln\frac{-z}{u(1-z)}\right]dz\\ &=\int_a^b \left[ \left(3+\tfrac1u\right)z +u\Bigl(\ln z-\ln(1-z)-\ln|u|\Bigr)\right]dz\\ &=\int_a^b \left[\left(3+\tfrac1u-u\ln|u|\right)z + 2u\ln z\right]dz\\ &= \tfrac12\left(3+\tfrac1u-u\ln|u|\right)\left[z^2\right]_a^b + 2u\Bigl[z\ln z-z\Bigr]_a^b\\ &= \left(3+\tfrac1u-u\ln|u|\right)\sqrt{u+\tfrac14} + 2u\Bigl[z\ln z \Bigr]_a^b - 4u \sqrt{u+\tfrac14}\\ } $$ ...to be continued.