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.