Hierarchical Bayesian inference aims to infer the population-wide hyperparameters that characterise the joint distribution of input parameters of a computational model M. Experimental observations y are used to condition population priors through a likelihood function that depends on M, typically requiring repeated evaluations of the forward model for multiple sets of (hyper-)parameters. This quickly becomes computationally prohibitive for high-dimensional parameter spaces or computationally expensive forward models. A common strategy to alleviate this issue is to replace the forward model with a deterministic surrogate, such as polynomial chaos expansions. However, this approach requires an accurate approximation of the forward model across the full input space, which is often intractable in high dimensions. This work presents an alternative approach based on stochastic emulators, which eliminates the need for accurate surrogates over the entire input space. Rather than approximating the forward model itself, we directly emulate the stochastic population response as a function of the hyperparameters of interest. The proposed framework implicitly marginalises over the individual-level parameters and allows direct evaluation of the marginal likelihood of the hyperparameters. This significantly reduces the dimensionality of the problem and enables efficient hierarchical Bayesian inference in settings that are otherwise computationally infeasible, particularly when the population-level parameters are targeted or where the underlying forward model is inherently stochastic.