Spectral energy distributions (SEDs) encode information about the stellar populations within galaxies. By investigating the properties of these stars, such as their ages, masses, and metallicities, we can gain insights into the underlying physical processes that drive the growth and transformation of galaxies over cosmic time. For this purpose, we explore an amortized implicit inference approach to estimate the posterior distribution of stellar mass, star formation rate, metallicity and redshift of galaxies, using ACS and NIRCam filters. Fed with the MILES stellar population models, we generate a sample of synthetic SEDs to train and test our model. We show that our approach is capable of reliably estimating the properties of an integrated stellar population with, crucially, well-calibrated uncertainties. Once trained, deriving the posteriors is five orders of magnitude times faster than classical MCMC sampling, being able to address a large number of galaxies, and to perform a thick sampling of the posteriors, estimating the deterministic trends and the inherent uncertainty of this highly degenerated inversion problem. As a preliminary work, we apply this method to resolved galaxies from JWST, fitting every pixel to study gradients in the stellar populations properties, showing a good generalization to data. We believe that this machine-learning-based implicit inference framework applied to SED fitting is remarkably promising to deal with the size and complexity of the LSST survey.