Recent work has shown that modified gravitational wave (GW) propagation can be a powerful probe of dark energy and modified gravity, specific to GW observations. We use the technique of Gaussian processes, that allows the reconstruction of a function from the data without assuming any parametrization, to measurements of the GW luminosity distance from simulated joint GW-GRB detections, combined with measurements of the electromagnetic luminosity distance by simulated DES data. For the GW events we consider both a second-generation LIGO/Virgo/Kagra (HVLKI) network, and a third-generation detector such as the Einstein Telescope. We find that the HVLKI network at target sensitivity, with O(15) neutron star binaries with electromagnetic counterpart, could already detect deviations from GR at a level predicted by some modified gravity models, and a third-generation detector such as ET would have a remarkable discovery potential. We discuss the complementarity of the Gaussian processes technique to the (Ξ0, n) parametrization of modified GW