Sampling fluctuations in stellar populations give rise to dispersion 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 multimodality 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 the interpretation of these quantities considerably in terms of age or evolutionary stages. These biases can only be assessed by Monte Carlo simulations. We find that biases are usually negligible when the effective number of stars, , which contribute to a given observable is larger than 10. Bimodal distributions may appear when is between 10 and 0.1. Predictions from any model of stellar population synthesis become extremely unreliable for small 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 Msolar, corresponds to about 105 Msolar (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.