It's useful to consider the following function, the Grand Partition function:
ξ(T,N) = ΣN Σs(N) exp((Nμ-εs)/T)
We can then note that if we define the probabilty as the following via examination of the Boltzmann factor:
P(N1,ε1) = exp((N1μ-ε1)/T)/ξ
If we sum P over s, we get ξ/ξ, which is indeed 1, and thus this is a valid probability function.
Then, using our knowledge of statistical averages, we can consider thermal averages, and calculate energy and various other properties using the Partition function. In general, given the thermodynamic property X, its thermodynamic property is:
X = <x> = (ΣN Σs(N) xsexp((Nμ-εs)/T)) / ξ