New representation of Wigner function in the form of path integral in canonical ensemble for finite temperatures is proposed. An explicit analytical expression for Wigner function is obtained in harmonic approximation. Monte-Carlo method is developed for calculating average values of quantum operators and thermodynamic quantities. The method is tested on different models of anharmonic oscillators. Average values of full, kinetic, and potential energies (as well as quantum momentum distribution functions) are calculated for nondegenerate non-ideal system of ions and electrons.