probability or statistics - Can Mathematica find an expression for the distribution of the median of N i.i.d. random variables?

I posted this question on Math, but there has been silence there since. So, I wonder if anyone here can get any closer to the answer to my question using Mathematica. Here is the question:

Suppose I draw $N$ random variables from independent but identical uniform distributions, where $N$ is an even integer. I now sort the drawn values and find the two middlemost of these. Finally, I calculate a simple average of these two middlemost values.

Is there a closed-form description of the progression of distributions that arise as $N$ increases from $N=2$ to $N=∞$ ? The first distribution is easily found to be Triangular, but what about the rest? Plots from simulations in MATLAB, with a uniform distribution on the range 0 to 1, provide the following illustrations:

Answer

Mathematica does make this pretty easy. The statistic of interest is the typical estimator of the median when the sample size is even. When the sample size is odd the sample median has a beta distribution:

OrderDistribution[{UniformDistribution[{0, 1}], n}, (n + 1)/2]
(* BetaDistribution[(1 + n)/2, 1 + 1/2 (-1 - n) + n] *)

Now for the case when $n$ is even. First find the joint distribution of the middle two order statistics. Then find the distribution of the mean of those two statistics.

n = 6;
od = OrderDistribution[{UniformDistribution[{0, 1}], n}, {n/2, n/2 + 1}];

md = TransformedDistribution[(x1 + x2)/2, {x1, x2} \[Distributed] od];

PDF[md, x]

Plot[Evaluate[PDF[md, x]], {x, 0, 1}]

To obtain the distribution for general $n$ when $n$ is even we have to use some other than TransformedDistribution. We need to integrate the joint density function and treat $0, $x=1/2$, and $1/2 separately.

fltOneHalf = 2 Integrate[(x1^(-1 + n/2) (1 - x2)^(-1 + n/2) n!)/((-1 + n/2)!)^2 /. 
  x2 -> 2 x - x1, {x1, 0, x}, Assumptions -> n > 1 && 0 < x < 1/2]
(* -((4 ((1 - 2 x) x)^(n/2) Gamma[n]*
  Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, x/(-1 + 2 x)])/((-1 + 2 x)*
  Gamma[n/2]^2)) *) 

fOneHalf = 2 Integrate[(x1^(-1 + n/2) (1 - x2)^(-1 + n/2) n!)/((-1 + n/2)!)^2 /.
  x2 -> 1 - x1, {x1, 0, 1/2}, Assumptions -> n > 1]

(* (2^(2 - n) n!)/((-1 + n) ((-1 + n/2)!)^2) *)

(* Because the density is symmetric, we'll take advantage of that *)
fgtOneHalf = FullSimplify[fltOneHalf /. x -> y /. y -> 1 - x]
(* (4 (-1 + (3 - 2 x) x)^(n/2) Gamma[n]*
  Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, (-1 + x)/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2) *)

Putting this together in a single function:

pdf[n_, x_] := 
 Piecewise[{{-((4 ((1 - 2 x) x)^(n/2)*Gamma[n] Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, 

       x/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2)), 0 < x < 1/2},
   {(2^(2 - n) n!)/((-1 + n) ((-1 + n/2)!)^2), x == 1/2},
   {(4 (-1 + (3 - 2 x) x)^(n/2) * Gamma[n]*
 Hypergeometric2F1[1 - n/2, n/2, (2 + n)/2, (-1 + x)/(-1 + 2 x)])/((-1 + 2 x) Gamma[n/2]^2), 
    1/2 < x < 1}}, 0]