We investigate whether the stellar initial mass function ( IMF ) is universal , or whether it varies significantly among young stellar clusters in the Milky Way . We propose a method to uncover the range of variation of the parameters that describe the shape of the IMF for the population of young Galactic clusters . These parameters are the slopes in the low and high stellar mass regimes , \gamma and \Gamma , respectively , and the characteristic mass , M _ { ch } . The method relies exclusively on the high mass content of the clusters , but is able to yield information on the distributions of parameters that describe the IMF over the entire stellar mass range . This is achieved by comparing the fractions of single and lonely massive O stars in a recent catalog of the Milky Way clusters with a library of simulated clusters built with various distribution functions of the IMF parameters . The synthetic clusters are corrected for the effects of the binary population , stellar evolution , sample incompleteness , and ejected O stars . Our findings indicate that broad distributions of the IMF parameters are required in order to reproduce the fractions of single and lonely O stars in Galactic clusters . They also do not lend support to the existence of a cluster mass-maximum stellar mass relation . We propose a probabilistic formulation of the IMF whereby the parameters of the IMF are described by Gaussian distribution functions centered around \gamma = 0.91 , \Gamma = 1.37 , and M _ { ch } = 0.41 M _ { \odot } , and with dispersions of \sigma _ { \gamma } = 0.25 , \sigma _ { \Gamma } = 0.60 , and \sigma _ { M _ { ch } } = 0.27 M _ { \odot } around these values .