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?

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 = \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$$, 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 = \sum_R P_R n_{i,R}$$
Expanding this out more, $$\left = \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)$$