Find the pdf of the Sample Median M when n is even. (Order Statistics) I have a problem about order statistics. I am trying to find the probability density function of the median $M=\frac{X_{\left(\frac{n}{2} \right)} + X_{\left(\frac{n}{2} +1\right)}}{2}$ where $n$ is even and  $X_1, X_2, ... ,X_n$ is a random sample from a Uniform$(\theta, \theta + 1)$ population with $\theta > 0$.
This was my attempt but I am not sure if I am in the right track. First, I find the joint of $M$ and $X_{(n/2)}$. The joint pdf of $V=M=\frac{X_{\left(\frac{n}{2} \right)} + X_{\left(\frac{n}{2} +1\right)}}{2}$ and $U=X_{(\frac{n}{2})}$
$$f_{U,V}(u,v) =2 \frac{n!}{\left( \frac{n}{2}-1 \right)!\left( \frac{n}{2}-1 \right)!}\left[(u-\theta)(\theta + 1-2v+u)\right]^{\frac{n}{2}-1}$$
where $ u<v<\frac{u+\theta+1}{2}$ and $\theta < u <\theta+1$. Then I was trying to find the marginal pdf of $V=M$ but the integral is too difficult. They give us the hint that the pdf involves sums and absolute value but this confuses me more. If someone can help me to proceed I will appreciated.
 A: Let's simplify the notation and the problem a little.  First, since the sample count is even, let's write it as $2n.$  Second, since $\theta$ appears as a location parameter, it suffices to solve the problem for a convenient value such as $\theta=0$ and then shift the result for general $\theta.$  Third, to avoid the fraction, let's find the distribution of twice the median and then scale it by a factor of $1/2$ at the end.  Finally, we only need to work out the result for the median between $0$ and $1/2,$ because its distribution is (obviously) symmetric about its midpoint.
In these terms, the joint density of the middle order statistics $(U,V) = (X_{(n)},X_{(n+1)})$ is
$$f_{U,V}(u,v) = \binom{2n}{n-1;1;1;n-1}\, u^{n-1}(1-v)^{n-1}\, \mathcal{I}(0\le u\le v\le 1).$$
More directly put, the density is nonzero only when $0\le u \le v \le 1$ and, in that region, has the values
$$n^2\binom{2n}{n}\, u^{n-1}(1-v)^{n-1}.$$
To find the distribution of $M=U+V,$ change variables to $(M, U)$  and "integrate out" $U$ to find the marginal distribution of $M.$  Before proceeding, note that $V = M-U,$ whence the absolute differential is
$$\mathrm{d}u\,\mathrm{d}v = \mathrm{d}(u)\mathrm{d}(m-u) = \mathrm{d}u\,\mathrm{d}m.$$
In terms of these new variables the differential element becomes
$$\begin{aligned} 
f_{M,V}(m,v)\mathrm{d}m\,\mathrm{d}v &= f_{U,V}\left(u, m-u\right)\,\mathrm{d}u\,\mathrm{d}v\\
&=n^2\binom{2n}{n}u^{n-1}(1-m+u)^{n-1}\,\mathrm{d}u\,\mathrm{d}m
\end{aligned}$$
wherever $0\le u \le m-u \le 1;$ that is,  $\max(0,1-m)\le u \le m/2$ and $0\le m \le 2.$
Let us solve this for $0\le m\le 1.$ We need to compute
$$\begin{aligned}
\int_0^{m/2} u^{n-1}(1-m+u)^{n-1}\,\mathrm{d}u &= \int_0^{m/2} u^{n-1}\sum_{j=0}^{n-1} \binom{n-1}{j}u^j(1-m)^{n-1-j}\,\mathrm{d}u\\
&= \sum_{j=0}^{n-1} (1-m)^{n-1-j}\binom{n-1}{j}\int_0^{m/2} u^{n-1+j}\,\mathrm{d}u\\
&= \sum_{j=0}^{n-1} (1-m)^{n-1-j}\binom{n-1}{j} \frac{(m/2)^{n+j}}{n+j}\\
&=h(m,n). 
\end{aligned}$$
This final sum (a function of $m$ and $n$) doesn't appear to simplify, so I called it $h(m,n)$ for future reference.
After accounting for all the initial simplifications, the density of the median $M$ for $\theta \le m \le \theta + 1/2$ is

$$f_M(m;\theta) = 2n^2\binom{2n}{n}h(2(m-\theta),n)$$

and for $\theta+1/2\le m\le \theta+1$ it is given by $f_M(2\theta+1-m;\theta).$
Here is an example of 40,000 simulated values of $M$ for $n=5$ (sample size $2n=10$) and $\theta=3.$  The plot is a histogram on which the graph of $f_M(;3)$ is shown in red.  It fits well.

What follows is the R code that implements $f_M$ as the function f.M.  It does it in two parts.  f implements the sum for $h(m,n)$ (times some constants) and then f.M applies the initial simplifying steps to reduce the calculation to f.  To avoid problems with overflowing values when $n$ is largish, the calculation is performed in terms of logarithms as much as possible.  The special case where $f_M$ is evaluated exactly at the median is computed separately (but that doesn't really matter, since $f_M$ can be defined any way one likes on such a set of measure zero).
f <- Vectorize(function(m, n) {
  if (isTRUE(m==1)) { # (Optional)
    y <- n^2 * choose(2*n, n) * (1/2)^(2*n-1) / (2*n-1)
    return(y)
  }
  g <- 2 * log(n) + lchoose(2*n, n) # A common factor, potentially very large
  j <- seq_len(n) - 1               # The sum indices from 0 through n-1
  x <- (n+j) * log(m / 2) - log(n+j) + lchoose(n-1, j) + (n-1-j) * log(1-m)
  sum(exp(g + x))
}, "m")

f.M <- function(m, n, theta=0) {
  x <- pmax(0, pmin(1, m - theta)) # Focus on the interval [0,1]
  x <- ifelse(x > 1/2, 1-x, x)     # Apply symmetry
  f(2*x, n)*2                      # Rescale by a factor of 2
}

