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$.