  Synthesis models predict the integrated properties of stellar populations.
Several problems exist in this field, mostly related to the fact that
integrated properties are distributed. To date, this aspect has been either
ignored (as in standard synthesis models, which are inherently deterministic)
or interpreted phenomenologically (as in Monte Carlo simulations, which
describe distributed properties rather than explain them). We approach
population synthesis as a problem in probability theory, in which stellar
luminosities are random variables extracted from the stellar luminosity
distribution function (sLDF). We derive the population LDF (pLDF) for clusters
of any size from the sLDF, obtaining the scale relations that link the sLDF to
the pLDF. We recover the predictions of standard synthesis models, which are
shown to compute the mean of the sLDF. We provide diagnostic diagrams and a
simplified recipe for testing the statistical richness of observed clusters,
thereby assessing whether standard synthesis models can be safely used or a
statistical treatment is mandatory. We also recover the predictions of Monte
Carlo simulations, with the additional bonus of being able to interpret them in
mathematical and physical terms. We give examples of problems that can be
addressed through our probabilistic formalism. Though still under development,
ours is a powerful approach to population synthesis. In an era of resolved
observations and pipelined analyses of large surveys, this paper is offered as
a signpost in the field of stellar populations.
