  In the last few years, there has been significant progress in the development
of machine learning methods tailored to astrophysics and cosmology. We have
recently applied one of these, namely, the neural network bundle method, to the
cosmological scenario. Moreover, we showed that in some cases the computational
times of the Bayesian inference process can be reduced. In this paper, we
present an improvement to the neural network bundle method that results in a
significant reduction of the computational times of the statistical analysis.
The novel aspect consists of the use of the neural network bundle method to
calculate the luminosity distance of type Ia supernovae, which is usually
computed through an integral with numerical methods. In this work, we have
applied this improvement to the Hu-Sawicki and Starobinsky $f(R)$ models. We
also performed a statistical analysis with data from type Ia supernovae of the
Pantheon+ compilation and cosmic chronometers. Another original aspect of this
work is the different treatment we provide for the absolute magnitude of type
Ia supernovae during the inference process, which results in different
estimates of the distortion parameter than the ones obtained in the literature.
We show that the statistical analyses carried out with our new method require
lower computational times than the ones performed with both the numerical and
the neural network method from our previous work. This reduction in time is
more significant in the case of a difficult computational problem such as the
ones addressed in this work.
