  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.
