It can be shown (see http://physics.stackexchange.com/a/101456/19976) that the quantum grand canonical partition function for a system of non-interacting, identical bosons is \begin{align} Z = \prod_{i\in I}\sum_{n=0}^\infty x_i^n , \qquad x_i := e^{-\beta(\epsilon_i-\mu)} \end{align} where $I$ is some index set whose elements label single particle energy states (elements of an orthonormal energy eigenbasis of the Hilbert space), $i$ is the index that labels these states, and $\epsilon_i$ is the energy of state $i$.

Now, suppose that for each state $i$, we define a quantity $Z_i$ as follows: \begin{align} Z_i = \sum_{n=0}^\infty x_i^n, \tag{$\star$} \end{align} which we interpret as the "partition function for the $i$th energy eigenstate," then the full partition function for the system can be written as a product of these $Z_i$'s; \begin{align} Z = \prod_{i\in I} Z_i. \end{align} Now, it can also be shown (again see http://physics.stackexchange.com/a/101456/19976) that the ensemble average occupancy $\langle n_i\rangle$ of state $i$ is given by \begin{align} \langle n_i \rangle = x_i \frac{\partial}{\partial x_i} \ln Z. \end{align} Now, notice what happens when we insert the expression for the full partition function $Z$ as a product of the $Z_i$'s into this expression; \begin{align} \langle n_i\rangle &= x_i \frac{\partial}{\partial x_i} \ln \left(\prod_{j\in I} Z_j\right) \\ &= x_i \frac{\partial}{\partial x_i} \sum_{j\in I} \ln Z_j \\ &= x_i \sum_{j\in I} \frac{\partial}{\partial x_i} \ln Z_j \\ &= x_i \frac{\partial}{\partial x_i} \ln Z_i \end{align} where in the last step, we used the fact that each $Z_j$ only depends on the corresponding $x_j$, its partial with respect to $x_i$ vanishes unless $i=j$. In other words, what we have shown is that

It can be shown (see https://physics.stackexchange.com/a/101456/19976) that the quantum grand canonical partition function for a system of non-interacting, identical bosons is \begin{align} Z = \prod_{i\in I}\sum_{n=0}^\infty x_i^n , \qquad x_i := e^{-\beta(\epsilon_i-\mu)} \end{align} where $I$ is some index set whose elements label single particle energy states (elements of an orthonormal energy eigenbasis of the Hilbert space), $i$ is the index that labels these states, and $\epsilon_i$ is the energy of state $i$.

Now, suppose that for each state $i$, we define a quantity $Z_i$ as follows: \begin{align} Z_i = \sum_{n=0}^\infty x_i^n, \tag{$\star$} \end{align} which we interpret as the "partition function for the $i$th energy eigenstate," then the full partition function for the system can be written as a product of these $Z_i$'s; \begin{align} Z = \prod_{i\in I} Z_i. \end{align} Now, it can also be shown (again see https://physics.stackexchange.com/a/101456/19976) that the ensemble average occupancy $\langle n_i\rangle$ of state $i$ is given by \begin{align} \langle n_i \rangle = x_i \frac{\partial}{\partial x_i} \ln Z. \end{align} Now, notice what happens when we insert the expression for the full partition function $Z$ as a product of the $Z_i$'s into this expression; \begin{align} \langle n_i\rangle &= x_i \frac{\partial}{\partial x_i} \ln \left(\prod_{j\in I} Z_j\right) \\ &= x_i \frac{\partial}{\partial x_i} \sum_{j\in I} \ln Z_j \\ &= x_i \sum_{j\in I} \frac{\partial}{\partial x_i} \ln Z_j \\ &= x_i \frac{\partial}{\partial x_i} \ln Z_i \end{align} where in the last step, we used the fact that each $Z_j$ only depends on the corresponding $x_j$, its partial with respect to $x_i$ vanishes unless $i=j$. In other words, what we have shown is that

For a given quantum statistical mechanical system with Hilbert space $\mathcal H$, hamiltonian $\hat H$, and number operator $\hat N$, the grand canonical partition function is defined as follows: \begin{align} Z(\beta, \mu) = \mathrm{tr}(e^{-\beta(\hat H - \mu \hat N)}) \end{align} where the trace is being taken over $\mathcal H$.

Suppose, for example, that we have a system of non-identical, non-interacting bosons. Suppose, for the sake of simplicity, that the single particle Hilbert space is two-dimensional and that the single-particle Hamiltonian has two distinct eigenvalues $E_1$ and $E_2$ with corresponding orthonormal eigenstates $|1\rangle$ and $|2\rangle$. Let $|0\rangle$ denote the state with no particles in it, then the Hilbert space of the system is spanned by the following energy eigenstates: \begin{align} \text{0-particles}:&\qquad|0\rangle\\ \text{1-particle}:&\qquad |1\rangle, |2\rangle\\ \text{2-particles}:&\qquad |1\rangle|1\rangle\\ &\qquad|1\rangle|2\rangle\\ &\qquad|2\rangle|1\rangle\\ &\qquad|2\rangle|2\rangle\\ \text{3-particles}: &\qquad |1\rangle|1\rangle|1\rangle\\ &\qquad |1\rangle|1\rangle|2\rangle\\ &\qquad |1\rangle|2\rangle|1\rangle\\ &\qquad |1\rangle|2\rangle|2\rangle\\ &\qquad |2\rangle|1\rangle|1\rangle\\ &\qquad |2\rangle|1\rangle|2\rangle\\ &\qquad |2\rangle|2\rangle|1\rangle\\ &\qquad |2\rangle|2\rangle|2\rangle\\ &\qquad\vdots \end{align} In fact, the $N$-particle subspace is spanned by $2^N$ distinct states. However, if the particles are identical bosons, then the Hilbert space is spanned only by the symmetric combinations of states in each $N$-particle subspace, so listing the states gives: \begin{align} \text{0-particles}:&\qquad|0\rangle\\ \text{1-particle}:&\qquad |1\rangle, |2\rangle\\ \text{2-particles}:&\qquad |1\rangle|1\rangle\\ &\qquad \tfrac{1}{\sqrt{2}}\big(|1\rangle|2\rangle+ |2\rangle|1\rangle\big)\\ &\qquad|2\rangle|2\rangle\\ \text{3-particles}: &\qquad |1\rangle|1\rangle|1\rangle\\ &\qquad \tfrac{1}{\sqrt{3}}\big(|1\rangle|1\rangle|2\rangle + |1\rangle|2\rangle|1\rangle + |2\rangle|1\rangle|1\rangle\big)\\ &\qquad \tfrac{1}{\sqrt{3}}\big(|1\rangle|2\rangle|2\rangle + |2\rangle|1\rangle|2\rangle + |2\rangle|2\rangle|1\rangle\big)\\ &\qquad |2\rangle|2\rangle|2\rangle\\ &\qquad\vdots \end{align} and it is easy to show that in general, the subspace with $N$ identical bosons is spanned by $N+1$ distinct states.

In general, when the single-particle Hilbert space is greater than two dimensional, the counting may get more difficult, but the moral is the same. You sum over less states when the particles are identical because the Hilbert space only includes appropriately symmetric or anti-symmetric states.

General quantum remarks - identical v. non-identical particles

For a given quantum statistical mechanical system with Hilbert space $\mathcal H$, hamiltonian $\hat H$, and number operator $\hat N$, the grand canonical partition function is defined as follows: \begin{align} Z(\beta, \mu) = \mathrm{tr}(e^{-\beta(\hat H - \mu \hat N)}) \end{align} where the trace is being taken over $\mathcal H$.

Suppose, for example, that we have a system of non-identical, non-interacting bosons. Suppose, for the sake of simplicity, that the single particle Hilbert space is two-dimensional and that the single-particle Hamiltonian has two distinct eigenvalues $E_1$ and $E_2$ with corresponding orthonormal eigenstates $|1\rangle$ and $|2\rangle$. Let $|0\rangle$ denote the state with no particles in it, then the Hilbert space of the system is spanned by the following energy eigenstates: \begin{align} \text{0-particles}:&\qquad|0\rangle\\ \text{1-particle}: &\qquad|1\rangle \\ &\qquad|2\rangle\\ \text{2-particles}:&\qquad |1\rangle|1\rangle\\ &\qquad|1\rangle|2\rangle\\ &\qquad|2\rangle|1\rangle\\ &\qquad|2\rangle|2\rangle\\ \text{3-particles}: &\qquad |1\rangle|1\rangle|1\rangle\\ &\qquad |1\rangle|1\rangle|2\rangle\\ &\qquad |1\rangle|2\rangle|1\rangle\\ &\qquad |1\rangle|2\rangle|2\rangle\\ &\qquad |2\rangle|1\rangle|1\rangle\\ &\qquad |2\rangle|1\rangle|2\rangle\\ &\qquad |2\rangle|2\rangle|1\rangle\\ &\qquad |2\rangle|2\rangle|2\rangle\\ &\qquad\vdots \end{align} In fact, the $N$-particle subspace is spanned by $2^N$ distinct states. However, if the particles are identical bosons, then the Hilbert space is spanned only by the symmetric combinations of states in each $N$-particle subspace, so listing the states gives: \begin{align} \text{0-particles}:&\qquad|0\rangle\\ \text{1-particle}:&\qquad |1\rangle \\ &\qquad |2\rangle\\ \text{2-particles}:&\qquad |1\rangle|1\rangle\\ &\qquad \tfrac{1}{\sqrt{2}}\big(|1\rangle|2\rangle+ |2\rangle|1\rangle\big)\\ &\qquad|2\rangle|2\rangle\\ \text{3-particles}: &\qquad |1\rangle|1\rangle|1\rangle\\ &\qquad \tfrac{1}{\sqrt{3}}\big(|1\rangle|1\rangle|2\rangle + |1\rangle|2\rangle|1\rangle + |2\rangle|1\rangle|1\rangle\big)\\ &\qquad \tfrac{1}{\sqrt{3}}\big(|1\rangle|2\rangle|2\rangle + |2\rangle|1\rangle|2\rangle + |2\rangle|2\rangle|1\rangle\big)\\ &\qquad |2\rangle|2\rangle|2\rangle\\ &\qquad\vdots \end{align} and it is easy to show that in general, the subspace with $N$ identical bosons is spanned by $N+1$ distinct states.

In general, when the single-particle Hilbert space is greater than two dimensional, the counting may get more difficult, but the moral is the same. You sum over less states when the particles are identical because the Hilbert space only includes appropriately symmetric or anti-symmetric states.

Bose-Einstein statistics - partition function

It can be shown (see http://physics.stackexchange.com/a/101456/19976) that the quantum grand canonical partition function for a system of non-interacting, identical bosons is \begin{align} Z = \prod_{i\in I}\sum_{n=0}^\infty x_i^n , \qquad x_i := e^{-\beta(\epsilon_i-\mu)} \end{align} where $I$ is some index set whose elements label single particle energy states (elements of an orthonormal energy eigenbasis of the Hilbert space), $i$ is the index that labels these states, and $\epsilon_i$ is the energy of state $i$.

Now, suppose that for each state $i$, we define a quantity $Z_i$ as follows: \begin{align} Z_i = \sum_{n=0}^\infty x_i^n, \tag{$\star$} \end{align} which we interpret as the "partition function for the $i$th energy eigenstate," then the full partition function for the system can be written as a product of these $Z_i$'s; \begin{align} Z = \prod_{i\in I} Z_i. \end{align} Now, it can also be shown (again see http://physics.stackexchange.com/a/101456/19976) that the ensemble average occupancy $\langle n_i\rangle$ of state $i$ is given by \begin{align} \langle n_i \rangle = x_i \frac{\partial}{\partial x_i} \ln Z. \end{align} Now, notice what happens when we insert the expression for the full partition function $Z$ as a product of the $Z_i$'s into this expression; \begin{align} \langle n_i\rangle &= x_i \frac{\partial}{\partial x_i} \ln \left(\prod_{j\in I} Z_j\right) \\ &= x_i \frac{\partial}{\partial x_i} \sum_{j\in I} \ln Z_j \\ &= x_i \sum_{j\in I} \frac{\partial}{\partial x_i} \ln Z_j \\ &= x_i \frac{\partial}{\partial x_i} \ln Z_i \end{align} where in the last step, we used the fact that each $Z_j$ only depends on the corresponding $x_j$, its partial with respect to $x_i$ vanishes unless $i=j$. In other words, what we have shown is that

The ensemble average occupancy of state $i$ can be computed entirely from the partition function $Z_i$ associated with that state.

Now, what does each single state partition function $Z_i$ look like? Well, going back to the definition, we notice that it's simply a geometric series; \begin{align} Z_i = \sum_{n=0}^\infty x_i^n = (e^{-\beta(\epsilon_i-\mu)})^0 + (e^{-\beta(\epsilon_i-\mu)})^1 + (e^{-\beta(\epsilon_i-\mu)})^2 + \cdots \end{align} which is exactly what you wrote. The physical interpretation of this is that the index $n$ represents the number of particles occupying the given energy state $i$ with energy $\epsilon_i$, so if we think of this state as its own little system, then its energy levels are $n\epsilon_i$ with no degeneracy, where $n=0,1,2,3, \cdots$.

Maxwell-Boltzmann statistics - the classical limit

When we take the classical limit of the quantum ensemble average occupancies, we obtain the Maxwell-Boltzmann ensemble average occupancy; \begin{align} n_i = e^{\beta(\epsilon_i-\mu)}. \end{align} Aside from the fact that the single-state partition function you wrote down for the Maxwell-Boltzmann distribution yields this limit of the quantum ensemble average occupancies, I'm not sure at present how one can argue that it's "correct." In fact, I'm not entirely sure there is another valid criterion by which to judge whether it's correct.