Solve this probability problem symbolically

Consider a unit square, Pick two points P and Q uniformly at random inside the square, What is the probability that |PQ|>1?

I tried solve this problem

Integrate[Boole[(x1-x2)^2+(y1-y2)^2>1],{x1,0,1},{y1,0,1},{x2,0,1},{y2,0,1}]

above code with NIntegrate given 0.025074 , but I want a symbolically result as except is $\frac{19}{6}-\pi$ .

• To show that the numeric result is equivalent to the analytic result: (NIntegrate[Boole[(x1 - x2)^2 + (y1 - y2)^2 > 1], {x1, 0, 1}, {y1, 0, 1}, {x2, 0, 1}, {y2, 0, 1}] + Pi // RootApproximant) - Pi evaluates to 19/6 - Pi Mar 23 '16 at 22:26
• @mathe thank you for accepting my answer. However, I think the other answers are more instructive and efficient, e.g. wolfies Mar 28 '16 at 22:45

This is not quick (includes J.M. comment):

pdf = UniformDistribution;
td = TransformedDistribution[(x - y)^2, {x, y} \[Distributed] pdf];
zd = TransformedDistribution[
a + b, {a \[Distributed] td, b \[Distributed] td}];

then

ans = 1 - FullSimplify[CDF[zd,1]]

yields the desired result. • Could also try Probability[ EuclideanDistance[{x1, y1}, {x2, y2}] < 1, {{x1, y1} \[Distributed] pdf, {x2, y2} \[Distributed] pdf}] but it's taking a fair while too... (also I think z should be zd in your CDF) Mar 23 '16 at 11:40
• You could use pdf = UniformDistribution for short. Mar 23 '16 at 13:49
• @J.M. thank you (as always for feedback). The other answers are much more instructive and I have commented to OP the same. I felt guilty so: ubpdqnmathematica.wordpress.com/2016/03/29/… Mar 29 '16 at 4:05

This is an interesting problem, because the difficulty is not the concept, but rather how to compute it (efficiently). Given points $$(X_i,Y_i)$$ distributed Uniformly on the unit square, we are interested in $$P\big[\sqrt{(X_2-X_1)^2 + (Y_2-Y_1)^2} \; > \; 1\big]$$

Let $$X = X_2 - X_1$$ denote the difference of two standard Uniform random variables, which is well-known to be a symmetric Triangular distribution on (-1,1). Similarly, let $$Y = Y_2-Y_1$$. By independence, the joint pdf of $$(X,Y)$$, say $$f(x,y)$$ is then:

f = (1-Abs[x]) (1-Abs[y]);       domain[f] = {{x,-1,1}, {y,-1,1}};

We seek:

All done. This takes just 2 seconds to evaluate, starting from a fresh kernel, where I am using the Prob function from the mathStatica package for Mathematica (how I roll, being one of the authors), but one can equally use in-built Mathematica functions to the same effect:

dist = ProbabilityDistribution[(1 - Abs[x]) (1 - Abs[y]), {x, -1, 1}, {y, -1, 1}];
Probability[Sqrt[x^2 + y^2] > 1, Distributed[{x, y}, dist]]
• Well. two seconds and, say, US$127 Mar 23 '16 at 18:49 • How much did you pay for your copy of Mathematica? Mar 23 '16 at 18:50 • As I noted above, one can follow the same approach with built-in functions. Here is some code to do that: dist = ProbabilityDistribution[(1 - Abs[x]) (1 - Abs[y]), {x, -1, 1}, {y, -1, 1}]; Probability[Sqrt[x^2 + y^2] > 1, Distributed[{x, y}, dist]] Mar 23 '16 at 19:08 • @woflies Nice! Do you mind adding it to your answer so I can upvote it? Mar 23 '16 at 19:10 • Sure. The point was never meant to be about actual code - but the conceptual approach, which works nicely either way here. Mar 23 '16 at 19:39 What's with all the heavy lifting and machinations? d = ProductDistribution[{TriangularDistribution[{-1, 1}], 2}]; Probability[a^2 + b^2 > 1, {a, b} \[Distributed] d]$\frac{19}{6}-\pi\$

Finishes in a few seconds on a netbook...

• After some testing of my own, that makes two netbooks. :) Mar 24 '16 at 0:13
• Looks like Wolfram could improve it code by replacing the difference of two uniforms by a triangular. The computation is much faster.
– A.G.
Mar 28 '16 at 23:30
• @ciao when I looked on my phone I didn't see who posted this...but I am not surprised...I expiated my sins :ubpdqnmathematica.wordpress.com/2016/03/29/… , and suggested OP look at all the other nicer answers...:) Mar 29 '16 at 4:09
• @ubpdqn - LOL - don't post here as much, busy herding a startup. Nice blog - I actually ran across it in the past (came up in a G-Search result), had it bookmarked as a cool place. Did not know it was you until now... though had I paid attention to the address it would have been obvious.
– ciao
Mar 30 '16 at 8:22
• @ciao best wishes for start up...blog is meaningless musings...play. Have fun :) Mar 30 '16 at 8:25

