  This work presents a Bayesian algorithm to fit the recombination and
collisionally excited line spectra of gas photoionized by clusters of young
stars. The current model consists in fourteen dimensions: two electron
temperatures, one electron density, the extinction coefficient, the optical
depth on the $HeI$ recombination lines and nine ionic species. The results are
in very good agreement with those previously published using the traditional
methodology. The probabilistic programming library PyMC3 was chosen to explore
the parameter space via a NUTs sampler. These machine learning tools provided
excellent convergence quality and speed. The primordial helium abundance
measured from a multivariable regression using oxygen, nitrogen and sulfur was
$Y_{P,\,O-N-S}=0.243\pm0.005$ in agreement with a standard Big Bang scenario.
