$$ \int_{0}^{1}\int^{1}_{0}\int_{0}^{1}\int^{1}_{0}\int_{0}^{1}\int^{1}_{0}Median(|x_{1} - x_{2}||y_{1} - y_{2}|, |x_{1} - x_{3}||y_{1} - y_{3}|, |x_{2} - x_{3}||y_{2} - y_{3}|), dx_{1}dx_{2}dx_{3}dy_{1}dy_{2}dy_{3} $$
Calculation attempts:
In Mathematica on a quad-core Macbook Pro, I can't get the integral to arrive at an answer. That doesn't surprise me, because it's a six dimensional integral and my attempts are probably pretty naive. When I attempt a Monte Carlo integration, even by setting the maximum number of points to 10^8 I can't reach a PrecisionGoal of more than 3 digits or so.
So the question is, what creative methods can I use to reduce this integral even further to something more tractable? I'm not sure where to go from here.
Update:
In a comment to this question, Federico Lebron helpfully suggested partitioning the integral regions such that we wouldn't need to use absolute value for calculating area. He also suggested a better way of expressing the median of three values. Likewise Dan Fulea's current answer sketches out a method of solving it by partitioning it into six cases of three regions each, then exploiting symmetries to combine these into the final answer.
I've taken inspiration from these methods and suggestions to partition the possible groups of rectangle areas into three per order relation on $(y_{1}, y_{2}, y_{3})$, of which there are 3! = 6. Then to the order relation $x_{1} > x_{2} > x_{3}$ there correspond six cases $j_{1}, ..., j_{6}$ which are the median of three area calculations each:
j1: y1 > y2 > y3
- (x1 - x2)(y1 - y2)
- (x1 - x3)(y1 - y3)
- (x2 - x3)(y2 - y3)
j2: y1 > y3 > y2
- (x1 - x2)(y1 - y2)
- (x1 - x3)(y1 - y3)
- (x2 - x3)(y3 - y2)
j3: y2 > y1 > y3
- (x1 - x2)(y2 - y1)
- (x1 - x3)(y1 - y3)
- (x2 - x3)(y2 - y3)
j4: y2 > y3 > y1
- (x1 - x2)(y2 - y1)
- (x1 - x3)(y3 - y1)
- (x2 - x3)(y2 - y3)
j5: y3 > y1 > y2
- (x1 - x2)(y1 - y2)
- (x1 - x3)(y3 - y1)
- (x2 - x3)(y3 - y2)
j6: y3 > y2 > y1
- (x1 - x2)(y2 - y1)
- (x1 - x3)(y3 - y1)
- (x2 - x3)(y3 - y2)
This removes the requirement to calculate absolute value by ensuring that each of the three areas in each case $j_{1}, ..., j_{6}$ is non-negative. To proceed I will integrate on each region $j_{1}, ..., j_{6}$ and finally, as there are also 3! = 6 order relations on $(x_{1}, x_{2}, x_{3})$, I will assume that the final result can be calculated using six cases instead of 36 cases via $6(j_{1} + j_{2} + j_{3} + j_{4} + j_{5} + j_{6})$.
Thus far I've only calculated closed solutions for four of the six regions in Mathematica, so I'll have to report back on whether or not this works correctly. In the meantime, I welcome further answers or criticism!
Final update:
The technique I just described in the first update above worked. In fact, there is actually a further symmetry which can be exploited such that you only have to calculate three cases instead of six. The median areas for $j_{1}$ and $j_{6}$, $j_{2}$ and $j_{5}$ and $j_{3}$ and $j_{4}$ are respectively equal. Therefore we need not calculate all six explicitly, it is sufficient to calculate $6(2j_{1} + 2j_{2} + 2j_{3}) = 12(j_{1} + j_{2} + j_{3})$.
I'm reluctant to answer my own question, so I'm going to mark Dan Fulea's answer as the solution. Thank you!
This is so far a partial answer, we fix a framework with "only" four (instead of six) variables, delimit six cases, and go till the end in five cases. One case remains. I need a better day and a better view for it.
Let $(X_1,Y_1)$, $(X_2,Y_2)$, $(X_3,Y_3)$ be the random variables, ( defined on some hidden probability space $\Omega$, that model the choice of the points. Then we set
$A$, $B$, $A+B$ with $0\le A\le B\le A+B\le 1$ to be the random variables corresponding to the ordered triple starting from the three random variables $|X_i-X_j|$, with $1\le i<j\le 3$,
$S$, $T$, $S+T$ with $0\le S\le T\le S+T\le 1$ to be the random variables corresponding to the ordered triple starting from the three random variables $|Y_i-Y_j|$, with $1\le i<j\le 3$.
We would like to compute the joined distribution of these variables, needed for the computation. Let us do this for $X_1,X_2,X_3$. We restrict to the $\Omega$-subspace $\Omega_{123}$ where $0\le X_1\le X_2\le X_3$. We use instead the variables $X_1$, $A=(X_2-X_1)$, $B=(X_3-X_2)$. The the quantity to be integrated does not depend on $X_1$, but only on $A,B$, and for $(A,B)=(a,b)$ with $(a,b)$ in the simplx $S$ determined by $0\le a,b\le a+b\le 1$, the possible values of $X_1$ are between $0$ and $(1-a-b)$, so $(A,B)$ has the distribution $(1-a-b)\; da\; db$ with $(a,b)$ in the simplex $S$, and so far both $a\le b$ and $b\le a$ are possible. As a check, we note that $$ \begin{aligned} \int_{(a,b)\in S}S (1-a-b)\; da\; db &= \int_{(a,b)\in S} 1\; da\; db - 2\int_{(a,b)\in S}S a\; da\; db \\ &=\frac 12-2\int_{a\in[0,1]}a(1-a)\; da \\ &=\frac 12-2\left(\frac 12-\frac 13\right)=\frac 12-\frac 13=\frac 16\ , \end{aligned} $$ which corresponds to the total probability of $\Omega_{123}$. On the other pieces $\Omega_{ijk}$ given by $0\le X_i\le X_j\le X_k\le 1$ we do the same, use $X_i$, $A=X_j-X_i$ and $B=X_k-X_j$. Finally we paste together a tuple random variable $(A,B)$ distributed on $S$ via $$ 6(1-a-b)\; da\; db\ . $$ There is a symmetry w.r.t. $a\leftrightarrow b$ of the above, so we consider instead an ordered tuple $(A,B)$, distributed on $S'$ defined by $S'=\{\ (a,b)\in[0,1]^2\ :\ 0\le a\le b\le 1-a\le 1\ \}$ via $$ 12(1-a-b)\; da\; db\ . $$ We do the same with $S,T$, and thus use a tuple random variable $(S,T)$, distributed on $S'$ via $$ 12(1-s-t)\; ds\; dt\ . $$ "Unfortunately", we have to split back the situation using $(A,B)$ on the spaces $\Omega_{ijk}$, $ijk$ permutation of $123$, and splitting parallely $$12(1-a-b)\; da\; db=2(1-a-b)\; da\; db+2(1-a-b)\; da\; db+ \dots+(1-a-b)\; da\; db\ , $$ since the order of $(Y_1,Y_2,Y_3)$ has to be "mixed" (equally probable) with any order of $(X_1,X_2,X_3)$. We fix the order of $(Y_1,Y_1,Y_3)$, leading to the fixed order $0\le S\le T\le 1$, but then we have to build the products of $(S,T,S+T)$ with each ordering of $(A,B,A+B)$, as shown in the following table: $$ \begin{array}{|c||c|c|c|c|c|c|} \hline & 1 & 2 & 3 & 4 & 5 & 6\\ \hline \hline \color{blue}S & A & A & B & B & A+B & A+B\\ \color{blue}T & B & A+B & A & A+B & A & B\\ \color{blue}{S+T} & A+B & B & A+B & A & B &A\\\hline \end{array} $$ (The blue side is fixed.) We have six cases, as the columns above are suggesting them.
In case 1, we have to build the median among the three random variables $SA$, $TB$, and $(S+T)(A+B)$. It is the only case with an immediate decision, $TB$ is the median, the corresponding integral is $$ \begin{aligned} J_1 &= \int_{(a,b)\in S'}\int_{(s,t)\in S'} bt\; 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt\ . \\ &=24\left(\int_{(a,b)\in S'}b\;(1-a-b)\; da\; db\right)^2 \\ &=24\left(\int_0^{1/2}da\int_a^{1-a}b\;(1-a-b)\; db\right)^2 \\ &=24\left(\int_0^{1/2}da\; \frac 16(2a-1)^2(a+1)\right)^2 \\ &=24\left( \frac 1{32}\right)^2 \\ &=\frac 3{128}\ . \end{aligned} $$
In case 2, we have to build the median among the three random variables $SA$, $T(A+B)$, and $(S+T)B$. Because the minimal value is $SA$, the median of the three is $\min(T(A+B),\ (S+T)B)=TB+\min(TA,\ SB)$. It may be a good idea to combine this case with the next case. We denote by $J_2$ the corresponding integral.
In case 3, we have to build the median among the three random variables $SB$, $TA$, and $(S+T)(A+B)$. Because the maximal value is $(S+T)(A+B)$, the median of the three is $\max(TA,\ SB)$. We denote by $J_3$ the corresponding integral. Then we have:
$$ \begin{aligned} J_2+J_3 &= \int_{(a,b)\in S'}\int_{(s,t)\in S'} (bt+\min(ta,sb)+\max(ta,sb))\; \\ &\qquad\qquad\qquad\cdot 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt\ . \\ &= \int_{(a,b)\in S'}\int_{(s,t)\in S'} (bt+ta+sb)\; 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt\ . \\ &=J_1+2\cdot 24 \left(\int_{(a,b)\in S'}a\;(1-a-b)\; da\; db\right) \left(\int_{(a,b)\in S'}b\;(1-a-b)\; da\; db\right) \\ &=\frac 3{128}+2\cdot 24\cdot\frac 1{96}\cdot\frac 1{32} =\frac 3{128}+\frac 1{64}= \frac 5{128}\ . \end{aligned} $$
In case 4, we have to build the median among the three random variables $SB$, $TA+TB$, and $TA+SA$. Because the maximal value is $TA+TB$, the median is $\max(SB,\ (S+T)A)$. We denote by $J_4$ the corresponding integral.
In case 5, we have to build the median among the three random variables $SA+SB$, $TA$, and $SB+TB$. Because the maximal value is $SB+TB$, the median is $\max(TA,\ S(A+B))$. We denote by $J_5$ the corresponding integral. The exchange $(a,b)\leftrightarrow(s,t)$ gives $J_4=J_5$. So we compute only one of the two integrals. Let us make algebraically explicit the condition $$ ta\ge sa+sb\ . $$ This is $(t-s)a\ge sb$, i.e. $s\left(\frac ba+1\right)\le t\le 1-s$, so this happens for fixed $a,b$ only when $$ 0\le s\le \left(\frac ba+2\right)^{-1}\ ,\qquad s\left(\frac ba+1\right)\le t\le 1-s\ . $$ The corresponding integral is $$ \begin{aligned} J_5' &= \int_{a=0}^{1/2} \int_{b=a}^{1-a} \int_{s=0}^{1/(b/a+2)} \int_{t=s(b/a+1)}^{1-s} ta\; 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt \\ &=\frac 1{12}\log\frac 23+\frac {11}{288}\approx 0.0044056854354307\dots \ . \end{aligned} $$ Computer force was used. The integral using $(sa+sb)$ instead of $ta$ on the same piece is $$ \begin{aligned} J_5'' &= \int_{a=0}^{1/2} \int_{b=a}^{1-a} \int_{s=0}^{1/(b/a+2)} \int_{t=s(b/a+1)}^{1-s} (sa+sb)\; 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt \\ &=\frac 16\log\frac 23+\frac 5{72}\approx 0.0018669264264170\dots \ . \end{aligned} $$ It is now enough to compute the integral of $(sa+sb)$ on the full piece $S'\times S'$, which is $$ \begin{aligned} J_5''' &= \int_{a=0}^{1/2} \int_{b=a}^{1-a} \int_{s=0}^{1/2} \int_{t=s}^{1-s} (sa+sb)\; 2\cdot 12\; (1-a-b)(1-s-t)\; da\; db\; ds\; dt \\ &=\frac 1{96} \ . \end{aligned} $$ This gives $$ J_4=J_5=J_5'+(J_5'''-J_5'')=\frac 1{12}\log \frac32-\frac1{48}\approx 0.012955425675680\dots \ . $$
In case 6, we have to build the median among the three random variables $S(A+B)=SA+SB$, $TB$, and $A(S+T)=SA+TA$. We denote by $J_6$ the corresponding integral.
Let us determine the median.
Note that this case contains some symmetry, since we pair $(A\le B\le A+B)$ with the reversed $(S+T\ge T\ge S)$. Or conversely, $(S\le T\le S+T)$ with $(A+B\ge B\ge A)$. Then we have either $TA\le SB$ or $SB\le TA$. It is enough to consider only the second case (twice), which is written algebraically $$ sb\le ta\ . $$ So we have the order $SA+SB\le SA+TA$. The value of $TB$ can now be in each of the three positions, before $SA+SB$, in between, or after $SA+TA$. This gives three domains with complicated boundaries.
In the case 6.1 we have to insure algebraically $$ sb\le ta\le tb \le sa+sb\ .$$ This can be reached by letting (in a first approach) $a\in [0,1/2]$, then $b\in[a,1-a]$, then we take $s\in[0, 1/2]$ and search for possible values for $t$, if any. Such a favorable $t$ is constrained to "only" $s\le sb/a\le t\le s(1+a/b)$ and $t\le 1-s$. This restricts first the domain for $(a,b)$, since we need $b/a\le 1+a/b$, and then for $s$ since we need $sb/a\le 1-s$.
The remained cases 6.2 ($sa+sb\le tb\le sa+ta$ and $sb\le ta$) and 6.3 ($sa+ta\le tb$ and $sb\le ta$) are also leading to domains delimited by algebraic curves.
I have to take a breath and a break here, need a better symmetry argument for attacking the case 6, if there is any. I will submit, and possibly come back with the explicit computation next days. (If nobody else does it, and i succeed to compute, against my today expectation...)