Vibro-acoustic analysis of complex systems at higher frequencies faces two challenges: how to compute the response without using an excessive number of degrees of freedom (DOFs), and how to quantify the uncertainty of the response due to small spatial variations in geometry, material properties, and boundary conditions, which have a wave scattering effect? In this study, a general method of analysis is presented that provides an answer to both questions while overcoming most limitations of statistical energy analysis. The fundamental idea is to numerically compute an artificial ensemble of realizations for the components of the built-up system that are highly sensitive to small random wave scatterers. This can be efficiently performed because their eigenvalue spacings and mode shapes conform to Gaussian orthogonal ensemble spacings and Gaussian random fields, respectively. The DOFs of the overall system are therefore limited to those of the deterministic components and the interface DOFs of the random components. The method is extensively validated by application to plate structures. Good agreement between the predicted response probability distributions and the results of detailed parametric probabilistic models is obtained, also for cases of low modal overlap, single point loading, and strong subsystem coupling.