# 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, 2016 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, 2016 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, 2016 at 11:40
• You could use pdf = UniformDistribution for short. Mar 23, 2016 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, 2016 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, 2016 at 18:49 • How much did you pay for your copy of Mathematica? Mar 23, 2016 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, 2016 at 19:08 • @woflies Nice! Do you mind adding it to your answer so I can upvote it? Mar 23, 2016 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, 2016 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, 2016 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, 2016 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, 2016 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, 2016 at 8:22
• @ciao best wishes for start up...blog is meaningless musings...play. Have fun :) Mar 30, 2016 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, 2016 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, 2016 at 7:41
• Why not post that as an interesting question here ... Mar 29, 2016 at 13:03
• @wolfies Thanks for the suggestion. Done. Mar 29, 2016 at 17:49
{"ft_lang_label":"__label__en","ft_lang_prob":0.7615861,"math_prob":0.9881363,"size":1186,"snap":"2023-40-2023-50","text_gpt3_token_len":404,"char_repetition_ratio":0.0820643,"word_repetition_ratio":0.0,"special_character_ratio":0.37942666,"punctuation_ratio":0.20895523,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.9988686,"pos_list":[0,1,2,3,4],"im_url_duplicate_count":[null,4,null,4,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2023-09-27T19:20:40Z\",\"WARC-Record-ID\":\"<urn:uuid:5e6cfaaf-833d-4d12-bfe5-8a2190a61229>\",\"Content-Length\":\"227679\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:a48572f8-ab30-4f5f-ac74-a7671ee10187>\",\"WARC-Concurrent-To\":\"<urn:uuid:52997a82-0aaf-4a1a-88c7-59b354decaaa>\",\"WARC-IP-Address\":\"104.18.11.86\",\"WARC-Target-URI\":\"https://mathematica.stackexchange.com/questions/110789/solve-this-probability-problem-symbolically\",\"WARC-Payload-Digest\":\"sha1:VBS4VIL5FWAXNHDU7RSJYUZ7HTT5RX22\",\"WARC-Block-Digest\":\"sha1:WYMAN753JDNXKKDFO4BU3NUNWVYSPA5I\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2023/CC-MAIN-2023-40/CC-MAIN-2023-40_segments_1695233510319.87_warc_CC-MAIN-20230927171156-20230927201156-00127.warc.gz\"}"}