Load-bearing elements composed of glass, as known, are often susceptible to buckling collapse mechanisms. This intrinsic characteristic (and thus potential limitation for design) typically derives from the use of relatively small thicknesses to cover large spans and surfaces, thus resulting in a multitude of columns, beams, or plates that are characterized by high slenderness. In the literature, accordingly, several and design propositions support of the definition of efficient calculation models to capture the typical buckling response of glass elements variably shaped, sized, restrained and loaded. In this study, the attention is focused on the buckling analysis of glass columns and on the assessment of uncertainties due to input random parameters. With the support of finite element numerical models, a total of 800 glass columns are investigated, by accounting for stochastic variations in the geometry (size and thickness), modulus of elasticity and density of glass, maximum amplitude of the imposed initial imperfection, material type. Based on the Monte Carlo simulation method, the final result takes the form of 2400 simulations, where the postprocessing analysis is spent on the derivation of empirical formulations for the correlation of the relevant buckling capacity indicators. From the global out-of-plane bending analysis, the input random parameters are observed to affect severely both long and short columns with different flexural stiffness. Besides, a stable linear correlation is found for some influencing indicators. The attention is thus focused on the sensitivity analysis of critical buckling load, ultimate failure configuration, deflection at collapse, buckling reduction coefficient.