We present BayeSED, a general purpose tool for Bayesian analysis of spectral energy distributions (SEDs) using pre-existing model SED libraries or their linear combinations. The artificial neural networks, principal component analysis, and multimodal-nested sampling (MultiNest) techniques are employed to allow the highly efficient sampling of posterior distribution and the calculation of Bayesian evidence. As a demonstration, we apply this tool to a sample of hyperluminous infrared galaxies (HLIRGs). The Bayesian evidence obtained for a pure starburst, a pure active galactic nucleus (AGN), and a linear combination of starburst+AGN models show that the starburst+AGN model has the highest evidence for all galaxies in this sample. The Bayesian evidence for the three models and the estimated contributions of starbursts and AGNs to infrared luminosity show that HLIRGs can be classified into two groups: one dominated by starbursts and the other dominated by AGNs. Other parameters and corresponding uncertainties about starbursts and AGNs are also estimated using the model with the highest Bayesian evidence. We find that the starburst region of the HLIRGs dominated by starbursts tends to be more compact and has a higher fraction of OB stars than that of HLIRGs dominated by AGNs. Meanwhile, the AGN torus of the HLIRGs dominated by AGNs tends to be more dusty than that of HLIRGs dominated by starbursts. These results are consistent with previous researches, but need to be tested further with larger samples. Overall, we believe that BayeSED could be a reliable and efficient tool for exploring the nature of complex systems such as dust-obscured starburst-AGN composite systems by decoding their SEDs.