We present a novel approach to compute reachable sets of dynamical systems with uncertain initial conditions or parameters, leveraging state-of-the-art statistical techniques. From a small set of samples of the true reachable function of the system, expressed as a function of initial conditions or parameters, we emulate such function using a Bayesian method based on Gaussian Processes. Uncertainty in the reconstruction is reflected in confidence bounds which, when combined with template polyhedra ad optimised, allow us to bound the reachable set with a given statistical confidence. We show how this method works straightforwardly also to do reachability computations for uncertain stochastic models.