This is a very cumbersome derivation. It's much easier in terms of 2nd quantization. Let's take free non-relativistic fermions (ideal gas). The grand canonical operator for thermal equilibrium (as derived from the maximum-entropy principle under the constraint of given mean energy and particle number) is given by
$$\hat{R}=\frac{1}{Z} \exp(-\beta \hat{H}+\alpha \hat{N}),$$
where I choose convenient Lagrange multipliers. Their relation to the usual thermodynamic quantities follow later when we've evaluated the partition sum. It will turn out that ##\beta=1/(k_\text{B} T)## and ##\alpha=\mu/T## (where ##T## is the temperature and ##\mu## the chemical potential of the gas).
We assume the particles to be in a cubic box of length ##L##. Assuming periodic boundary conditions, the possible single-particle momentum eigenvalues are ##\vec{p} \in \frac{\hbar}{2 \pi L} \mathbb{Z}^3##. In addition we have the ##2s+1## spin states with ##\sigma_z \in \{-s,-s+1,\ldots,s \}##, where ##s \in \{1/2,3/2,\ldots \}## is the spin of the fermions. The possible occupation numbers for each momentum-spin particle states are only ##0## or ##1## (because we've fermions, and the field operators thus obey canonical equal-time anti-commutation relations). Thus the partion sum is
$$Z=\mathrm{Tr} \exp(-\beta \hat{H}+\alpha \hat{N})=\prod_{\vec{p},\sigma_z} \sum_{N(\vec{p},\sigma_z)=0}^{1} \exp[-N(\vec{p},\sigma_z) E(\vec{p})+\alpha N(\vec{p},\sigma_z)]= \prod_{\vec{p},\sigma_z} [1+\exp(-\beta E(\vec{p})+\alpha)].$$
Thus we have
$$\hat{R}=\frac{\exp(-\beta \hat{H}+\alpha \hat{N})}{\prod_{\vec{p},\sigma_z} [1+\exp(-\beta E(\vec{p})+\alpha)]}.$$
Now the mean occupation number of each single-particle state is given by
$$\mathrm{Tr}[\hat{R} \hat{N}(\vec{p},\sigma_z)]=\frac{\exp[-\beta E(\vec{p})+\alpha]}{1+\exp(-\beta E(\vec{p})+\alpha)} = \frac{1}{\exp[\beta E(\vec{p})-\alpha]+1},$$
which is the Fermi-Dirac distribution.