  Dark Matter (DM) direct detection experiments usually assume the simplest
possible 'Standard Halo Model' for the Milky Way (MW) halo in which the
velocity distribution is Maxwellian. This model assumes that the MW halo is an
isotropic, isothermal sphere, hypotheses that are unlikely to be valid in
reality. An alternative approach is to derive a self-consistent solution for a
particular mass model of the MW (i.e. obtained from its gravitational
potential) using the Eddington formalism, which assumes isotropy. In this paper
we extend this approach to incorporate an anisotropic phase-space distribution
function. We perform Bayesian scans over the parameters defining the mass model
of the MW and parameterising the phase-space density, implementing constraints
from a wide range of astronomical observations. The scans allow us to estimate
the precision reached in the reconstruction of the velocity distribution (for
different DM halo profiles). As expected, allowing for an anisotropic velocity
tensor increases the uncertainty in the reconstruction of f(v) but the
distribution can still be determined with a precision of a factor of 4-5. The
mean velocity distribution resembles the isotropic case, however the amplitude
of the high-velocity tail is up to a factor of 2 larger. Our results agree with
the phenomenological parametrization proposed in Mao et al. (2013) as a good
fit to N-body simulations (with or without baryons), since their velocity
distribution is contained in our 68% credible interval.
