N-body simulations have unveiled several apparently universal properties of dark matter halos, including a cusped density profile, a power-law pseudo-phase-space density rho/sigma(3)(r), and a linear beta-gamma relation between the density slope and the velocity anisotropy. We present a family of self-consistent phase-space distribution functions (DFs) F(E, L), based on the Dehnen-McLaughlin Jeans models, that incorporate these universal properties very accurately. These DFs, derived using a quadratic programming technique, are analytical, positive, and smooth over the entire phase space and are able to generate four-parameter velocity anisotropy profiles beta(r) with arbitrary asymptotic values beta(0) and beta(infinity). We discuss the orbital structure of six radially anisotropic systems in detail and argue that, apart from its use for generating initial conditions for N-body studies, our dynamical modeling provides a valuable complementary approach to understand the processes involved in the formation of dark matter halos.