The question of considering every macroscopic body in the universe is formidable, mostly because on smaller scales (~sub-galactic) you need to worry about more than just gravity. Also, are you only interested in the velocity distribution at the present epoch? The universe is far from static and so it is natural to consider a time dependent distribution.
Neglecting the Hubble expansion, it is possible to answer your question on a cosmological scale. If you have ever seen cosmological n-body simulations of structure formation, such simulations are only possible because we know how to compute a realistic initial position and velocity distribution for the simulation particles. These initial position and velocity fields are then evolved forward in time with gravity, and the velocity distribution is known at all times.
Edit: I'll expand on this answer to be more concrete. According to our best theories of cosmological structure formation, the universe was amazingly homogeneous following the big bang. However, there were small deviations which were amplified by gravity over time to form the massive structures (galaxies, galaxy clusters, filaments, voids, etc.) that we observe today. Fortunately, the universe has gifted us an incredible source of information about the early universe: the Cosmic Microwave Background.
Now, the CMB itself is incredibly uniform with deviations of only 1 part in 100,000. We are able to measure these deviations (in the form of temperature fluctuations relative to the mean temperature) using instruments like WMAP, Planck, etc. The useful information that we extract from the CMB is in the form of a power spectrum, $P(k)$, which tells us how large the fluctuations are over different angular scales. Taking the temperature deviation field and decomposing it into spherical harmonics (since the CMB is projected on the sphere of the sky), we compute how much power there is in each mode.
Roughly speaking, cosmological simulations use the information from the CMB to generate a distribution of particles with the same statistical properties.
Suppose we start with a perfect 3-D lattice of particles. We want to perturb the particles so they have the same statistical properties as the CMB. Our first goal is to compute the dimensionless density fluctuation
\begin{equation}
\delta(\vec{x})\equiv \frac{\rho(\vec{x})-\bar{\rho}}{\bar{\rho}},
\end{equation}
where $\bar{\rho}$ is the mean density.
Starting with a field of Gaussian white noise, $\xi(\vec{x})$, the density fluctuation at each point is given by the convolution of the white noise field with what is known as the transfer function $T(\vec{x})$:
\begin{equation}
\delta(\vec{x})=(\xi * T)(\vec{x})=\int d^3 y \;\xi(\vec{y})T(|\vec{x}-\vec{y}|)
\end{equation}
Note: In practice, these expression are all discrete since we are working on a lattice, not continuous space. The basic ideas are the same.
The transfer function encodes the important statistical properties of the CMB. It is related to the power spectrum: $T(k)\equiv [(2\pi/L)^3P(k)]^{1/2}$.
Once we know the value of the density fluctuations at each point, we want to compute the position and velocity of each particle. We begin by solving for the gravitational field of the density distribution. One way to do this is to differentiate in Fourier space
\begin{equation}
\hat{\Phi}(\vec{k})=\frac{i\vec{k}}{k^2} \hat{\delta}(\vec{k}).
\end{equation}
Then, we utilize the Zeldovich approximation
\begin{equation}
\begin{split}
\vec{x} &= \vec{q} + \Phi(\vec{x}) \\
\dot{\vec{x}} &= \dot{\Phi}(\vec{q}),
\end{split}
\end{equation}
where $\vec{q}$ labels the points of the lattice. This gives both a position and a velocity distribution for particles perturbed by Gaussian density fluctuations with the same power spectrum as the CMB. Now we just press PLAY and watch gravity evolve the system forward. Eventually we will watch filamentary structure and galaxy clusters forming, and we will know the velocity of each particle in the system at all times.
In summary, we have examined how to compute a velocity distribution for the largest objects in our universe, those of cosmological scale.
