  The mass of the Milky Way is a critical quantity which, despite decades of
research, remains uncertain within a factor of two. Until recently, most
studies have used dynamical tracers in the inner regions of the halo, relying
on extrapolations to estimate the mass of the Milky Way. In this paper, we
extend the hierarchical Bayesian model applied in Eadie & Juri\'c (2019) to
study the mass distribution of the Milky Way halo; the new model allows for the
use of all available 6D phase-space measurements. We use kinematic data of halo
stars out to $142~{\rm kpc}$, obtained from the H3 Survey and $\textit{Gaia}$
EDR3, to infer the mass of the Galaxy. Inference is carried out with the
No-U-Turn sampler, a fast and scalable extension of Hamiltonian Monte Carlo. We
report a median mass enclosed within $100~{\rm kpc}$ of $\rm M(<100 \; kpc) =
0.69_{-0.04}^{+0.05} \times 10^{12} \; M_\odot$ (68% Bayesian credible
interval), or a virial mass of $\rm M_{200} = M(<216.2_{-7.5}^{+7.5} \; kpc) =
1.08_{-0.11}^{+0.12} \times 10^{12} \; M_\odot$, in good agreement with other
recent estimates. We analyze our results using posterior predictive checks and
find limitations in the model's ability to describe the data. In particular, we
find sensitivity with respect to substructure in the halo, which limits the
precision of our mass estimates to $\sim 15\%$.
