Across all fields of science, statistical modeling often involves simplifying assumptions of functional forms in order to make problems tractable. The advent of modern deep learning techniques, however, has begun to alter this approach by replacing overly simplistic functions with universal function approximators that are fast to train and evaluate. In this work, we exhibit three seemingly disparate applications united under the same framework of neural probabilistic modeling.We will begin by using a neural density estimator - known as a normalizing flow - to model intrinsic quasar continua near Lyman-$\alpha$ given the redward spectrum. We use these predictions to estimate the neutral fraction of hydrogen in the spectrum of two z>7 quasars and apply constraints to the timeline of the Epoch of Reionization. Secondly, we show how to use normalizing flows for identifying stellar streams in data from the Gaia telescope. We use anomaly detection techniques developed for High Energy Physics with limited astrophysical assumptions to re-discover GD-1. Finally, we will demonstrate how to use the approximate Bayesian techniques of simulation-based inference to efficiently sample pMSSM models from an experimentally constrained parameter space. Interestingly, the majority of such models are just outside of current experimental bounds.