Geometrical imperfections are ubiquitous in load-bearing structures, including beams, columns, and shells. Fabrication processes of structural members most often create geometrical imperfections of random size and shape, which lead to non-deterministic load-carrying capacity. This study investigates the statistics of the buckling load of a beam with a random initial imperfection profile that rests on a nonlinear elastic foundation. The geometrical imperfection is represented by a zero-mean Gaussian random field, generated using the Karhunen–Loève expansion. The spatial distribution of the random imperfection is characterized by the probability distribution of the local imperfection magnitude and a spatial autocorrelation function. A finite-difference scheme is used to solve the governing equilibrium equation for a given initial imperfection profile, from which the buckling load is determined. Through a set of Monte Carlo simulations, the mean and variance of the buckling load are determined. The simulations reveal the influence of different length scales on the statistics of the buckling load, including the beam length and the autocorrelation length of the geometrical imperfection. The size effects predicted with the simplified model have implications for reliability-based structural design.