Machine learning (ML) has often been applied to space weather (SW) problems in recent years. SW originates from solar perturbations and is comprised of the resulting complex variations they cause within the systems between the Sun and Earth. These systems are tightly coupled and not well understood. This creates a need for skillful models with knowledge about the confidence of their predictions. One example of such a dynamical system is the thermosphere, the neutral region of Earth's upper atmosphere. Our inability to forecast it has severe repercussions in the context of satellite drag and collision avoidance operations for objects in low Earth orbit. Even with (assumed) perfect driver forecasts, our incomplete knowledge of the system results in often inaccurate neutral mass density predictions. Continuing efforts are being made to improve model accuracy, but density models rarely provide estimates of uncertainty. In this work, we propose two techniques to develop nonlinear ML models to predict thermospheric density while providing calibrated uncertainty estimates: Monte Carlo (MC) dropout and direct prediction of the probability distribution, both using the negative logarithm of predictive density (NLPD) loss function. We show the performance for models trained on local and global datasets. This shows that NLPD provides similar results for both techniques but the direct probability method has a much lower computational cost. For the global model regressed on the SET HASDM density database, we achieve errors of 11% on independent test data with well-calibrated uncertainty estimates. Using an in-situ CHAMP density dataset, both techniques provide test error on the order of 13%. The CHAMP models (on independent data) are within 2% of perfect calibration for all prediction intervals tested. This model can also be used to obtain global predictions with uncertainties at a given epoch.