The Gaia satellite will soon deliver magnitudes, parallaxes, and proper motions of billions of stars, offering an unprecedented detailed view of the structure and dynamics of our Galaxy. Correctly exploiting Gaia is challenging due to the colossal amount of data, the complicated selection effects and noise, and the inaccuracy of current stellar and 3D Galactic models with respect to the precision of the data. I will demonstrate how those issues can be addressed with modern statistical methods, specifically hierarchical data-driven probabilistic models, which allow us to bridge the gap between the statistical power of the data and the physics of interest. I will showcase results obtained from the first Gaia data release, including precise models of the color-magnitude diagram, red-clump stars, and unresolved double and triple systems, and discuss upcoming data releases and the prospects of constructing detailed three-dimensional models the Galaxy.