Abstract We present a novel infrared spectral energy distribution (SED) modeling methodology that uses likelihood-based weighting of the model fitting results to construct probabilistic Hertzsprung–Russell diagrams (pHRD) for X-ray-identified, intermediate-mass (2–8 M ⊙), pre-main-sequence young stellar populations. This methodology is designed specifically for application to young stellar populations suffering strong, differential extinction (ΔA V > 10 mag), typical of Galactic massive star-forming regions. We pilot this technique in the Carina Nebula Complex (CNC) by modeling the 1–8 μm SEDs of 2269 likely stellar members that exhibit no excess emission from circumstellar dust disks at 4.5 μm or shorter wavelengths. A subset of ∼100 intermediate-mass stars in the lightly obscured Trumpler 14 and 16 clusters have available spectroscopic T eff, measured from the Gaia-ESO survey. We correctly identify the stellar temperature in 85% of cases, and the aggregate pHRD for all sources returns the same peak in the stellar age distribution as obtained using the spectroscopic T eff. The SED model parameter distributions of stellar mass and evolutionary age reveal significant variation in the duration of star formation among four large-scale stellar overdensities within the CNC and a large distributed stellar population. Star formation began ∼10 Myr ago and continues to the present day, with the star formation rate peaking ≲3 Myr ago when the massive Trumpler 14 and 16 clusters formed. We make public the set of 100,000 SED models generated from standard pre-main-sequence evolutionary tracks and our custom software package for generating pHRDs and mass–age distributions from the SED fitting results.