  The variational stochastic self-consistent harmonic approximation is combined
with the calculation of third-order anharmonic coefficients within
density-functional perturbation theory and the "$2n+1$" theorem to calculate
anharmonic properties of crystals. It is demonstrated that in the perturbative
limit the combination of these two methods yields the perturbative phonon
linewidth and frequency shift in a very efficient way, avoiding the explicit
calculation of fourth-order anharmonic coefficients. Moreover, it also allows
calculating phonon lifetimes and inelastic neutron scattering spectra in solids
where the harmonic approximation breaks down and a non-perturbative approach is
required to deal with anharmonicity. To validate our approach, we calculate the
anharmonic phonon linewidth in the strongly anharmonic palladium hydrides. We
show that due to the large anharmonicity of hydrogen optical modes the
inelastic neutron scattering spectra are not characterized by a Lorentzian
line-shape, but by a complex structure including satellite peaks.
