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.

The primordial helium abundance $Y_P$ is calculated using sulphur as metallicity tracer in the classical methodology (with $Y_P$ as an extrapolation of $Y$ to zero metals). The calculated value, $Y_{P,\,S}=0.244\pm0.006$, is in good agreement with the estimate from the Planck experiment, as well as, determinations in the literature using oxygen as the metallicity tracer. The chemical analysis includes the sustraction of the nebular continuum and of the stellar continuum computed from simple stellar population synthesis grids. The $S^{+2}$ content is measured from the near infrared $\left[SIII\right]\lambda\lambda9069\AA,9532\AA$ lines, while an $ICF\left(S^{3+}\right)$ is proposed based on the $Ar^{3+}/Ar^{2+}$ fraction. Finally, we apply a multivariable linear regression using simultaneously oxygen, nitrogen and sulphur abundances for the same sample to determine the primordial helium abundance resulting in $Y_{P-O,\,N,\,S}=0.245\pm0.007$