How to obtain the Bose-Einstein distribution from the canonical ensemble? 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.
 A: The proof is here https://en.wikipedia.org/wiki/Bose%E2%80%93Einstein_statistics#Derivation_in_the_canonical_approach.
They also have a proof in the setting of the grand canonical ensemble which I think is a much simpler derivation.
A: The expression:
$$\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}$$ can be rewrite in a usual for statphysics manner:
$$\bar{n}_{i}=-\frac{1}{\beta}\frac{\partial\ Ln {\ Z(N)}}{\partial E_i}$$
We denoted denominator as $Z(N)$ -canonical partition function of the system (do not confuse it with $Z_i(N)$ -  partition function with "punctured" i-th one particle state. Then, denoting $z_i=e^{ \beta \mu}$ -"punctured" fugacity, we can write:
$$Z(N)=\sum_{n_i=0}^ {n_i=N}\ e^{-\beta\, n_iE_i} \, Z_i(N-n_i) = \sum_{n_i=0}^ {n_i=N}\ (e^{-\beta\,E_i})^{n_i} (\frac{1}{z_i})^{N-n_i} =$$
 $$=z_i^{-N}\sum_{n_i=0}^ {n_i=N}\ (z_ie^{-\beta\,E_i})^{n_i} $$ Then, sequentially we  write down the sum formulae for a finite geometric series,  take the logarithm,   tend $N$ to infinity, and differentiate by $E_i$. Finally, acting in the wiki spirit, we replace  $\mu_i$ =$\mu$
PS. After I wrote this answer, I looked at the link in Wikipedia, it turns out the derivation of the B-E distribution by a similar method is in the same link ((Reif 1965) on p. 342
