Understanding the computation of integrals of PDF's

Question

I have been having problems with the computation of integrals of PDF's.
For example:
The joint probability density function of $n$ order statistics is $$f_{X_{(1)},...,X_{(n)}}(x_1,..., x_n)=n!f(x_1)\cdots f(x_n)$$ To find the density function of $X_{(i)}$ I know that it is possible to use a sort of binomial distribution argument by saying that, in order for $X_{(i)}=x$ to be the $i$th smallest number, $i-1$ of the other random variables have to be smaller than $x$ and $n-i$ of the remaining variables have to be greater than $x$. Because all r.v.'s are independent and identically distributed, the probability any one of them is greater than, or smaller than x is $1-F(x)$ and $F(x)$ respectively. Independence says therefore that the probability should be $\binom {n}{i-1} \binom {n-i+1}{n-i}(1-F(x))^{n-i}F(x)^{i-1}f(x)$
I understand this explanation, however, I would like to derive the formula by integrating the PDF. The problem is that the formula should be $$ f_{X_{(i)}}(x)= \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty}f_{X_{(1)},...,X_{(n)}}(x_1,..., x_n) \mathrm{d}^{n-1} x_{j \ne i}$$ and I think this integral reduces to $$1^{n-1}f(x)$$ (because of the independence of the variables and the fact that PDF's must integrate to $1$ over their whole interval). Obviously, this is wrong, but I do not know how to carry out the integration algorithmically (the way a computer would do it) to arrive at the correct result.
I have run into several other problems involving integrals of PDFs that I have not been able to solve because of issues with the limits and other nuances. As a particular example, I have $$P(X_{(n)}-X_{(1)} <a)= n \int_{-\infty}^{\infty} (F(x_1+a)-F(x_1))^{n-1}f(x_1)\mathrm{d}x_1$$ Where the random variables are independent and uniformly distributed over $(0,1)$. The above integral should reduce to: $$n\int_0^1 (x_1+a-x_1)^{n-1}\mathrm{d}x_1=na^{n-1}$$ However, the text says: $$n\int_0^{1-a}a^{n-1}\mathrm{d}x_1+n\int_{1-a}^1 (1-x_1)^{n-1}\mathrm{d}x_1=n(1-a)a^{n-1} +a^n$$ I do not understand why the integral was split into two pieces, and I do not even know where the integrand of the second piece comes from $(1-x_1)^{n-1}$
My question is: how do we compute these integrals without physical intuition, using only formal rules (as if we were writing a formal proof)?

Answer 1

(Note: We will be using $\chi_i$ to visually distinguish the function argument from the bound variables of integration.)

Because these are order statistics, we require: $x_1<x_2<\ldots<\chi_i<\ldots <x_{n-1}<x_n$, so then the integral of interest is actually:

$$f_{X_{(i)}}(\chi_i) = n! \cdot f(\chi_i) \cdot{\int_{\chi_i}^{\infty}\!\int_{\chi_i}^{x_{n}}\!\iiint_{\chi_i}^{x_{i+2}} f(x_{i+1})\cdots f(x_{n})\operatorname d x_{i+1}\cdots\operatorname d x_{n} \\ \cdot \int_{-\infty}^{\chi_i}\!\int_{x_{1}}^{\chi_i}\!\iiint_{x_{i-2}}^{\chi_i} f(x_{1})\cdots f(x_{i-1})\operatorname d x_{i-1}\cdots\operatorname d x_{1} }$$

Which would be horrendous to compute.   However, since the samples are independent and identically distributed, we can argue from symmetry that:

$$\begin{align}\iiint\limits_{x_1\leqslant x_2\leqslant\ldots\leqslant\chi_i} f(x_1)\cdots f(x_{k-1})\operatorname d x_{k-1}\cdots\operatorname d x_1 &= \frac 1{(k-1)!} \iiint\limits_{x_1\leqslant\chi_i, x_2\leqslant\chi_i, \ldots x_{k-1}\leqslant\chi_i} f(x_1)\cdots f(x_{k-1})\operatorname d x_{k-1}\cdots\operatorname d x_1 \\[1ex] &= \frac 1{(k-1)!} \left(\int_{s\leqslant\chi_i} f(s)\operatorname d s\right)^{k-1} \\[1ex] & = \frac 1{k!} F(\chi_i)^{k-1}\end{align}$$

Also similarly argued for the other factor, thus obtaining the much cleaner expression:

$$f_{X_{(i)}}(\chi_i) = \dfrac{n!}{(n-i)!(i-1)!}\,f(\chi_i)\,F(\chi_i)^{i-1}\,(1-F(\chi_i))^{n-i}$$

---

The argument by symmetry is that as $$\{(a,b): a\leqslant b\leqslant c\}\cup\{(a,b):b\leqslant a\leqslant c\} =\{(a,b): a\leqslant c \land b\leqslant c\}$$

$\tiny(\textrm{... and that }\{(a,b): a=b, b\leqslant c\}\textrm{ is a null set})$

Then it follows that $$ \int_{a\leqslant c}\int_{b\leqslant c} g(a,b)\operatorname d b\operatorname d a=\iint_{a\leqslant b\leqslant c} g(a,b)\operatorname d (a,b) + \iint_{b\leqslant a\leqslant c} g(a,b)\operatorname d (a,b) $$

And further if $g(a,b)=g(b,a)$ then:

$$ \int_{a\leqslant c}\int_{b\leqslant c} g(a,b)\operatorname d b\operatorname d a=2!\iint_{a\leqslant b\leqslant c} g(a,b)\operatorname d (a,b)$$

And so it could be shown by induction that the result we used holds.