Cyber-physical system models often feature stochastic behaviour that itself depends on uncertain parameters (e.g., transition rates). For these systems, verifying reachability amounts to computing a range of probabilities depending on how uncertainty is resolved. In general, this is a hard problem for which rigorous solutions suffer from the well-known curse of dimensionality. In this paper we focus on hybrid systems with random parameters whose distribution is subject to nondeterministic uncertainty. We show that for these systems the reachability probability is a smooth function of the nondeterministic parameters, and thus Gaussian processes can be used to approximate the reachability probability function itself very efficiently over its entire domain. Furthermore, we introduce a novel approach that exploits rigorous probability enclosures for training Gaussian processes. We apply our approaches to non-trivial hybrid systems case studies, and we empirically demonstrate their advantages with respect to standard statistical model checking.