Three practical methods for computing sums and densities of states of fully coupled anharmonic vibrations are compared. All three methods are based on the standard perturbation theory expansion for the vibrational energy. The accuracy of the perturbation theory expansion is tested by comparisons with computed eigenvalues and/or experimental vibrational constants taken from the literature for three- and four-atom molecules. For a number of examples, it is shown that the Xij terms in the perturbation theory expansion account for most of the anharmonicity, and the Yijk terms also make a small contribution; contributions from the Zijkl terms are insignificant. For molecules containing up to ∼4 atoms, the sums and densities of states can be computed by using nested DO-loops, but this method becomes impractical for larger species. An efficient Monte Carlo method published previously is both accurate and practical for molecules containing 3-6 atoms but becomes too slow for larger species. The Wang-Landau algorithm is shown to be practical and reasonably accurate for molecules containing ∼4 or more atoms, where the practical size limit (with a single computer processor) is currently on the order of perhaps 50 atoms. It is shown that the errors depend mostly on the average number of stochastic samples per energy bin. An automated version of the Wang-Landau algorithm is described. Also described are the effects of Fermi resonances and procedures for deperturbation of the anharmonicity coefficients. Computer codes based on all three algorithms are available from the authors and