In heavily \ensuremath{\delta}-doped semiconductor superlattices, the quantizing local space charge potential V(z,r) varies with the position r parallel to the layers, due to the statistical two-dimensional (2D) configuration of dopants. Assuming linear screening, the probability distribution of this fluctuating potential profile is analyzed using a Monte Carlo technique. As a first approximation of the resulting electronic structure, the 1D quantum problem corresponding to each sampled z profile is solved independently, assuming constant electron and hole quasi-Fermi levels throughout the sample. The subband energy fluctuations, the local wave functions, and the resulting luminescence spectra are calculated for different models. The simulated lateral fluctuations of the band edges are strongly non-Gaussian in the vicinity of the doping layers, leading to an exponential low-energy tail in the luminescence. The recombination of carriers populating tail states is drastically suppressed due to increased confinement in the model, which includes the local perturbation of the wave functions. The theoretical results are compared with electroluminescence spectra of a \ensuremath{\delta}-doped GaAs n-i-p-i superlattice. Both the expected non-Gaussian tail and the wave-function shrinkage of the tail states are confirmed by the experiments. Quantitative agreement is achieved without using fitting parameters.