We present a novel infrared spectral energy distribution (SED) modeling methodology that uses likelihood-based weighting of the model fitting results to construct probabilistic H-R diagrams (pHRD) for X-ray identified, intermediate-mass (2-8 $M_{\odot}$), pre-main sequence young stellar populations. This methodology is designed specifically for application to young stellar populations suffering strong, differential extinction ($\Delta 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 $\mu$m SEDs of 2269 likely stellar members that exhibit no excess emission from circumstellar dust disks at 4.5 $\mu$m or shorter wavelengths. A subset of ${\sim}100$ intermediate-mass stars in the lightly-obscured Trumpler 14 and 16 clusters have available spectroscopic $T_{\rm eff}$, measured from the Gaia-ESO survey. We correctly identify the stellar temperature in 70% of cases, and the aggregate pHRD for all sources returns the same peak in the stellar age distribution as obtained using the spectroscopic $T_{\rm 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 ${\sim}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.