Several general trends have been identified for equilibrated, self-gravitating collisionless systems, such as density or anisotropy profiles. These are integrated quantities which naturally depend on the underlying velocity distribution function (VDF) of the system. We study this VDF through a set of numerical simulations, which allow us to extract both the radial and the tangential VDFs. We find that the shape of the VDF is universal, in the sense that it depends only on two things: namely the dispersion (radial or tangential) and the local slope of the density. Both the radial and the tangential VDFs are universal for a collection of simulations, including controlled collisions with very different initial conditions, radial infall simulation, and structures formed in cosmological simulations.