We devise a framework to compute accurate phonons in molecular crystals even in case of strong quantum anharmonicity. Our approach is based on the calculation of the static limit of the phononic Matsubara Green's function from path integral molecular dynamics simulations. Our method enjoys a remarkably low variance, which allows one to compute accurate phonon frequencies after a few picoseconds of nuclear dynamics, and it is further stabilized by the use of appropriate constrained displacement operators. We applied it to solid hydrogen at high pressure. For phase III, our predicted infrared (IR) and Raman active vibrons agree very well with experiments. We then characterize the crystalline symmetry of phase IV by direct comparison with vibrational data and we determine the character of its Raman and IR vibron peaks.