This paper presents a method for optimizing the geometry of lightning receptors of the external lightning protection system (LPS) of the high voltage air-insulated substation (AIS), which is based on the Bayesian optimization with Gaussian processes. This is a statistical optimization method, which is here additionally interfaced with a stochastic approach to the analysis of lightning exposure of AIS to direct lightning strikes. Bayesian optimization assumes that there is no closed-form expression for the objective function, only a stochastic process that uses certain inputs (disposition and geometry of lightning receptors) to produce a noise-corrupted outputs (i.e. level of shielding), while the process itself can be regarded as a black box model. Namely, the analysis of shielding of LPS is seen as an iterative stochastic process that converges to a certain criterion, rather than an exact analytical calculation. Lightning shielding analysis is based on the Monte Carlo method and the electrogeometric model. The proposed Bayesian optimization method will be applied on a concrete high voltage substation test case.