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

$\endgroup$
  • $\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 '16 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 '16 at 22:45
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

$\endgroup$
  • 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 '16 at 11:40
  • 1
    $\begingroup$ You could use pdf = UniformDistribution[2] for short. $\endgroup$ – J. M. will be back soon Mar 23 '16 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 '16 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]]
$\endgroup$
  • 2
    $\begingroup$ Well. two seconds and, say, US$ 127 $\endgroup$ – Dr. belisarius Mar 23 '16 at 18:49
  • $\begingroup$ How much did you pay for your copy of Mathematica? $\endgroup$ – wolfies Mar 23 '16 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 '16 at 19:08
  • $\begingroup$ @woflies Nice! Do you mind adding it to your answer so I can upvote it? $\endgroup$ – Dr. belisarius Mar 23 '16 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 '16 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...

$\endgroup$
  • 1
    $\begingroup$ After some testing of my own, that makes two netbooks. :) $\endgroup$ – J. M. will be back soon Mar 24 '16 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 '16 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 '16 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 '16 at 8:22
  • $\begingroup$ @ciao best wishes for start up...blog is meaningless musings...play. Have fun :) $\endgroup$ – ubpdqn Mar 30 '16 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 - π *)
$\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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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