Quantum statistics from the (anti)commutation relations of the operators? From a QFT point of view, the difference between bosons and fermions is that their creation/annihilation operators ($a^{\dagger}$, $a$ and $c^{\dagger}$, $c^{\dagger}$ respectively) obey the following relations:
$$ [a_i, a_j^{\dagger}] = \delta_{ij},$$
$$ \{c_i, c_j^{\dagger}\} = \delta_{ij}.$$
How can I derive/relate this microscopic and single-particle view, with the statistical mechanical Bose-Einstein and Fermi-Dirac distributions?
 A: For a single species of boson/fermion with no interactions, the Hamiltonian is
\begin{align}
  H &=\sum_k \omega_k a_k^\dagger a_k \hskip1cm\text{(boson)} \\
\\
  H &=\sum_k \omega_k c_k^\dagger c_k \hskip1cm\text{(fermion)}
\tag{1}
\end{align}
with
\begin{align}
 a_ia_j - a_j a_i  = 0
\hskip1cm
 a_ia_j^\dagger - a_j^\dagger a_i 
  &=\delta_{ij} \hskip1cm\text{(boson)} \\
\\
 c_ic_j + c_j c_i = 0
\hskip1cm
 c_ic_j^\dagger + c_j^\dagger c_i 
  &=\delta_{ij} \hskip1cm\text{(fermion)}.
\tag{2}
\end{align}
The vacuum state $|0\rangle$, with zero particles, satisfies
\begin{align}
 a_k|0\rangle &=0 \hskip1cm\text{(boson)} \\
\\
 c_k|0\rangle &=0 \hskip1cm\text{(fermion)}
\tag{3}
\end{align}
for all modes $k$. 
Each application of $a_k^\dagger$ or $c_k^\dagger$ to the
vacuum state creates a particle in mode $k$.
The operator
\begin{align}
  N_k &= a_k^\dagger a_k \hskip1cm\text{(boson)} \\
\\
  N_k &= c_k^\dagger c_k \hskip1cm\text{(fermion)}
\tag{4}
\end{align}
counts the number of particles in the $k$-th mode, because
a state $|\psi\rangle$ that satisfies 
$$
N_k|\psi\rangle=n_k|\psi\rangle
\tag{5}
$$
has $n_k$ particles in the $k$-th mode. 
To see this, use equations (2) to deduce
\begin{align}
  N_k a_j^\dagger &= a_j^\dagger (N_k+\delta_{jk}) \hskip1cm\text{(boson)} \\
\\
  N_k c_j^\dagger &= c_j^\dagger (N_k+\delta_{jk}) \hskip1cm\text{(fermion)}.
\tag{5b}
\end{align}
A state that satisfies
$$
 H|\psi\rangle=E_\psi|\psi\rangle
\tag{6}
$$ 
has total energy $E_\psi$.
The adjoint of the lower-left equation
in (2) implies $(c_k^\dagger)^2=0$, so 
$n_k\in\{0,1\}$ for fermions. 
The boson version of equation (2) does
not impose any such restriction,
so $n_k\in\{0,1,2,3,...\}$ for bosons.
Here's how this is used in statistical mechanics.
If the boson/fermion system is in thermal
equilibrium with some other (unmodelled) system, then 
the expectation value of any observable $X$ associated
with the boson/fermion system is
$$
   \rho(X) = \frac{1}{Z}\sum_\psi e^{-\beta E_\psi}
 \frac{\langle \psi|X|\psi\rangle}{
  \langle \psi|\psi\rangle}
\hskip2cm
 Z\equiv \sum_\psi e^{-\beta E_\psi}
\tag{7}
$$
where the sum is over states satisfying (6).
For photons, the sum is over all states satisfying (6).
For a system of matter bosons (or fermions), the sum is typically restricted
to states with a given total number of particles.
The Bose-Einsten and Fermi-Dirac distributions are obtained by
using (7) to calculate $\rho(N_k)$,
the average occupation number in a given mode.
This calculation can be done using the operator identity
$$
   H = \sum_k\omega_k N_k
\tag{8}
$$
to get
$$
  E_\psi = \sum_k\omega_k n_k,
\tag{9}
$$
where $E_\psi$ and $n_k$ are defined by equations (5)-(6). Use this in (7)
to get
$$
   \rho(N_k) = \frac{\sum_\psi e^{-\beta E_\psi}n_k}{
 \sum_\psi e^{-\beta E_\psi}}
\tag{10}
$$
which can also be written
$$
   \rho(N_k) = -\beta^{-1}\frac{\partial}{\partial\omega_k} \log Z
\tag{11}
$$
with the partition function $Z$ defined in (7) regarded as a function
of the energy-coefficients $\omega_k$.
Derivations of the
Bose-Einsten and Fermi-Dirac distributions in 
typical statistical-mechanics books,
such as chapter 9 in Reif's Statistical and Thermal Physics,
start with these ingredients:


*

*equations (9) and (11), which are equations (9.2.1) and (9.2.5) in Reif,
respectively;

*the fact that $n_k$ is unrestricted for bosons,
and restricted to $n_k\in\{0,1\}$ for fermions
(because of equation (2), as mentioned above), which are equations
(9.2.13) and (9.2.15) in Reif;

*the constraint (if any) on the total number of particles $\sum_k n_k$,
which is equation (9.2.14) and (9.2.16) in Reif.
This constraint leads to the "chemical potential," usually denoted $\mu$.
The derivation from this point on is standard, so I won't repeat it here.
