Irrespective of the dark matter (DM) candidate, several potentially observable signatures derive from the velocity distribution of DM in halos, in particular in the Milky Way (MW) halo. Examples include direct searches for weakly-interacting massive particles (WIMPs), p-wave suppressed or Sommerfeld-enhanced annihilation signals, microlensing events of primordial black holes (PBHs), etc. Most current predictions are based on the Maxwellian approximation which is not only theoretically inconsistent in bounded systems, but also not supported by cosmological simulations. A more consistent method sometimes used in calculations for direct WIMP searches relies on the so-called Eddington inversion method, which relates the DM phase-space distribution function (DF) to its mass density profile and the total gravitational potential of the system. Originally built upon the isotropy assumption, this method can be extended to anisotropic systems. We investigate these inversion methods in the context of Galactic DM searches, motivated by the fact that the MW is a strongly constrained system, and should be even more so with the ongoing Gaia survey. We still draw conclusions that apply to the general case. In particular, we illustrate how neglecting the radial boundary of the DM halo leads to theoretical inconsistencies. We also show that several realistic configurations of the DM halo and the MW baryonic content entail ill-defined DFs, significantly restricting the configuration space over which these inversion methods can apply. We propose consistent solutions to these issues. Finally, we compute several observables inferred from constrained Galactic mass models relevant to DM searches (WIMPs or PBHs), e.g. moments and inverse moments of the DM speed and relative speed distributions.