Probability of a random point inside the intersection of two triangles

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]

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

And to confirm that our guess is not too far:

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

Which looks suspiciously nice.

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

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

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