A new method is described for the calculation of molecular vibrational partition functions and thermal effects on molecular properties including an explicit account of anharmonicity. The approach is based on the vibrational self-consistent-field method. Partition functions and thermal averages of the energies calculated with the new method are generally in good agreement with the result of more accurate methods. At lower temperatures the method gives in addition good results for thermal averages of dipole moments and polarizabilities. The new method is much more efficient than explicit sum-over-states approaches previously used for calculation of thermal averages. Unlike the standard sum-over-states approach, the newly developed method is feasible for larger systems despite the formal exponential increase in the number of states with the size of the system. Thus, it is presently the only practical way for including an explicit treatment of anharmonicity in vibrational wave function based calculations of molecular vibrational partition functions and thermally averaged properties of larger molecules.