Deriving PDF for order statistics

Question

In Sheldon Ross' book, First Course in Probability, he introduces order statistics and goes on to derive a formula for their distribution using the following reasoning:
First, for the order statistics $X_{(1)},...,X_{(n)}$ to take values $x_1\le ... \le x_n$ the random variables $(X_1,...,X_n)$ must equal $(x_{i_1},...,x_{i_n})$ for some permutation $(i_1,...,i_n)$ of $(1,...,n)$.
Then we must calculate $$p=P\{ x_{i_1}-\frac \epsilon 2 < X_1 <x_{i_1}+\frac \epsilon 2,...,x_{i_n}-\frac \epsilon 2 <X_n<x_{i_n}+\frac \epsilon 2\} \approx \epsilon^n f(x_1)\cdots f(x_n)$$ It follows that, because there are $n!$ different ways that the random variables can equal the given values, we have that $$p=n!\epsilon^n f(x_1)\cdots f(x_n)$$ This far, I understand, however, it is when he says:
"Dividing by $\epsilon^n$ and letting $\epsilon \to 0$"
What bewilders me. I understand why we make the approximation through the derivative of the CDF and that we can multiply the probabilities out because the random variables are independent. However, I do not understand why we divide by $\epsilon^n$ to obtain the PDF of the order statistics. Why do we do it?

Answer 1

Let $\nu(dx)=f(x)\,\lambda_1(dx)$ be a probability measure and $\nu_n(d\mathbf{x})= f(x_1)\cdots f(x_n)\,\lambda_n(d\mathbf{x})$, the product measure. The map $T:\mathbf{x}=(x_1,\ldots,x_n)\mapsto (x_{(1)},\ldots,x_{(n)})$, where $x_{(1)}\leq\ldots x_{(k)}\leq\ldots x_{(n)}$, is called the order statistic map.

If $B=\{[y_1,\ldots,y_n]: y_1<\ldots < y_n\}$, then for each permutation $\sigma$ on $\{1,\ldots,n\}$, there is one to one linear map $P_{\sigma}: B\rightarrow \mathbb{R}^n$ such that $T\circ P_{\sigma}=I_n$ the identity on $B$. It is easy to check that $\lambda_n(\mathbb{R}^n\setminus\bigcup_{\sigma\in\Sigma_n}P_\sigma(B))=0$. The maps $P_\sigma$ are represented by matrices with exactly one $1$ in each row and each column and $\det(P_\sigma)=\pm1$. Then, $\lambda_n\circ T^{-1}\ll\lambda_n$ and $$ \begin{align} \frac{\lambda_n\circ T^{-1}}{d\lambda_n}(\mathbf{y}):=f_T(y_1,\ldots,y_n)=\sum_{\sigma\in\Sigma_n}|\det(P_\sigma)| f(y_{\sigma(1)})\cdots f(y_{\sigma(n)})&= n! f(y_1)\cdots f(y_n) \mathbb{1}_{B}(\mathbf{y})\tag{1}\label{order} \end{align} $$