I don't really see the answer in the other answer so let me do the calculation here. Your general Boltzmann Ansatz says that the probability of a state $n$ depends on its energy as
$$ p_n = C \exp(-\beta E_n) $$
where $\beta = 1/kT$. Fermions are identical particles that, for each "box" or one-particle state they can occupy (given e.g. by $nlms$ in the case of the Hydrogen atom-like states), admit either $N=0$ or $N=1$ particles in it. Higher numbers are forbidden by the Pauli exclusion principle. The energies of the multi-particle state with $N=1$ and $N=0$ in a particular one-particle state $nlms$ differ by $\epsilon$. Consequently,
$$ \frac{p_1}{p_0} = \frac{C\exp(-\beta (E+\epsilon))}{\exp(-\beta E)} = \exp(-\beta \epsilon) $$
where I used the Boltzmann distribution. However, the probabilities that the number of particles in the given one-particle state is equal to $N=0$ or $N=1$ must add to one,
$$ p_0 + p_1 = 1.$$
These conditions are obviously solved by
$$ p_0 = \frac{1}{1+\exp(-\beta\epsilon)}, \qquad p_1 = \frac{\exp(-\beta\epsilon)}{1+\exp(-\beta\epsilon)}, $$
which implies that the expectation value of $n$ is equal to the right formula for the Fermi-Dirac distribution:
$$\langle N \rangle = p_0\times 0 + p_1 \times 1 = p_1= \frac{1}{\exp(\beta\epsilon)+1} $$
The calculation for bosons is analogous except that the Pauli exclusion principle doesn't restrict $N$. So the number of particles (indistinguishable bosons) in the given one-particle state may be $N=0,1,2,\dots $. For each such number $N$, we have exactly one distinct state (because we can't distinguish the particles). The probability of each such state is called $p_n$ where $n=0,1,2,\dots$.
We still have
$$\frac{p_{n+1}}{p_n} = \exp(-\beta\epsilon) $$
and
$$ p_0 + p_1 + p_2 + \dots = 1 $$
These conditions are solved by
$$ p_n = \frac{\exp(-n\beta\epsilon)}{1+\exp(-\beta\epsilon)+\exp(-2\beta\epsilon)+\dots } $$
Note that the ratio of the adjacent $p_n$ is what it should be and the denominator was chosen so that all the $p_n$ from $n=0,1,2\dots$ sum up to one.
The expectation value of the number of particles is
$$ \langle N \rangle = p_0 \times 0 + p_1 \times 1 + p_2\times 2 + \dots $$
because the number of particles, an integer, must be weighted by the probability of each such possibility. The denominator is still inherited from the denominator of $p_n$ above; it is equal to a geometric series that sums up to
$$ \frac{1}{1-q} = \frac{1}{1-\exp(-\beta\epsilon)} $$
Don't forget that $1-\exp(-\beta\epsilon)$ is in the denominator of the denominator, so it is effectively in the numerator.
However, the numerator of $\langle N \rangle$ is tougher and contains the extra factor of $n$ in each term. Nevertheless, the sum is analytically calculable:
$$ \sum_{n=0}^\infty n \exp(-n \beta\epsilon) = - \frac{\partial}{\partial (\beta\epsilon)} \sum_{n=0}^\infty \exp(-n \beta\epsilon) =\dots$$
$$\dots = - \frac{\partial}{\partial (\beta\epsilon)} \frac{1}{1-\exp(-\beta\epsilon)} = \frac{\exp(-\beta\epsilon)}{(1-\exp(-\beta\epsilon))^2} $$
This result's denominator has a second power. One of the copies gets cancelled with the denominator before and the result is therefore
$$ \langle N \rangle = \frac{\exp(-\beta\epsilon)}{1-\exp(-\beta\epsilon)} = \frac{1}{\exp(\beta\epsilon)-1} $$
which is the Bose-Einstein distribution.
You could also obtain another version of the Boltzmann distribution for distinguishable particles by a similar calculation. For such particles, $N$ could take the same values as it did for bosons. However, the multiparticle state with $N$ particles in the one-particle state would be degenerate because the particles are distinguishable. There would actually be $N!$ multiparticle states with $N$ particles in them. The sum would yield a Taylor expansion for the same exponential.
Note added later: the derivation above was for $\mu=0$. When the chemical potential is nonzero, all appearances of $\epsilon$ have to be replaced by $(\epsilon-\mu)$. Of course that one may only talk about a well-defined value of $\mu$ when we deal with a grand canonical potential; it is impossible to derive a formula depending on $\mu$ from one that contains no $\mu$ and assumes it's ill-defined. The derivation above was meant to show that the difficult $1/(\exp\pm 1)$ structures do appear from a simpler $\exp(-\beta E)$ Ansatz because I think it's the only nontrivial thing to be shown while discussing the relations between the Boltzmann and BE/FD distributions. If that derivation proves the same link as the textbook does, then I apologize but I think there is "nothing else" of a similar kind to be proven.
The expression:
$$\begin{equation}
\bar{n}_{i}=\frac{\displaystyle \sum_{n_i} n_{i}\,e^{-\beta\, n_iE_i} \, Z_i(N-n_i)}{\displaystyle \sum_{n_i}\,e^{-\beta\, n_iE_i} \, Z_i(N-n_i)}
\end{equation}$$ can be rewrite in a usual for statphysics manner:
$$\bar{n}_{i}=-\frac{1}{\beta}\frac{\partial\ Ln {\ Z(N)}}{\partial E_i}$$
We denoted denominator as $Z(N)$ -canonical partition function of the system (do not confuse it with $Z_i(N)$ - partition function with "punctured" i-th one particle state. Then, denoting $z_i=e^{ \beta \mu}$ -"punctured" fugacity, we can write:
$$Z(N)=\sum_{n_i=0}^ {n_i=N}\ e^{-\beta\, n_iE_i} \, Z_i(N-n_i) = \sum_{n_i=0}^ {n_i=N}\ (e^{-\beta\,E_i})^{n_i} (\frac{1}{z_i})^{N-n_i} =$$
$$=z_i^{-N}\sum_{n_i=0}^ {n_i=N}\ (z_ie^{-\beta\,E_i})^{n_i} $$ Then, sequentially we write down the sum formulae for a finite geometric series, take the logarithm, tend $N$ to infinity, and differentiate by $E_i$. Finally, acting in the wiki spirit, we replace $\mu_i$ =$\mu$
PS. After I wrote this answer, I looked at the link in Wikipedia, it turns out the derivation of the B-E distribution by a similar method is in the same link ((Reif 1965) on p. 342
Best Answer
The derivation of the Fermi-Dirac distribution using the density matrix formalism proceeds as follows:
The setup.
We assume that the single-particle hamiltonian has a discrete spectrum, so the single-particle energy eigenstates are labeled by an index $i$ which runs over some finite or countably infinite index set $I$. A basis for the Hilbert space of the system is the occupation number basis \begin{align} |\mathbf n\rangle = |n_0, n_1, \dots\rangle \end{align} where $n_i$ denotes the number of particles occupying the single-particle energy eigenstate $i$. For a system of non-interacting identical fermions, the set $\mathscr N_-$ of admissible occupation sequences $\mathbf n$ consists of those sequences with each $n_i$ equal to either $0$ or $1$. Let $H$ be the hamiltonian for such a system, and let $N$ be the number operator, then we have \begin{align} H|\mathbf n\rangle = \left(\sum_{i\in I}n_i\epsilon_i\right)|\mathbf n\rangle, \qquad N|\mathbf n\rangle = \left(\sum_{i\in I} n_i\right) |\mathbf n\rangle \end{align} where $\epsilon_i$ is the energy of eigenstate $i$. We can also define an observable $N_i$ which tells us the occupation number of the $i^\mathrm{th}$ single-particle energy state; \begin{align} N_i|\mathbf n\rangle = n_i|\mathbf n\rangle \end{align}
Note that we are attempting to determine the ensemble average occupation number of the $j^\mathrm{th}$ energy eigenstate. In the density matrix formalism, this is given by \begin{align} \langle n_j\rangle =\mathrm{tr}(\rho N_i) \end{align} where \begin{align} \rho = \frac{e^{-\beta(H-\mu N)}}{Z}, \qquad Z = \mathrm {tr}\big(e^{-\beta(H-\mu N)}\big) \end{align}
The proof.