  We present a hierarchical Bayesian method for estimating the total mass and
mass profile of the Milky Way Galaxy. The new hierarchical Bayesian approach
further improves the framework presented by Eadie, Harris, & Widrow (2015) and
Eadie & Harris (2016) and builds upon the preliminary reports by Eadie et al
(2015a,c). The method uses a distribution function $f(\mathcal{E},L)$ to model
the galaxy and kinematic data from satellite objects such as globular clusters
(GCs) to trace the Galaxy's gravitational potential. A major advantage of the
method is that it not only includes complete and incomplete data simultaneously
in the analysis, but also incorporates measurement uncertainties in a coherent
and meaningful way. We first test the hierarchical Bayesian framework, which
includes measurement uncertainties, using the same data and power-law model
assumed in Eadie & Harris (2016), and find the results are similar but more
strongly constrained. Next, we take advantage of the new statistical framework
and incorporate all possible GC data, finding a cumulative mass profile with
Bayesian credible regions. This profile implies a mass within $125$kpc of
$4.8\times10^{11}M_{\odot}$ with a 95\% Bayesian credible region of
$(4.0-5.8)\times10^{11}M_{\odot}$. Our results also provide estimates of the
true specific energies of all the GCs. By comparing these estimated energies to
the measured energies of GCs with complete velocity measurements, we observe
that (the few) remote tracers with complete measurements may play a large role
in determining a total mass estimate of the Galaxy. Thus, our study stresses
the need for more remote tracers with complete velocity measurements.
