Hyper physics has a page for the energy distribution functions (here), they say that each of the distributions are the probabilites that a particle has a certain energy state E, but other websites like this one say that Fermi-Dirac provides the probability. I interpret this as you can use the distribution function to get the probability, and the distribution itself, $\bar n_{FD}={1\over e^{\beta (\epsilon - \mu)} +1}$, is not the probability. Is this the case? This question is the same for the other two types of distributions.
1 Answer
Start with the grand canonical partition function $Y$ and the microstate $r=(n_{p_1},n_{p_2},...)=\{n_p\}$: \begin{align} Y&=\sum_r\exp\left(-\beta\left(E_r\left(V_rN_r\right)-\mu N_r\right)\right)\\ &=\sum_{n_{p1}=0}^\infty\exp\left(-\beta\left(\epsilon_{p_1}-\mu\right)n_{p_1}\right)\cdot\sum_{n_{p2}=0}^\infty\exp\left(-\beta\left(\epsilon_{p_2}-\mu\right)n_{p_2}\right)\cdot...\\ &=\frac{1}{1-\exp\left(-\beta\left(\epsilon_{p_1}-\mu\right)\right)}\cdot\frac{1}{1-\exp\left(-\beta\left(\epsilon_{p_2}-\mu\right)\right)}\cdot...\\ &=\prod_p\frac{1}{1-\exp\left(-\beta\left(\epsilon_{p}-\mu\right)\right)} \end{align}
Now calculate the mean occupation number with momentum $p_i$:
\begin{align} \bar{n_{p_1}}&=\frac{1}{Y}\sum_rn_{p_1}\exp\left(-\beta\left(E_r-\mu N_r\right)\right)\\ &=...\\ &=\frac{1}{\exp\left(\beta\left(\epsilon_{p_i}-\mu\right)\right)-1} \end{align}
This is the Bose-Einstein statistic, which gives you an idea which occupation number is to be expected with a given energy / momentum. The Fermi-Dirac statistic can be derived accordingly and gives as well a relation between the energy and the mean occupation number. I would not consider this a probability, rather a distribution from which probability can be derived.