In terms of statistical fluctuations, stellar population synthesis models are only asymptotically correct in the limit of a large number of stars, where sampling errors become asymptotically small. When dealing with stellar clusters, starbursts, dwarf galaxies or stellar populations within pixels, sampling errors introduce a large dispersion in the predicted integrated properties of these populations. We present here an approximate but generic statistical formalism which allows a very good estimation of the uncertainties and confidence levels in any integrated property, bypassing extensive Monte Carlo simulations, and including the effects of partial correlations between different observables. Tests of the formalism are presented and compared with proper estimates. We derive the minimum mass of stellar populations which is required to reach a given confidence limit for a given integrated property. As an example of this general formalism, which can be included in any synthesis code, we apply it to the case of young (t <= 20 Myr) starburst populations. We show that, in general, the UV continuum is more reliable than other continuum bands for the comparison of models with observed data. We also show that clusters where more than 105 Msun have been transformed into stars have a relative dispersion of about 10% in Q(He+) for ages smaller than 3 Myr. During the WR phase the dispersion increases to about 25% for such massive clusters. We further find that the most reliable observable for the determination of the WR population is the ratio of the luminosity of the WR bump over the Hβ luminosity. A fraction of the observed scatter in the integrated properties of clusters and starbursts can be accounted for by sampling fluctuations.