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-zspectra 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 H0, H2, 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 Hiiregions.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 H2(~105M⊙) which is predominantly CO-dark and traced by C+rather than C0. Nevertheless, more than 90% of the [Cii] emission is associated with the neutral atomic gas. Our models confirm the necessity to include an X-ray source with an inferred luminosity in good agreement with direct X-ray observations. Finally, we investigate the escape fraction of ionizing photons for different energy ranges. While the escape fraction for the main Hiiregion lies around 50–65%, we show that most of the soft X-ray photons are able to escape and may play a role in the ionization and heating of the circumgalactic or intergalactic medium.Conclusions.Multicomponent ISM models associate a complex enough distribution of matter and phases with a simple enough topological description to be constrained with probabilistic frameworks. Despite ignoring effects such as reflected light, the diffuse radiation field, or ionization by several non-cospatial sources, they remain well adapted to individual Hiiregions and to star-forming galaxies dominated by one or a few Hiiregions, and the improvement due to the combination of several components largely compensates for other secondary effects.