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.