The Milky Way (MW) mass profile is crucial to the Galactic dynamics and many astrophysical studies at various scales. Here we estimate the mass and concentration of the MW halo using the kinematic data of its satellite galaxies, including the latest measurements from Gaia DR2; we also demonstrate how to verify the galaxy formation model in cosmological simulations with observed satellite kinematics. With a simulation-based 6D phase-space distribution function (DF) model of satellite kinematics, we can infer halo properties efficiently, and handle the selection function and measurement errors rigorously. Applying our DF from the EAGLE simulation to 28 satellites, we obtain an MW halo mass of
and a concentration of
. The inferred mass profile is consistent with previous measurements, but with better precision and reliability due to the improved methodology and data. Using the DF model and best-fit MW potential, we provide much more precise estimates of kinematics for those satellites with uncertain measurements. Compared to the EAGLE DF, which matches the MW satellite kinematics very well, the DF from the semi-analytical model (SAM) based on the dark-matter-only simulation over-represents satellites with small radii and velocities. We attribute this difference to the inadequate disruption of satellites with small pericenter distances in the SAM. Our DF method can also be extended to other tracers or other galaxy groups and clusters.