  Sampling fluctuations in stellar populations give rise to dispersions in
observables when a small number of sources contribute effectively to the
observables. This is the case for a variety of linear functions of the spectral
energy distribution (SED) in small stellar systems, such as galactic and
extragalactic HII regions, dwarf galaxies or stellar clusters. In this paper we
show that sampling fluctuations also introduce systematic biases and
multi-modality in non-linear functions of the SED, such as luminosity ratios,
magnitudes and colours. In some cases, the distribution functions of rational
and logarithmic quantities are bimodal, hence complicating considerably the
interpretation of these quantities in terms of age or evolutionary stages.
These biases can be only assessed by Monte Carlo simulations. We find that
biases are usually negligible when the effective number of stars, Neff, which
contribute to a given observable is larger than 10. Bimodal distributions may
appear when Neff is between 10 and 0.1. Predictions from any model of stellar
population synthesis become extremely unreliable for small Neff values,
providing an operational limit to the applicability of such models for the
interpretation of integrated properties of stellar systems. In terms of stellar
masses, assuming a Salpeter Initial Mass Function in the range 0.08 -- 120 Mo,
Neff=10 corresponds to about 10**5 Mo (although the exact value depends on the
age and the observable). This bias may account, at least in part, for claimed
variations in the properties of the stellar initial mass function in small
systems, and arises from the discrete nature of small stellar populations.
