I 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 Msun), pre-main sequence young stellar populations. This methodology is designed specifically for application to young stellar populations suffering significant differential extinction that can exceed 10 visual magnitudes, 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 lack 4.5 µm excess from circumstellar dust disks. Over 100 stars from our sample in the lightly-obscured Trumpler 14 and 16 clusters have available spectroscopic effective temperatures 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 spectroscopic temperatures. The model parameter distributions of stellar mass and evolutionary age reveal significant variation in the star formation history among four different sub-region groups within the CNC and across an older, 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. This methodology will next be applied to ~30 other Galactic giant molecular clouds as part of a program to produce a new empirical calibration for obscured star formation rate diagnostics. Future large, high-resolution near-infrared datasets obtained with JWST and WFIRST will extend the application of this methodology to denser and/or more distant young massive stellar clusters.