This is a "please don't reinvent the wheel" type of answer. From Michael Trott, here

p[l_] := Piecewise[{
{2 l (l^2 - 4 l + π), 0 <= l <= 1},
{2 l (4 Sqrt[l^2 - 1] - (l^2 + 2 - π) - 4 ArcTan[Sqrt[l^2 - 1]]), 1 < l <= Sqrt}}]

With[{n = 1}, Integrate[ p[l], {l, 1, Sqrt@2}]] // Simplify

(* 19/6 - π *)

The problem can also be solved in a generalized form. What is the probability that the distance between the two points is greater than t ?

The PDF for the distance of two random variables equally distributed between 0 and 1 is

f[u_] := 1 - u Sign[u]

Hence the probability in question is

p[t_] = Integrate[f[u] f[v] Boole[u^2 + v^2 > t^2], {v, -1, 1}, {u, -1, 1},
Assumptions -> t > 0] // Simplify

(* Because Latex has unfavorable line breaks I write down the Mathematica expressions for the two regions of t *)

Simplify[p[t], 0 < t < 1]

(* Out= 1 - π t^2 + (8 t^3)/3 - t^4/2 *)

Simplify[p[t], 1 < t < Sqrt]

(* Out= 1/12 (8 + 24 t^2 - 3 π t^2 + 6 t^4 -
16 Sqrt[-1 + t^2] - 32 t^2 Sqrt[-1 + t^2] - 18 t^2 ArcCsc[t] +
30 t^2 ArcTan[Sqrt[-1 + t^2]]) *)

(* Here's Latex *)

$$\begin{array}{cc} -\frac{t^4}{2}+\frac{8 t^3}{3}-\pi t^2+1 & 0<t\leq 1 \\ \frac{1}{12} \left(6 t^4-32 \sqrt{t^2-1} t^2-3 \pi t^2+24 t^2-16\sqrt{t^2-1}+30 t^2 \tan^{-1}\left(\sqrt{t^2-1}\right)-18 t^2 \csc^{-1}(t)+8\right)& 1<t<\sqrt{2} \\ \end{array}$$

We recover

p

(* Out= 19/6 - π *)

The graph of p[t] is

Plot[p[t], {t, 0, Sqrt}, AxesLabel -> {"t", "p(t)"},
PlotLabel ->
"Probability that the distance between two points\nrandomly chosen in the \
unit square is greater than t"] • another nice answer +. I felt so guilty with little effort I had to : ubpdqnmathematica.wordpress.com/2016/03/29/… Mar 29 '16 at 4:02
• @ubpdqn Thanks for the hint to your beautiful page. BTW I started to study the same question for the unit circle, but haven't yet found an analytic solution. Mar 29 '16 at 7:41
• Why not post that as an interesting question here ... Mar 29 '16 at 13:03
• @wolfies Thanks for the suggestion. Done. Mar 29 '16 at 17:49