In this paper we use The Quijote simulations in order to extract the cosmological parameters through Bayesian Neural Networks. This kind of models has a remarkable ability of estimating the associated uncertainty, which is one of the ultimate goals in the precision cosmology era. We demonstrate the advantages of BNNs for extracting more complex output distributions and non-Gaussianities information from the simulations.