How to derive non-degenerate Bose-Einstein statistics in the canonical ensemble? The derivation in Wikipedia assumes high degeneracy:

Let $w(n,g)$ be the number of ways of distributing $n$ particles among the $g$ sublevels of an energy level.
(...)
The number of ways that a set of occupation numbers $n_i$ can be realized is the product of the ways that each individual energy level can be populated:
$$W=\prod_i w(n_i,g_i)=\prod_i{\frac{(n_i+g_i-1)!}{n_i!(g_i-1)!}}\approx \prod_i{\frac {(n_i+g_i)!}{n_i!(g_i-1)!}}$$
where the approximation assumes that $n_i\gg 1$.

I think the "$n_i\gg 1$" is a typo and they actually meant $g_i\gg 1$. If $g_i = 1$, for example, then $W = 1$.
 A: You have to remember that Bose-Einstein statistics only arises from the grand canonical ensemble, i.e. for a system where the energy and the number of particles are not fixed ; only their mean energy $\langle E \rangle $ and mean number of particles $\langle N \rangle $ are fixed respectively by the temperature $T$ and the chemical potential $\mu$ of the thermostat.
However, one can show that, at the thermodynamic limit, at thermodynamic equilibrium, all statistical ensembles are equivalent since the relative particle number fluctuations $\Delta N$ drop to zero as $N\rightarrow +\infty$.
So what you are trying to do here is to derive the Bose-Einstein statistics in the canonical ensemble, which is only possible in the thermodynamic limit. This is the reason of the approximation "$n_i\gg 1$", which actually stands for "at the thermodynamic limit".

EDIT : How to perform properly the approximation 

As stated in the wikipedia article, we have :
$$
W=\prod_i\frac{(n_i+g_i-1)!}{n_i!(g_i-1)!}
$$
In the case $\forall\,i,\;g_i=1$, we have :
$$
W=\prod_i\frac{n_i!}{n_i!0!}=1
$$
with the convention $0!=1$.
This obviously stays true in the limit $n_i\gg 1$. You always have to take the limit $n_i\gg 1$ at the end of all calculations, i.e. after taking $g_i=1$. 
