Predict the distribution of the residuals. This will involve modelling a contribution from noise within the data and a contribution from uncertainty in the neural network weights. You can model the former with negative log-likelihood loss. And you can model the latter either by training multiple of the same model and taking the average and standard deviation of their predictions, or by using Monte Carlo dropout.

For more info, read this nice paper. This one is good too, and there's examples of code as well.