We use data of ∼13-=000 stars from the Sloan Digital Sky Survey/Apache Point Observatory Galactic Evolution Experiment survey to study the shape of the bulge metallicity distribution function (MDF) within the region |ℓ| ≤ 11° and |b| ≤ 13°, and spatially constrained to RGC ≤ 3.5 kpc. We apply Gaussian mixture modelling and non-negative matrix factorization decomposition techniques to identify the optimal number and the properties of MDF components. We find that the shape and spatial variations of the MDF (at [Fe/H] ≥ -1 dex) are well represented as a smoothly varying contribution of three overlapping components located at [Fe/H] = +0.32, -0.17, and -0.66 dex. The bimodal MDF found in previous studies is in agreement with our trimodal assessment once the limitations in sample size and individual measurement errors are taken into account. The shape of the MDF and its correlations with kinematics reveal different spatial distributions and kinematical structure for the three components co-existing in the bulge region. We confirm the consensus physical interpretation of metal-rich stars as associated with the secularly evolved disc into a boxy/peanut X-shape bar. On the other hand, metal-intermediate stars could be the product of in-situ formation at high redshift in a gas-rich environment characterized by violent and fast star formation. This interpretation would help us to link a present-day structure with those observed in formation in the centre of high-redshift galaxies. Finally, metal-poor stars may correspond to the metal-rich tail of the population sampled at lower metallicity from the study of RR Lyrae stars. Conversely, they could be associated with the metal-poor tail of the early thick disc.