Context. Spectroscopic observations of high-redshift galaxies slowly reveal the same complexity of the interstellar medium (ISM) as expected from resolved observations in nearby galaxies. While providing, in principle, a wealth of diagnostics concerning galaxy evolution, star formation, or the nature and influence of compact objects, such high-z spectra are often spatially and spectrally unresolved, and inferring reliable diagnostics represents a major obstacle. Bright, nearby, unresolved galaxies observed in the optical and infrared domains provide many constraints to design methods to infer ISM properties, but they have so far been limited to deterministic methods and/or with simple topological assumptions (e.g., single 1D model). Aims. It is urgent to build upon previous ISM multiphase and multicomponent methods by using a probabilistic approach that makes it possible to derive probability density functions for relevant parameters while also enabling a large number of free parameters with potential priors. The goal is to provide a flexible statistical framework that is agnostic to the model grid and that considers either a few discrete components defined by their parameter values and/or statistical distributions of parameters. In this paper, we present a first application with the objective to infer probability distributions of several physical parameters (e.g., the mass of H 0 , H 2 , escape fraction of ionizing photons, and metallicity) for the star-forming regions of the metal-poor dwarf galaxy I Zw 18 in order to confirm the low molecular gas content and high escape fraction of ionizing photons from H ii regions. Methods. We present a Bayesian approach to model a suite of spectral lines using a sequential Monte Carlo method provided by the Python package PyMC which combines several concepts such as tempered likelihoods, importance sampling, and independent Metropolis-Hastings chains. The algorithm, provided by the associated code MULTIGRIS, accepts a few components which can be represented as sectors around one or several stellar clusters, or continuous (e.g., power-law, normal) distributions for any given parameter. We applied this approach to a grid of models calculated with the photoionization and photodissociation code Cloudy in order to produce topological models of I Zw 18.Results. The statistical framework we present makes it possible to consider a large number of spectroscopic tracers, with the extinction and systematic uncertainties as potential additional random variables. We applied this technique to the galaxy I Zw 18 in order to reproduce and go beyond previous topological models specifically tailored to this object. While our grid is designed for global properties of low-metallicity star-forming galaxies, we were able to calculate accurate values for the metallicity, number of ionizing photons, masses of ionized and neutral hydrogen, as well as the dust mass and the dust-to-gas mass ratio in I Zw 18. We find a relatively modest amount of H 2 (∼ 10 5 M ) which is predominantly CO-dark and trac...