The next generation 21 cm surveys open a new window onto the early stages of cosmic structure formation and provide new insights about the Epoch of Reionization (EoR). However, the non-Gaussian nature of the 21 cm signal along with the huge amount of data generated from these surveys will require more advanced techniques capable to efficiently extract the necessary information to constrain the Reionization History of the Universe. In this paper we present the use of Bayesian Neural Networks (BNNs) to predict the posterior distribution for four astrophysical and cosmological parameters. Besides achieving state-of-the-art prediction performances, the proposed methods provide accurate estimation of parameters uncertainties and infer correlations among them. Additionally, we demonstrate the advantages of Normalizing Flows (NF) combined with BNNs, being able to model more complex output distributions and thus capture key information as non-Gaussianities in the parameter conditional density distribution for astrophysical and cosmological dataset. Finally, we propose novel calibration methods employing Normalizing Flows after training, to produce reliable predictions, and we demonstrate the advantages of this approach both in terms of computational cost and prediction performances.