Average number of particles in a certain energy level in the Canonical Ensemble A quantum system has $r$ discrete energy levels $\varepsilon_1,\varepsilon_2,\varepsilon_3,...,\varepsilon_r$ and $N$ particles distributed in these levels, with the number of particles at each level denoted by $n_1,n_2,n_3,...,n_r$. I'm trying to find the average number of particles in the $i$-th energy level, $\left\langle n_i\right\rangle$, and the fluctuation of this average, $\left\langle(\Delta n_i)^{2}\right\rangle$, using the Canonical Ensemble.
My attempt
The average energy of the system at the state $R$ determined by the occupation numbers $(n_1,n_2,n_3,...,n_r)_R$ can be computed by
$$
\langle E\rangle=\left\langle E_{R}\right\rangle=\sum_{R} P_{R} E_{R}
=\frac{1}{Z}\sum_{R} E_{R} e^{-\beta E_{R}}
=-\frac{1}{Z}\bigg(\frac{\partial Z}{\partial \beta}\bigg)_{N, V}
=-\bigg(\frac{\partial \ln Z }{\partial \beta}\bigg)_{N, V}
$$
With a similar process, keeping in mind that $E_{R} = \sum_{r} n_r \varepsilon_{r}$, one gets that
$$\langle n_i\rangle = \sum_{R} P_{R} n_i
=\frac{1}{Z}\sum_{R} n_r e^{-\beta \sum_{r} n_i \varepsilon_{r}}
=-\frac{1}{\beta}\bigg(\frac{\partial \ln Z}{\partial \varepsilon_i}\bigg)_{N, V}
$$
Which is supposed to be the correct result. However, I am not sure that this $\langle n_i \rangle = \sum_{R} P_{R} n_i$ is valid for this average since $P_r$ is the probability that the system is in the $R$-state, not that the $r$-th energy level has a certain number of particles...
Is the procedure I have performed in this correct?
 A: This isn't quite right.  A microstate of your system is defined by the $r$-tuple $R=(n_1,n_2,\ldots,n_r)$ which gives the occupation numbers of each energy level.  Each $r$-tuple has a corresponding energy given by $E_R=\sum_{i=1}^r n_{i,r} \epsilon_i$ (where $n_{i,R}$ is the occupation number of the $i^{th}$ energy level in microstate $R$) and the probability that the system occupies each microstate is $P_R = e^{-\beta E_R}/Z$, where $Z$ is the partition function.
It makes sense to compute the average energy of the system via this probability distribution:
$$\left<E\right> = \sum_R P_R E_R = \frac{\sum_R E_R e^{-\beta E_R}}{Z} = -\frac{\partial}{\partial \beta} \log(Z)$$
It doesn't make sense to talk about $\left<E_R\right>$, however. For each microstate $R$, $E_R$ is a fixed number.

The expected number of particles in energy level $i$ can be computed precisely the same way.  We're averaging over all possible microstates, weighted by the probability of that microstate being inhabited by the system:
$$\left<n_i\right> = \sum_R P_R n_{i,R}$$
Expanding this out more,
$$\left<n_i\right> = \frac{1}{Z}\sum_R \exp\left[-\beta \sum_j n_{j,R} \epsilon _j\right]n_{i,R}= \frac{1}{Z}\sum_R -\frac{1}{\beta}\frac{\partial}{\partial \epsilon_i}\exp\left[-\beta\sum_j n_{j,R} \epsilon_j\right]$$
$$= -\frac{1}{\beta}\frac{\partial}{\partial \epsilon_i} \log(Z)$$
