15
$\begingroup$

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$ .

Improve this question
$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – Bob Hanlon
    Mar 23, 2016 at 22:26
  • $\begingroup$ @mathe thank you for accepting my answer. However, I think the other answers are more instructive and efficient, e.g. wolfies $\endgroup$
    – ubpdqn
    Mar 28, 2016 at 22:45

5 Answers 5

Reset to default
12
$\begingroup$

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

pdf = UniformDistribution[2];
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.

enter image description here

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ 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) $\endgroup$
    – Quantum_Oli
    Mar 23, 2016 at 11:40
  • 1
    $\begingroup$ You could use pdf = UniformDistribution[2] for short. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Mar 23, 2016 at 13:49
  • $\begingroup$ @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/… $\endgroup$
    – ubpdqn
    Mar 29, 2016 at 4:05
14
$\begingroup$

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:


(source: tri.org.au)

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]]
Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ Well. two seconds and, say, US$ 127 $\endgroup$
    – Dr. belisarius
    Mar 23, 2016 at 18:49
  • $\begingroup$ How much did you pay for your copy of Mathematica? $\endgroup$
    – wolfies
    Mar 23, 2016 at 18:50
  • 2
    $\begingroup$ 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]] $\endgroup$
    – wolfies
    Mar 23, 2016 at 19:08
  • $\begingroup$ @woflies Nice! Do you mind adding it to your answer so I can upvote it? $\endgroup$
    – Dr. belisarius
    Mar 23, 2016 at 19:10
  • 1
    $\begingroup$ Sure. The point was never meant to be about actual code - but the conceptual approach, which works nicely either way here. $\endgroup$
    – wolfies
    Mar 23, 2016 at 19:39
13
$\begingroup$

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...

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ After some testing of my own, that makes two netbooks. :) $\endgroup$
    – J. M.'s lack of A.I. ♦
    Mar 24, 2016 at 0:13
  • $\begingroup$ Looks like Wolfram could improve it code by replacing the difference of two uniforms by a triangular. The computation is much faster. $\endgroup$
    – A.G.
    Mar 28, 2016 at 23:30
  • 1
    $\begingroup$ @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...:) $\endgroup$
    – ubpdqn
    Mar 29, 2016 at 4:09
  • $\begingroup$ @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. $\endgroup$
    – ciao
    Mar 30, 2016 at 8:22
  • $\begingroup$ @ciao best wishes for start up...blog is meaningless musings...play. Have fun :) $\endgroup$
    – ubpdqn
    Mar 30, 2016 at 8:25
8
$\begingroup$

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[2]}}]

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

(* 19/6 - π *)
Improve this answer
$\endgroup$
7
$\begingroup$

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[230]= 1 - π t^2 + (8 t^3)/3 - t^4/2 *)

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

(* Out[231]= 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[1]

(* Out[222]= 19/6 - π *)

The graph of p[t] is

Plot[p[t], {t, 0, Sqrt[2]}, AxesLabel -> {"t", "p(t)"}, 
 PlotLabel -> 
  "Probability that the distance between two points\nrandomly chosen in the \
unit square is greater than t"]

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ another nice answer +. I felt so guilty with little effort I had to : ubpdqnmathematica.wordpress.com/2016/03/29/… $\endgroup$
    – ubpdqn
    Mar 29, 2016 at 4:02
  • $\begingroup$ @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. $\endgroup$
    – Dr. Wolfgang Hintze
    Mar 29, 2016 at 7:41
  • $\begingroup$ Why not post that as an interesting question here ... $\endgroup$
    – wolfies
    Mar 29, 2016 at 13:03
  • $\begingroup$ @wolfies Thanks for the suggestion. Done. $\endgroup$
    – Dr. Wolfgang Hintze
    Mar 29, 2016 at 17:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.