In this article we perform a critical assessment of different known methods for the analytical continuation of bosonic functions, namely the maximum entropy method, the non-negative least-square method, the non-negative Tikhonov method, the Pade approximant method, and a stochastic sampling method. Three functions of different shape are investigated, corresponding to three physically relevant scenarios. They include a simple two-pole model function and two flavours of the non-interacting Hubbard model on a square lattice, i.e. a single-orbital metallic system and a two-orbitals insulating system. The effect of numerical noise in the input data on the analytical continuation is discussed in detail. Overall, the stochastic method by Mishchenko et al. [Phys. Rev. B textbf{62}, 6317 (2000)] is shown to be the most reliable tool for input data whose numerical precision is not known. For high precision input data, this approach is slightly outperformed by the Pade approximant method, which combines a good resolution power with a good numerical stability. Although none of the methods retrieves all features in the spectra in the presence of noise, our analysis provides a useful guideline for obtaining reliable information of the spectral function in cases of practical interest.