We use data of \sim 13.000 stars from the SDSS/APOGEE survey to study the shape of the MDF within the region |l| \leq 11 ^ { \circ } and |b| \leq 13 ^ { \circ } , and spatially constrained to { R _ { GC } \leq 3.5 } kpc .
We apply Gaussian Mixture Modeling ( GMM ) and Non-negative Matrix Factorization ( NMF ) decomposition techniques to identify the optimal number and the properties of MDF components in different spatial locations and under different sampling conditions .
We find the shape and spatial variations of the MDF ( at { [ Fe / H ] \geq - 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 profiles for the three chemical components co-existing in the bulge region .
We confirm the consensus physical interpretation of metal-rich stars as associated with the boxy/peanut X-shape bar , originating from the secular evolution of the early disc .
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 to link a present-day structure with those observed in formation in the center of high redshift galaxies .
We refrain from associating the metal-poor stars with any particular formation mechanism .
They seem to be inconsistent with being thick disc or halo stars , but may be the metal-rich tail of the population currently being characterized 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 .