The mass of the dark matter halo of the Milky Way can be estimated by fitting analytical models to the phase-space distribution of dynamical tracers. We test this approach using realistic mock stellar halos constructed from the Aquarius N-body simulations of dark matter halos in the $Lambda$CDM cosmology. We extend the standard treatment to include a Navarro-Frenk-White (NFW) potential and use a maximum likelihood method to recover the parameters describing the simulated halos from the positions and velocities of their mock halo stars. We find that the estimate of halo mass is highly correlated with the estimate of halo concentration. The best-fit halo masses within the virial radius, $R_{200}$, are biased, ranging from a 40% underestimate to a 5% overestimate in the best case (when the tangential velocities of the tracers are included). There are several sources of bias. Deviations from dynamical equilibrium can potentially cause significant bias; deviations from spherical symmetry are relatively less important. Fits to stars at different galactocentric radii can give different mass estimates. By contrast, the model gives good constraints on the mass within the half-mass radius of tracers even when restricted to tracers within 60kpc. The recovered velocity anisotropies of tracers, $beta$, are biased systematically, but this does not affect other parameters if tangential velocity data are used as constraints.