  We present results of the investigations of the statistical properties of a
joint density and velocity divergence probability distribution function (PDF)
in the mildly non-linear regime. For that purpose we use both perturbation
theory results, extended here for a top-hat filter, and numerical simulations.
  In particular we derive the quantitative (complete as possible up to third
order terms) and qualitative predictions for constrained averages and
constrained dispersions -- which describe the nonlinearities and the
stochasticity properties beyond the linear regime -- and compare them against
numerical simulations. We find overall a good agreement for constrained
averages; however, the agreement for constrained dispersions is only
qualitative. Scaling relations for the Omega-dependence of these quantities are
satisfactory reproduced.
  Guided by our analytical and numerical results, we finally construct a robust
phenomenological description of the joint PDF in a closed analytic form. The
good agreement of our formula with results of N-body simulations for a number
of cosmological parameters provides a sound validation of the presented
approach.
  Our results provide a basis for a potentially powerful tool with which it is
possible to analyze galaxy survey data in order to test the gravitational
instability paradigm beyond the linear regime and put useful constraints on
cosmological parameters. In particular we show how the nonlinearity in the
density--velocity relation can be used to break the so-called Omega-bias
degeneracy in cosmic density-velocity comparisons.
