In the wikipedia page
https://en.wikipedia.org/wiki/Fermi%E2%80%93Dirac_statistics#Canonical_ensemble
the Fermi-Dirac distribution is obtained from the canonical ensemble in the following manner:
The average occupation number is given by
\begin{equation} \bar{n}_{i}=\frac{\displaystyle \sum_{n_1,n_2,...}n_{i}\,e^{-\beta(n_1E_1+n_2E_2+...)}}{\displaystyle \sum_{n_1,n_2,...}e^{-\beta(n_1E_1+n_2E_2+...)}} \end{equation}
where the summation is over the sets of $n_r$ that satisfy $\sum\,n_r=N$, being $N$ the total number of particles in the system. The derivation followed in Wikipedia then separates the state of energy $E_i$ and calls
\begin{equation} \bar{n}_{i}=\frac{\displaystyle \sum_{n_i} n_{i}\,e^{-\beta\, n_iE_i} \, Z_i(N-n_i)}{\displaystyle \sum_{n_i}\,e^{-\beta\, n_iE_i} \, Z_i(N-n_i)} \end{equation}
where $Z_{i}(N-n_i) = \sum_{n_1,n_n,...} \,e^{-\beta\,(n_1E_1+n_2E_2+...)}$. The summation here does not consider the state $n_i$. Since for the Fermi-Dirac case $n_i=0,1$, it is possible to write the summation easily (it has only two terms) to obtain
\begin{equation} \bar{n}_{i}=\frac{\displaystyle 1}{\displaystyle \bigg[Z_i(N)/Z_i(N-1)\bigg]\,e^{-\beta\,E_i} +1} \end{equation}
which will, in turn, result in the Dirac-Fermi distribution once the ratio between the partition functions is related to the chemical potential.
Now, I'm trying to follow the same procedure to obtain the Bose-Einstein distribution but in this case $n_i$ can vary between $0$ and $\infty$ and the summation cannot be done in any direct way I could think of. How can I approach this?
Thank you very much.
