The performance of different approximate algorithms for computing anharmonic features in vibrational spectra is analysed and compared on model and more realistic systems that present relevant nuclear quantum effects. The methods considered combine approximate sampling of the quantum thermal distribution with classical time propagation and include Matsubara dynamics, path integral dynamics approaches, linearized initial value representation and the recently introduced adaptive quantum thermal bath. A perturbative analysis of these different methods enables to account for the observed numerical performance on prototypes for overtones and combination bands and to draw qualitatively correct trends for the numerical results obtained for Fermi resonances. Our results prove that the unequal performances of these approaches often derive from the method employed to sample initial conditions and not, as usually assumed, from the lack of coherence in the time propagation. Furthermore, as confirmed by the analysis reported in J. Chem. Phys. 130, 194510 (2021), we demonstrate, both via the perturbative approach and numerically, that path integral dynamics methods fail to reproduce the intensities of these anharmonic features and follow purely classical trends with respect to their temperature behaviour. Finally, the remarkably accurate performance of the adaptive quantum thermal bath approach is documented and motivated.