This study presents a data-driven, nonparametric framework for modeling the Initial Mass Function (IMF) in young massive clusters (YMCs), emphasizing the role of unresolved binaries and observational biases in shaping the stellar mass spectrum. While recent Bayesian hierarchical approaches ([22]) have demonstrated the importance of propagating measurement uncertainties to constrain IMF slopes across cluster environments, such models often rely on explicit likelihood functions. Here, a distribution-free framework integrating Kernel Density Estimation (KDE) and Approximate Bayesian Computation with Markov Chain Monte Carlo (ABC–MCMC) are implemented to infer stellar mass distributions directly from Gaia DR3 observations without assuming any parametric form. The methodology incorporates an escape-mass criterion to distinguish between single and unresolved binary systems, followed by KDE-based mass density estimation and ABC–MCMC simulation for reconstructing the IMF. Validation through the Friedman–Rafsky multivariate runs test and Linear-Time Maximum Mean Discrepancy confirms close agreement between simulated and observed distributions. The approach yields a higher binary mass fraction (88%) compared to previous parametric studies, offering a flexible, uncertainty-aware alternative for IMF reconstruction. The framework complements Bayesian hierarchical models by providing a scalable, likelihood-free route for nonparametric inference of stellar mass functions in young clusters.