  Since the early 1970s, stellar population modelling has been one of the basic
tools for understanding the physics of unresolved systems from observation of
their integrated light. Models allow us to relate the integrated spectra (or
colours) of a system with the evolutionary status of the stars of which it is
composed and hence to infer how the system has evolved from its formation to
its present stage. On average, observational data follow model predictions, but
with some scatter, so that systems with the same physical parameters (age,
metallicity, total mass) produce a variety of integrated spectra. The fewer the
stars in a system, the larger is the scatter. Such scatter is sometimes much
larger than the observational errors, reflecting its physical nature. This
situation has led to the development in recent years (especially since 2010) of
Monte Carlo models of stellar populations. Some authors have proposed that such
models are more realistic than state-of-the-art standard synthesis codes that
produce the mean of the distribution of Monte Carlo models.
  In this review, I show that these two modelling strategies are actually
equivalent, and that they are not in opposition to each other. They are just
different ways of describing the probability distributions intrinsic in the
very modelling of stellar populations. I show the advantages and limitations of
each strategy and how they complement each other. I also show the implications
of the probabilistic description of stellar populations in the application of
models to observational data obtained with high-resolution observational
facilities. Finally, I outline some possible developments that could be
realized in stellar population modelling in the near future.
  Open your window and take a look at the night sky on a clear night.....
