How to derive Fermi-Dirac and Bose-Einstein distribution using canonical ensemble?

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 "Boltzmann distribution" is a classical function $e^{-\beta E}$ on the classical phase space. The classical phase space is not a quantum notion, but a classical one. The quantum state space for a collection of non-interacting Bosonic particles is given by first writing down a basis for one particle states, in terms of the collection $|k\rangle$ of all possible momenta (you can work in a periodic box where this collection is discrete, to avoid the complications of continuously labelled states). Then the most general basis state of the bosonic field is

$$|\psi\rangle = |n(k)\rangle$$

where $n(k)$ is an integer for each allowed $k$, telling you how many particles have this momentum. The canonical ensemble is the statement that each such state occurs with probability proportional to

$$ e^{-\beta E(n(k))} = e^{-\beta \sum_k E(k) n(k)} $$

subject to the constraint of fixed particle number

$$ \sum_k n(k) = N $$

In order to mathematically analyze the constraint of fixed particle number, it is most convenient to introduce a chemical potential, which acts as a Lagrange multiplier in the large $N$ limit, to fix the number of particles. But it is not strictly necessary, because you can consider all the n's restricted by the sum condition.

For fermions, $n(k)$ is either 0 or 1, but all the remaining relations are unchanged. The physics is entirely different--- for large $\beta$, the $N$ particles are forced to occupy the states where $E(k)<E_f$ where $E_f$ is adjusted to make the particle count come out right.

The description of multiparticle states in this way gives the Fock space. In Fock space for free particles, the classical expression $e^{-\beta E}$ gives the probability of being in any one of the orthogonal particle occupation number basis states.