Grand partition function substate average number of particles Let's say I have a grand partition function with two states $\epsilon_1$ and $\epsilon_2$:
\begin{equation}
Z = \sum\limits_i \exp\left[-\beta\left(\epsilon_1 n_1+\epsilon_2 n_2 - \mu n_i\right)\right],
\end{equation}
where $n_1+n_2 = n_i$. Better formulation is in the comment of @By Symmetry.
I am interested in obtaining the average $n_1$.
WRONG APPROACH:
I can rewrite my partition function as:
\begin{equation}
Z = \sum\limits_i \exp\left[-\beta\left(\epsilon_1 n_{1i}+\epsilon_2 \delta^{-1} n_{1i} - \mu (1+\delta^{-1}) n_{1i}\right)\right],
\end{equation}
where I used $\delta = n_1/n_2$.
If I take partial derivative over $\mu(1+\delta^{-1})$ this does NOT give me the average number of particles in the state $\epsilon_1$, because there is $\delta^{-1}$ in front of the $\epsilon_2$ term. Is there a way around it?
EDIT: $\delta$ is a bad variable as @Roger Vadim mentioned
 A: Add the multiplicity to your equation, given  $n = n_1 + n_2$:
\begin{align*}
Z &= \sum_{n=0}\,\,\,\sum_{n1=0, n2=n-n1}^n \frac{n!}{n_1! n_2!} e^{-\beta\left\{\epsilon_1 n_1+\epsilon_2 n_2 - \mu (n_1+n_2)\right\}}.\\
&= \sum_n \,\sum_{n_1=0}^n \frac{n!}{n_1! n_2!} e^{-\left\{\beta\epsilon_1 \,n_1+\beta\epsilon_2\, n_2 - \beta\mu (n_1+n_2)\right\}}\\
&= \sum_n\,\sum_{n1, n2} \frac{n!}{n_1! n_2!} e^{-\left\{\beta (\epsilon_1 - \mu) \,n_1+\beta (\epsilon_2-\mu)\, n_2\right\}}\\
&= \sum_{n=0}^{\infty}\left\{ e^{-\beta (\epsilon_1 - \mu) } +  e^{-\beta (\epsilon_2 - \mu) }\right\}^n = \sum_n z_1^n = \frac{1}{1-z_1} 
\end{align*}
where $z_1 = e^{-\beta (\epsilon_1 - \mu) } +  e^{-\beta (\epsilon_2 - \mu) }$, and $\ln Z = - \ln \left(1 - z_1 \right).$
The average number $\bar n_1$
\begin{align*}
\bar n_1 &= \frac{1}{Z} \sum_n \, \sum_{n_1=0}^n n_1\,\frac{n!}{n_1! n_2!} e^{-\beta\left\{\epsilon_1 n_1+\epsilon_2 n_2 - \mu (n_1+n_2)\right\}}\\
&= \frac{-1}{Z}\frac{\partial}{\partial (\beta \epsilon_1)}\sum_n \, \sum_{n_1=0}^n \frac{n!}{n_1! n_2!} e^{-\beta\left\{\epsilon_1 n_1+\epsilon_2 n_2 - \mu (n_1+n_2)\right\}} \\
&= \frac{-1}{Z}\frac{\partial Z}{\partial (\beta \epsilon_1)} 
= -\frac{\partial \ln Z}{\partial (\beta \epsilon_1)}\\
&= + \frac{\partial}{\partial (\beta \epsilon_1)}  \ln \left(1 - z_1 \right)\\
&=  \frac{\partial}{\partial (\beta \epsilon_1)}\ln \left\{1- e^{-\beta (\epsilon_1 - \mu) } -  e^{-\beta (\epsilon_2 - \mu) }\right\}\\
&= \frac{e^{-\beta (\epsilon_1 - \mu) } }{1 -e^{-\beta (\epsilon_1 - \mu) } -  e^{-\beta (\epsilon_2 - \mu) } }\\
\end{align*}
