  We study the distribution function (DF) of dark matter particles in haloes of
mass range 10^{14}--10^{15}\Msun. In the numerical part of this work we measure
the DF for a sample of relaxed haloes formed in the simulation of a standard
\LambdaCDM model. The DF is expressed as a function of energy E and the
absolute value of the angular momentum L, a form suitable for comparison with
theoretical models. By proper scaling we obtain the results that do not depend
on the virial mass of the haloes. We demonstrate that the DF can be separated
into energy and angular momentum components and propose a phenomenological
model of the DF in the form
f_{E}(E)[1+L^{2}/(2L_{0}^{2})]^{-\beta_{\infty}+\beta_{0}}L^{-2\beta_{0}}. This
formulation involves three parameters describing the anisotropy profile in
terms of its asymptotic values (\beta_{0} and \beta_{\infty}) and the scale of
transition between them (L_{0}). The energy part f_{E}(E) is obtained via
inversion of the integral for spatial density. We provide a straightforward
numerical scheme for this procedure as well as a simple analytical
approximation for a typical halo formed in the simulation. The DF model is
extensively compared with the simulations: using the model parameters obtained
from fitting the anisotropy profile, we recover the DF from the simulation as
well as the profiles of the dispersion and kurtosis of radial and tangential
velocities. Finally, we show that our DF model reproduces the power-law
behaviour of phase space density Q=\rho(r)/\sigma^{3}(r).
