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. Comment: 16 pages, 14 figures. Accpeted for publication in MNRAS