In the last decades enormous advances have been made possible for modelling complex (physical) systems by mathematical equations and computer algorithms. To deal with very long running times of such models a promising approach has been to replace them by stochastic approximations based on a few model evaluations. In this paper we focus on the often occuring case that the system modelled has two types of inputs x = (xc, xe) with xc representing control variables and xe representing environmental variables. Typically, xc needs to be optimised, whereas xe are uncontrollable but are assumed to adhere to some distribution. In this paper we use a Bayesian approach to address this problem: we specify a prior distribution on the underlying function using a Gaussian process and use Bayesian Monte Carlo to obtain the objective function by integrating out environmental variables. Furthermore, we empirically evaluate several active learning criteria that were developed for the deterministic case (...