Speed distribution of objects in the universe If we pick a reference frame where, for example, our planet is standing still (although it should be irrelevant), what would the speed distribution of all objects in the universe look like? For the sake of simplification, we can limit ourselves to macroscopic objects (planets, stars, black holes etc.), and not consider particles (especially those rotating with great speeds, which mess with the distribution) or nebulae.
The first idea that popped in mind was the Maxwell-Boltzmann distribution, similar to gas particles, but it doesn't account for gravity (which is analogous to intermolecular forces in the kinetic theory of gasses) and all the motion resulting from it, among other things, so it's out of the game. 
Moreover, does it even make sense to talk about a speed distribution in our universe which is expanding?
Note: I found a similar question on PSE, where the author initially asked a similar thing, but finally resorted to answering a different problem. So the question is still open.
EDIT: Since it's likely that the expansion of the universe makes the problem much more complicated and in fact makes the speed distribution innormalizable, we may give a further constraint in terms of locality of the solution, i.e. limit ourselves to observing motion only within our own galaxy.
 A: There is a great deal of effort to try an answer this question in an observational way. Or let be more precise: there is lots of effort to establish what the peculiar velocity field is in our "neighbourhood". The peculiar velocity $v_p$ is the velocity that an object has in addition to the velocity we expect to see because the object is part of the "Hubble flow". On way of defining this is that it is the velocity with respect to the cosmic microwave background, which defines a local frame of reference for the co-moving volume.
i.e.:
$$ v_p = cz - H_0 d,$$
where $z$ is the redshift, $c$ is the speed of light, $d$ is the distance and $H_0$ is the Hubble parameter. 
There are several problems with this approach including: (1) Although $z$ is reasonably easy to measure, $d$ is a nightmare! (2) This only gives you the line of sight peculiar velocity. At present the tangential motion of all but the very nearest galaxies are not measurable, so one has to make assumptions and adopt some statistical reconstruction procedures to get the 3D field. (3) Uncertainties in distance translate fairly directly to uncertainties in peculiar velocity and this eventually swamps the signal. However, such measurements are important because they test ideas of cosmic structure formation - the velocity field is sensitive to the dark matter distribution for instance - and can yield estimates of cosmological parameters.
There are some spectacular visualisations of the available data in Courtois et al. (2013), including a video presentation that deals with both the positions and peculiar velocites of Galaxies within a redshift range of $\leq 5000$ km/s, corresponding to a distance of $\leq 60$ Mpc. Well worth a look!
Here is one visualisation from the paper. It shows peculiar velocities (with respect to the Hubble flow) for galaxies within about 20 Mpc. The view is confined to the 2D X,Y plane defined by our Galactic plane. Each point represents a galaxy and the axes are its position expressed in terms of a velocity (with $H_0 = 74$ km/s per Mpc). The arrows represent the amplitude of the line of sight peculiar velocity. Blue is towards us and red is away. The length of the arrow gives the amplitude of the velocity on the same scale as the axes.
You can immediately see there is loads of structure on these scales. The velocity field is most definitely not a straightforward Maxwellian or power law. The local group of galaxies is moving at around 600 km/s with respect to the CMB, there are bulk flows and various "attractors". Such flows and anisotropies appear to exist on scales at least out to 100 Mpc (e.g. Lavaux et al. 2010). Typical peculiar velocities are $\pm 300-500$ km/s.

Then here is an even more awesome view, which is a reconstructed view of the 3D velocitiy field (in the X,Y plane of our galaxy). Now each galaxy has an estimate of its 3D velocity shown and reveals the flows that are present in a more persuasive image. The paper also contains similar images in the X,Z and Y,Z planes.

Now, changing scales completely, your question also asks what the local velocity distribution of stars looks like. This is much easier to answer; it is possible to estimate the distances to stars quite precisely and their motions can be measured in the line of sight and also tangentially using proper motions. Here we can think about the problem much more "thermodynamically". To first order the velocity distributions of stars with respect to the local standard of rest, i.e. with respect to the mean motions of stars in the solar neighbourhood, can be approximated as Gaussian. However, the dispersions of these Gaussians increase with the age of a star - the so-called phenomenon of "disc heating". The absolutely classic work on this is by Wielen (1977), who quotes the following numbers (that are still roughly correct). In terms of Galactic velocity coordinate ($U$ towards the Galactic centre, $V$ in the direction of the Sun's rotation around the Galactic centre and $W$ out of the Galactic plane), the Gaussian velocity dispersions ($\sigma_U$, $\sigma_V$, $\sigma_W$) rise from around $(14,11,8)$ km/s for 1 billion year-old stars to $(34,21,21)$ km/s at 5 billion years. Roughly speaking the the total velocity dispersion increases as:
$$\sigma^2 = 100 + 6\times 10^{-7} \tau,$$
where velocities are in km/s, ages $\tau$ are in years, 
and the age distribution of stars in the solar neighbourhood is reasonably uniform between a few hundred million and 10 billion years old. There are also a small fraction (1%) of high velocity interlopers with speeds above 100 km/s from the Galactic halo.
A: 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.
A: It seems strange to me that no one mentioned the Hubble Law. 
Basically, all galaxies are receding from us with "velocity" proportional to the proper distance. 
$$v=H_0 D$$
Where $H_0$ is the Hubble constant and the "velocity" is the derivative of the proper distance with respect to cosmological time. There are some subtleties with this definition of velocity (it can be superluminal for instance).
This is the main confirmation that the universe is expanding. Earth of course has no special place in the universe, from every place in the universe you would see a similar behavior.
PS: Reading the comments to the previous question, I guess that this is not the answer that are you looking for. Anyway I think that every possible non trivial speed distribution would be strongly suppressed by these effects of universe's expansion.
A: My guess, and I hope you get a better answer, is that gravity is the only force involved, and that does not have complicated laws, unless it's a many body problem, which I don't think you are referring to here, as you want to treat cosmological objects as distinct, and far away from each other.

Moreover, does it even make sense to talk about a speed distribution in our universe which is expanding?

There is a kind of speed distribution law, in that the further the objects are away from each other, the faster their relative velocity, in general terms.

You can see the speed distribution based on Doppler shifting of spectral lines, in the diagram above.
This does not answer your question, but it is interesting to read and related to it.
Universe Computer simulation.  I have a feeling there is not such a law as you are looking for, basically because of missing data/obscured objects, but possibly extrapolations based on simulations such as the above, might resolve things
