The estimation of system failure probabilities in presence of uncertainties may be a difficult task when the values involved are very small, so that sampling-based Monte Carlo methods may become computationally impractical, especially if the computer codes used to model the system response require large computational efforts, both in terms of time and memory. In this work, we propose to exploit the Bayesian Monte Carlo (BMC) approach to the estimation of definite integrals for developing a new, efficient algorithm for estimating small failure probabilities. The Bayesian framework allows an effective use of all the information available, i.e. the computer code evaluations and the input uncertainty distributions, and, at the same time, the analytical formulation of the Bayesian estimator guarantees the construction of a computationally lean algorithm. The proposed method is first satisfactorily tested with reference to an analytic, two-dimensional case study of literature, offering satisfactory results; then, it is applied to a realistic case study of a natural convection-based cooling system of a gas-cooled fast reactor, operating under a post-loss-of-coolant accident (LOCA), showing performances comparable to those of other efficient alternative methods of literature.