  A novel approach is suggested for the statistical description of quantum
systems of interacting particles. The key point of this approach is that a
typical eigenstate in the energy representation (shape of eigenstates, SE) has
a well defined classical analog which can be easily obtained from the classical
equations of motion. Therefore, the occupation numbers for single-particle
states can be represented as a convolution of the classical SE with the quantum
occupation number operator for non-interacting particles. The latter takes into
account the wavefunctions symmetry and depends on the unperturbed energy
spectrum only. As a result, the distribution of occupation numbers $n_s$ can be
numerically found for a very large number of interacting particles. Using the
model of interacting spins we demonstrate that this approach gives a correct
description of $n_s$ even in a deep quantum region with few single-particle
orbitals.
