11
$\begingroup$

I hope I can explain this problem clearly. It has been bothering me for sometime and I do not know how to solve it even with Mathematica's help. I finally stumbled upon an idea, see my answer below. The old adage, "even a bllnd squirrel finds an acorn now and then" sure fits.

Right isosceles triangle E (the large red triangle) has base = 1 and height = 1 and the base and height are parallel to the x and y axis.

Inside E there are 2 smaller isosceles right triangles. A ( the yellow triangle ) and B ( the blue triangle ), both of these smaller triangles have a base that is 1 / 2 and a height of 1 / 2 and their bases and heights are parallel to the x and y axes.

If triangles A and B are randomly placed inside E and are wholly inside E, and a random point is uniformly chosen in E, find the probability that the point lies inside the intersection of A and B using Mathematica.

enter image description here

My simulation which is slow but simple is

n = 10000; (* Warning slow code ahead, keep this small *)
t = Triangle[{{0, 0}, {1, 0}, {0, 1}}];
h = {1/2, 0};
g = {0, 1/2};

ans =
  Table[
   {r1, r2, r3, r4} = RandomReal[{0, 1}, 4];
   d = Sqrt[r1] (1 - r2) h + r2 Sqrt[r1] g;
   i = Sqrt[r3] (1 - r4) h + r4 Sqrt[r3] g;
   t1 = Triangle[{d, d + g, d + h}];
   t2 = Triangle[{i, i + h, i + g}];
   p1 = Point[RandomPoint[t]];
   RegionIntersection[t1, t2, p1], {n}];

(n - Count[ans, EmptyRegion[2]])/n // N

this yields a probability of about 1 / 10. The answer given by the poser of the problem is 5 / 36 so I am not even sure about my simulation anymore. I would like to see a faster simulation than mine and one that gets the right answer, or maybe clever use of the Probability command. Is this possible?

I now suggest that you view the answer 5 / 36 with suspicion as pointed out to me by Dr. belisarius.

Improve this question
$\endgroup$
4
  • $\begingroup$ I haven't followed your code in detail yet ... but are you sure that the small triangles are inside the big one? $\endgroup$
    – Dr. belisarius
    Feb 25, 2016 at 18:00
  • $\begingroup$ @Dr. belisarius I have generated lots of them and I am pretty sure. At this point I am not certain about anything. $\endgroup$
    – bobbym
    Feb 25, 2016 at 18:05
  • 1
    $\begingroup$ To clarify: are the corners of the smaller triangles supposed to be also uniformly drawn from within the triangle, or do they have a different distribution altogether? $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 25, 2016 at 18:28
  • $\begingroup$ @J.M. Good point +1. They are uniformly drawn from within the triangle. $\endgroup$
    – bobbym
    Feb 25, 2016 at 18:30

6 Answers 6

Reset to default
7
$\begingroup$

It seems you're right.

In fact you don't need to run the Monte Carlo (point shooting) simulation for getting the result. It's enough with calculating (also via simulation) the mean quotient between the area of the intersection and the area of the large triangle.

The following "calculates" the mean area of the intersection as of Mma V9:

Graphics`Mesh`MeshInit[];
big = {{0, 0}, {1, 0}, {0, 1}};
small = {{0, 0}, {1/2, 0}, {0, 1/2}};
rand := Module[{a},
               While[ManhattanDistance[{0, 0}, a = RandomReal[{0, 1/2}, 2]] > .5,]; a]
randTrans := TranslationTransform[rand]

Table[PolygonArea @@ Graphics`Mesh`PolygonIntersection[{Polygon@randTrans[small], 
     Polygon@randTrans[small]}], {10000}] // Mean

(* 0.0504425 *)

which is 1/10 of the big triangle area

The same, done with DirichletDistribution, just in case

