  The interpretation of cosmological observations relies on a notion of an
average Universe, which is usually considered as the homogeneous and isotropic
Friedmann-Lema\^itre-Robertson-Walker (FLRW) model. However, inhomogeneities
may statistically bias the observational averages with respect to FLRW, notably
for distance measurements, due to a number of effects such as gravitational
lensing and redshift perturbations. In this article, we review the main known
theoretical results on average distance measures in cosmology, based on
second-order perturbation theory, and we fill in some of their gaps. We then
comprehensively test these theoretical predictions against ray tracing in a
high-resolution dark-matter $N$-body simulation. This method allows us to
describe the effect of small-scale inhomogeneities deep into the non-linear
regime of structure formation on light propagation up to $z=10$. We find that
numerical results are in remarkably good agreement with theoretical predictions
in the limit of super-sample variance. No unexpectedly large bias originates
from very small scales, whose effect is fully encoded in the non-linear power
spectrum. Specifically, the directional average of the inverse amplification
and the source-averaged amplification are compatible with unity; the change in
area of surfaces of constant cosmic time is compatible with zero; the biases on
other distance measures, which can reach slightly less than $1\%$ at high
redshift, are well understood. As a side product, we also confront the
predictions of the recent finite-beam formalism with numerical data and find
excellent agreement.
