  As a step towards a more accurate modelling of redshift-space distortions in
galaxy surveys, we develop a general description of the probability
distribution function of galaxy pairwise velocities within the framework of the
so-called streaming model. For a given galaxy separation $\vec{r}$, such
function can be described as a superposition of virtually infinite local
distributions. We characterize these in terms of their moments and then
consider the specific case in which they are Gaussian functions, each with its
own mean $\mu$ and dispersion $\sigma$. Based on physical considerations, we
make the further crucial assumption that these two parameters are in turn
distributed according to a bivariate Gaussian, with its own mean and covariance
matrix. Tests using numerical simulations explicitly show that with this
compact description one can correctly model redshift-space distorsions on all
scales, fully capturing the overall linear and nonlinear dynamics of the galaxy
flow at different separations. In particular, we naturally obtain
Gaussian/exponential, skewed/unskewed distribution functions, depending on
separation as observed in simulations and data. Also, the recently proposed
single-Gaussian description of redshift-space distortions is included in this
model as a limiting case, when the bivariate Gaussian is collapsed to a
two-dimensional Dirac delta function. We also show how this description
naturally allows for the Taylor expansion of $1+\xi_S(\vec{s})$ around
$1+\xi_R(r)$, which leads to the Kaiser linear formula when truncated to second
order, expliciting its connection with the moments of the velocity distribution
functions. More work is needed, but these results indicate a very promising
path to make definitive progress in our program to improve RSD estimators.