(t = PolygonArea /@ Graphics`Mesh`PolygonIntersection /@ 
     Map[Polygon[TranslationTransform[#][small]] &, 
      RandomVariate[DirichletDistribution[{1, 1, 1}], {10000, 2}]/ 2, {2}]) // Mean

(* 0.050252 *)

(same result)

Moreover, the distribution of values of the areas is:

Histogram[t, 50]

Mathematica graphics

Which looks suspiciously like a Beta Distribution. So I tried:

fdp = FindDistributionParameters[t/.125, NoncentralBetaDistribution[a, b, c]]; 
Plot[{PDF[SmoothKernelDistribution[t/.125], x], 
      PDF[NoncentralBetaDistribution[a, b, c] /. fdp, x]}, {x, 0, 1}, 
         Filling -> Axis]]

Mathematica graphics

And to confirm that our guess is not too far: 

QuantilePlot[t/ .125, NoncentralBetaDistribution[a, b, c] /. fdp]

Mathematica graphics

Which looks suspiciously nice.

Improve this answer
$\endgroup$
15
  • 1
    $\begingroup$ The best way to approach this problem is to specify the location of the right-angle corner of each small triangle, ${\bf p}_1 = (x_1,y_1)$ and ${\bf p}_2 = (x_2,y_2)$, and calculate the area of overlap as a function of ${\bf p}_1$ and ${\bf p}_2$. Then, specify the ranges for these points (constrained by the large triangle) and then perform a four-dimensional integral over these ranges, all properly normalized. No simulation at all. $\endgroup$
    – David G. Stork
    Feb 25, 2016 at 20:03
  • 1
    $\begingroup$ @DavidG.Stork But I don't have the patience for trying to get the intersection formula :) Be my guest! $\endgroup$
    – Dr. belisarius
    Feb 25, 2016 at 20:14
  • 1
    $\begingroup$ @bobby, you need to make the appropriate replacements: Graphics`PolygonUtils`PolygonIntersection[] and Graphics`PolygonUtils`PolygonArea[]. They were moved to a different context in version 10. $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 25, 2016 at 20:43
  • 1
    $\begingroup$ @Dr. belisarius, clever use of the Manhattan Distance. $\endgroup$
    – bobbym
    Feb 25, 2016 at 21:04
  • 1
    $\begingroup$ Dammit bel, don't keep me in suspense! :P What does Expectation[x, x \[Distributed] (NoncentralBetaDistribution[a, b, c] /. fdp)] return? $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 26, 2016 at 5:09
5
$\begingroup$

The exact answer is 13/120 = 0.108333.

This problem can be solved without using Mathematica. All you have to do is to follow David G. Stork - divide the whole problem into 8 cases via x1,x2 and y1,y2 relative positions and specify the corresponding integration limits.

But if you want to use Mathamatica for this, then why not use the fact that Integrate knows how to handle Max and Min automatically?

The following code gives the right answer:

Integrate[
     1/2*(Min[-y1 + 1/2 + x2 + y2, -y2 + 1/2 + x1 + y1] - Max[x1, x2])
      *(Min[-x1 + 1/2 + x2 + y2, -x2 + 1/2 + x1 + y1] - Max[y1, y2])
      *1/(1/8)*1/(1/8),
     {x1, 0, 1/2},
     {y1, 0, -x1 + 1/2},
     {x2, 0, 1/2},
     {y2, 0, -x2 + 1/2}
 ]/(1/2)

The first two lines of the integral are the lengths of the resulting intersection triangle sides, and the third line is the multiplication of two pdfs of a uniformly distributed point (the left-bottom corner of a small triangle) in a triangle with sides 1/2 and 1/2.

You can further expand this integral by substituting Min and Max with the corresponding integration limits.

EDIT: Sorry, first I wrote 4 cases, but it's actually 8 cases:

1) x1 <= x2
1.1) y1 <= y2
1.1.a) -x1 + 1/2 + x2 + y2 <= -x1 + 1/2 + x2 + y2
1.1.b) -x1 + 1/2 + x2 + y2 > -x1 + 1/2 + x2 + y2
1.2) y1 > y2
1.2.a) -x1 + 1/2 + x2 + y2 <= -x1 + 1/2 + x2 + y2
1.2.b) -x1 + 1/2 + x2 + y2 > -x1 + 1/2 + x2 + y2

