An important part of cosmological model fitting relies on correlating distance indicators of objects (for example type Ia supernovae) with their redshift, often illustrated on a Hubble diagram. Comparing the observed correlation with a homogeneous model is one of the key pieces of evidence for dark energy. The presence of cosmic structures introduces a bias and scatter, mainly due to gravitational lensing and peculiar velocities, but also due to smaller non-linear relativistic contributions which are more difficult to account for. For the first time we perform ray tracing onto halos in a relativistic N-body simulation, taking into account all the relevant effects from general relativity, both in the evolution of structure and in the projection of that structure onto the observer's past light cone. We show that the mean of the bias in the Hubble diagram is indeed as small as expected from perturbation theory, but the distribution of sources is significantly skewed with a very long tail of highly magnified objects. We illustrate how important this is in the context of supernovae forecasts, where we show that a flat background in the true cosmology would be ruled out at high significance when a Gaussian likelihood in the Hubble diagram is assumed. We show that this bias cannot be reduced by increasing the number of sources and instead requires a proper modeling of non-linear effects.