The current contribution shows how the stochastic distribution of the buckling load of a cylindrical shell induced by random geometric imperfections can be determined with Monte Carlo simulations at low computational cost by using reduced order models. Once the buckling modes and Koiter’s b-factors of a structure have been determined, nonlinear buckling analyses can be carried out very efficiently using a reduced order model. Given a set imperfection data, the current contribution shows how these can be described as a superposition of buckling modes. Thereby, the basis to generate random realizations of geometric imperfections is the same as the basis used for the reduced order approach. This allows determining the buckling load for any random realization within a Monte Carlo simulation at extremely low cost.