and the same for

2) x1 > x2

EDIT 2: Forgot to add 1/2 into the integral (as the area of a right triangle is 1/2*a*b) and another 1/2, which we divide the whole integral by (as we need the probability of a point being inside the intersection, and it will be equal to the area of interest divided by the whole area).

Improve this answer
$\endgroup$
2
  • $\begingroup$ You can assume one triangle is to the right of the other and not have to consider both. Nevertheless +1. $\endgroup$
    – David G. Stork
    Feb 26, 2016 at 0:02
  • $\begingroup$ You are right, as the problem will be symmetrical in this case. However, then we have to explicitly sat that x2 > x1 (or vice versa) in one of the integration limits and to change the pdf of the second point. $\endgroup$
    – Andrew S.
    Feb 26, 2016 at 0:04
3
$\begingroup$

Maybe a little Monte Carlo can help:

n = 1*^5; (* number of sample points *)
big = {{0, 0}, {1, 0}, {0, 1}};
SeedRandom[42, Method -> "Legacy"]; (* for reproducibility *)
While[t1 = {#2, 1 - #1 - #2} & @@ 
      RandomVariate[DirichletDistribution[{1, 1, 1}]];
      ! (-1/2 <= #1 < 1/2 && 0 <= #2 <= 1/2 - #1) & @@ t1];
While[t2 = {#2, 1 - #1 - #2} & @@ 
      RandomVariate[DirichletDistribution[{1, 1, 1}]];
      ! (-1/2 <= #1 < 1/2 && 0 <= #2 <= 1/2 - #1) & @@ t2];
ts1 = {t1, t1 + {1/2, 0}, t1 + {0, 1/2}};
ts2 = {t2, t2 + {1/2, 0}, t2 + {0, 1/2}};
tsi = Graphics`PolygonUtils`PolygonIntersection[Polygon[ts1], Polygon[ts2]][[1, 1]];
rmf = RegionMember[Triangle[tsi]];
pts = {#2, 1 - #1 - #2} & @@@ RandomVariate[DirichletDistribution[{1, 1, 1}], n];
inside = Select[pts, rmf];
outside = Complement[pts, inside];
Graphics[{{Red, Triangle[big]}, {Orange, Point[outside]},
          {Opacity[1/2, Yellow], Triangle[ts1]}, {Opacity[1/2, Blue], Triangle[ts2]},
          {Green, Point[inside]}},
         PlotLabel -> Row[{"p=", N[Length[inside]/n]}]]

who's feeling lucky?

I might edit this post later to show results of multiple simulations...

Improve this answer
$\endgroup$
11
  • $\begingroup$ This is much faster than mine but using it I am not getting anything close to the predicted answer of 1 / 10. $\endgroup$
    – bobbym
    Feb 25, 2016 at 20:36
  • $\begingroup$ I'm not sure why either. If you'll look at the code, the sampling of the points is completely uniform over the triangle (via DirichletDistribution[{1, 1, 1}]), so everything seems to be within the assumptions of your problem. I'll investigate later. $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 25, 2016 at 20:41
  • $\begingroup$ I liked the use of that distribution very much. $\endgroup$
    – bobbym
    Feb 25, 2016 at 20:47
  • $\begingroup$ I've redone my calc using DirichletDistribution (nice idea!). Same result $\endgroup$
    – Dr. belisarius
    Feb 26, 2016 at 2:25
  • $\begingroup$ @Dr. bel, I've actually been trying to figure out an appropriate NProbability[] expression involving DirichletDistribution[] for the past hour or so, but haven't had much luck. I'll need to think about it more... $\endgroup$
    – J. M.'s eventual burnout ♦
    Feb 26, 2016 at 2:27
3
$\begingroup$

PARTIAL ANSWER: anyone should feel free to continue the derivation. (I will, when I have more time.)

Our general approach is to solve the integrals for the area of the yellow region exactly, rather than to do a Monte Carlo or other simulation.

Consider this figure:

enter image description here

We will need the equations for the hypotenuses:

  1. The equation of the hypotenuse of the red triangle can be derived from realizing it is a straight line (of slope $m = -1$) going through the point $(0,1)$. Thus the equation of the red hypotenuse is: $y = -x + 1$.
  2. The equation of the hypotenuse of the blue triangle can be derived from realizing it is a straight line (of slope $m = -1$) going through the point $(x_1, y_1 + 1/2)$. Thus the equation of the blue hypotenuse is: $y = -x + (x_1 + y_1 + 1/2)$.

Without loss of generality, we can assume that triangle 2 (green) is to the right of triangle 1 (blue), i.e., $x_2 > x_1$. [Later we will multiply by a factor of 2 to include the reverse case.] You can see there are just two cases:

  • Case A: $y_2 > y_1$ (shown in the figure)
  • Case B: $y_2 < y_1$

Assume that both small triangles lie within the large red triangle. (We will place limits on the integration variables to ensure this, below.)

Case A:

The position of the top of the yellow triangle is the point on this hypotenuse corresponding to the $x$ value of $x_2$, that is, the point $(x_2, -x_2 + (x_1 + y_1 + 1/2))$. Thus the height of the yellow triangle is $-x_2 + x_1 + y_1 + 1/2 - y_2$. Thus the area of the yellow triangle is half the area of a square of that given height, that is: $A = 1/2 (-x_2 + x_1 + y_1 + 1/2 - y_2)^2$.

Now we want to find the expected area of the yellow triangle and must integrate over all of Case A conditions, and for that we have to be careful about the limits on integration variables. Our integral will be over four variables, $x_1$, $y_1$, $x_2$, and $y_2$ with the following conditions, which must all hold:

  • $0 < x_1 < 1/2$ (the blue triangle may be within the red one)
  • $0< y_1 < 1/2 - x_1$ (the blue triangle is within the red one)
  • $x_2 > x_1$ (the green triangle is to the right of the blue)
  • $y_1< y_2 < 1/2 - x2$ (our Case A condition)

So now we integrate the area under these conditions:

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

(* 1/7680 *)

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

(* 1/20 *)

Case B:

enter image description here

Now the conditions are:

  • $0 < x_1 < 1/2$ (the blue triangle does not extend beyond the red one)
  • $x_2 > x_1$ (the green is to the right of the blue)
  • $x_1 < x_2 < 1/2$ (the blue triangle is within the red one)
  • $0 < y_2 < y_1$ (our Case B condition)

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Re "I'M POSTING AS I DERIVE:", that's cool, that's like Vangelis composing and recording in real time (see YouTube documentary). $\endgroup$
    – alancalvitti
    Feb 25, 2016 at 22:00
  • $\begingroup$ @alancalvitti Or like Schubert composing his eighth symphony $\endgroup$
    – Dr. belisarius
    Feb 25, 2016 at 23:38
  • $\begingroup$ Think this is more of a pure mathematics problem than Mathematica one. Not more, but rather IS. $\endgroup$
    – Andrew S.
    Feb 25, 2016 at 23:51
  • $\begingroup$ The question poser explicitly wanted help with software to address his problem, and included Mathematica code to address it. One could easily solve it stochastically using $C$ or $C++$. One of the many strengths of Mathematica is of course its symbolic algebra and calculus. Thus if the problem can be cast into such a form, it can be solved by Mathematica (but not by $C$ or $C++$...). Yes, it is pure math to cast the problem into the proper form, but we do this all the time on this site, just as we use clever Mathematica functions rather than the "obvious" assumed ones. $\endgroup$
    – David G. Stork
    Feb 25, 2016 at 23:57
2
$\begingroup$

As J.M. and others have observed DirichletDistribution can be used. In the following simulation the random variate is used to generate lower left points of randomly placed triangles. The intersection and its area are calculated with built-in functions. This is not efficient but I present v[250] and stat[10000].

v[n_] := Table[
  Module[{d = RandomVariate[DirichletDistribution[{1, 1, 1}], 2]/2, 
    b = {{0, 0}, {1/2, 0}, {0, 1/2}}, tg, t1, t2, i},
   tg = Function[x, x + # & /@ b];
   {t1, t2} = Triangle /@ (tg /@ d);
   i = RegionIntersection[t1, t2];
   Graphics[{Triangle[2 b], Blue, t1, Red, t2, Yellow, i}, 
    PlotLabel -> Row[{"Probability: ", 2 Area[i]}]]
   ], {n}]
pr[n_] := 
 Table[Module[{d = 
     RandomVariate[DirichletDistribution[{1, 1, 1}], 2]/2, 
    b = {{0, 0}, {1/2, 0}, {0, 1/2}}, tg, t1, t2, i},
   tg = Function[x, x + # & /@ b];
   {t1, t2} = Triangle /@ (tg /@ d);
   i = RegionIntersection[t1, t2];
   2 Area[i]
   ], {n}]
stat[n_] := With[{dat = pr[n]}, Column[{
     Histogram[dat, Automatic, Frame -> True], 
     Row[{"Probability: ", Mean[dat]}]}, Frame -> True]];

enter image description here

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ Looks so good! +1 $\endgroup$
    – bobbym
    Jun 15, 2016 at 16:04
1
$\begingroup$

We can start by trying to improve the simulation I gave in the OP.

AbsoluteTiming[
n = 1000;
ans =
Table[    
{a, b} = RandomVariate[DirichletDistribution[{1, 1, 1}]]/2;
{c, d} = RandomVariate[DirichletDistribution[{1, 1, 1}]]/2;
t1 = Polygon[{{a, b}, {a, b + 1/2}, {a + 1/2, b}}];
t2 = Polygon[{{c, d}, {c, d + 1/2}, {c + 1/2, d}}];
Area[RegionIntersection[t1, t2]], {n}] // Mean]

This does not look very promising, for one thing it is much slower than the one I gave in the OP. But this one can be sped up. We can solve analytically for the area:

ClearAll[a, b, c, d];
t1 = Polygon[{{a, b}, {a, b + 1/2}, {a + 1/2, b}}];
t2 = Polygon[{{c, d}, {c, d + 1/2}, {c + 1/2, d}}];
p = Assuming[{a \[Element] Reals, b \[Element] Reals, 
c \[Element] Reals, d \[Element] Reals, b != d, a != c}, 
Area[RegionIntersection[t1, t2]]];

The new simulation is much faster:

AbsoluteTiming[n = 10000; 
ans =
Table[    
{a, b} = RandomVariate[DirichletDistribution[{1, 1, 1}]]/2;
{c, d} = RandomVariate[DirichletDistribution[{1, 1, 1}]]/2;
p, {n}] // Mean]

The last bonus is that we can use built in commands to get an analytical answer:

ClearAll[a, b, c, d];
p2 = p /. a -> a/2 /. b -> b/2 /. c -> c/2 /. d -> d/2;

That was necessary because we want to draw points from only a small piece of the larger triangle. (the two smaller triangles are half the height and width)

Expectation[p2, {{a, b} \[Distributed] 
DirichletDistribution[{1, 1, 1}], {c, d} \[Distributed] 
DirichletDistribution[{1, 1, 1}]}]

This yields 1 / 20, so the probability of the point being in the intersection of A and B is (1 / 20 ) / ( 1 / 2 ) = 1 / 10.

This was rather surprising to me, I went against the simulation evidence and this and posted at Math.SE. Luckily Joriki was kind enough to provide an analytical solution that backed this answer up.

Improve this answer
$\endgroup$

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.