La cuestión es que Maxwell-Boltzmann hace para cualquier sistema. Las partículas cuánticas, sin embargo, no sólo son idénticas, sino que completamente indistinguibles. Para tales partículas, la distribución se caracteriza únicamente por el número de ocupación.
Consideremos un sistema de partículas idénticas que no interactúan. Supongamos que las partículas individuales tienen un Hamiltoniano conocido $H$ con estados propios $|k\rangle$ tal que
$$H|k\rangle=E_k|k\rangle.$$
Si tenemos un montón de decir, $N\gg1$ de estas partículas, escribimos el estado cuántico de todo el sistema en una base diferente (la base de Fock), dada por el número de partículas en estado $k, n_k$ :
$$|n_1,n_2,...\rangle,$$
Y esta base caracteriza completamente nuestro estado cuántico. Se trata de estados propios del Hamiltoniano de muchos cuerpos, con
$$H|n_1,n_2,...\rangle=\left(\sum\limits_kn_kE_k\right)|n_1,n_2,...\rangle.$$
Supongamos entonces que nuestro conjunto de partículas tiene una probabilidad de estado
$$p(n_1,n_2,...).$$
La estadística Maxwell-Boltzmann se genera maximizando la Entropía de Gibbs: $$ S=-\sum\limits_{n_1,n_2,...}p(n_1,n_2,...)\log p(n_1,n_2,...),$$
Dónde $\log$ es el logaritmo natural. Entonces podemos hacerlo sujeto a las restricciones de normalización, energía total y número total de partículas:
$$\sum\limits_{n_1,n_2,...}p(n_1,n_2,...)=1,\quad\sum\limits_{n_1,n_2,...}p(n_1,n_2,...)\sum\limits_kn_kE_k=U,\quad\sum\limits_{n_1,n_2,..}p(n_1,n_2,...)\sum\limits_kn_k=N.$$
La condición de maximización corresponde entonces a (utilizando multiplicadores de Lagrange, el cálculo es un poco complicado)
$$p(n_1,n_2,...)=\prod\limits_k\frac{e^{-\beta n_k(E_k-\mu)}}{\sum\limits_{n_k'}e^{-\beta n_k'(E_k-\mu)}},$$
donde $\beta$ y $\mu$ son parámetros correspondientes a la temperatura inversa, $(k_BT)^{-1}$ y potencial químico. Si tomamos la distribución marginal para un único número de ocupación, digamos, $n_j$ sumando todos los demás valores posibles para todos los demás números de ocupación, obtenemos la distribución de probabilidad para el número de ocupación en cuestión:
$$p_j(n)=\frac{e^{-\beta n(E_j-\mu)}}{\sum\limits_{n'}e^{-\beta n'(E_j-\mu)}}.$$
El paso siguiente dado por las estadísticas de Bose-Einstein y Fermi-Dirac consiste ahora en calcular el número de ocupación medio o esperado del estado cuántico $j$ :
$$f_j=\langle n_j\rangle=\sum\limits_nnp_j(n).$$
Un buen truco computacional es escribir $f_j$ como sigue:
$$f_j=\frac{\sum\limits_nne^{-\beta n(E_j-\mu)}}{\sum\limits_{n'}e^{-\beta n'(E_j-\mu)}}$$
$$f_j=\frac{1}{E_j-\mu}\frac{\sum\limits_nn(E_j-\mu)e^{-\beta n(E_j-\mu)}}{\sum\limits_{n'}e^{-\beta n'(E_j-\mu)}}$$ $$f_j=-\frac{1}{E_j-\mu}\frac{\partial}{\partial\beta}\log\left(\sum\limits_ne^{-\beta n(E_j-\mu)}\right).$$
Como tal, si podemos calcular ese logaritmo a la derecha, podemos encontrar $f_j$ . Aquí, tenemos que hacer una suposición.
Bosones:
Para bosones, $n$ puede ser cualquier número entero comprendido entre $0$ a $\infty$ . Como tal, escribimos
$$\sum\limits_{n=0}^\infty e^{-\beta n(E_j-\mu)}=\frac{1}{1-e^{-\beta(E_j-\mu)}},$$ una serie geométrica convergente. Ahora introducimos esto en nuestra expresión para $f_j$ y obtener
$$f_j=\frac{1}{e^{\beta(E_j-\mu)}-1},$$ que es la estadística de Bose-Einstein.
Fermiones
Para los fermiones, $n$ puede ser 0 ó 1. A continuación, escribimos la suma sobre $n$ como
$$1+e^{-\beta(E_j-\mu)},$$
Lo que nos da
$$f_j=\frac{1}{e^{\beta(E_j-\mu)}+1},$$
El número de ocupación de la estadística de Fermi-Dirac.
EDITAR
Me acabo de dar cuenta de que me he metido hasta el cuello en las matemáticas y se me ha olvidado responder a la pregunta. La intuición detrás de todos estos cálculos.
La cuestión es que, en esencia, la estadística de Maxwell-Boltzmann se ocupa de los estados del sistema, y la de Bose-Einstein/Fermi-Dirac se ocupa de los estados de las partículas. Es fácil confundir ambos conceptos, pero son fundamentalmente diferentes. Mientras que M-B nos da la probabilidad de nuestro sistema (todos nuestros $N$ partículas) para estar en una configuración específica, B-E/F-D nos da la número medio de partículas que se encuentran en un estado determinado.
Como, por definición, se trata de promedios, tenemos que promediar todas las configuraciones del sistema. Y eso es exactamente lo que hemos hecho aquí.
También es bueno mencionar que podemos calcular $\mu$ resolviendo la ecuación
$$\sum\limits_jf_j=N.$$
