A variant of variational perturbation theory is developed for a calculation of the quantum-statistical density matrix in thermal equilibrium. For the diagonal elements, which represent the particle distribution, higher-order corrections are calculated, and show excellent agreement with numerical results from Schr"odinger theory.