We investigate the error properties of certain galaxy luminosity function (GLF) estimators. Using a cluster expansion of the density field, we show how, for both volume- and flux-limited samples, the GLF estimates are covariant. The covariance matrix can be decomposed into three pieces: a diagonal term arising from Poisson noise, a sample variance term arising from large-scale structure in the survey volume and an occupancy covariance term arising due to galaxies of different luminosities inhabiting the same cluster. To evaluate the theory one needs the mass function and bias of clusters, and the conditional luminosity function (CLF). We use a semi-analytic model (SAM) galaxy catalogue from the Millennium Run N-body simulation and the CLF of Yang et al. to explore these effects. The GLF estimates from the SAM and the CLF qualitatively reproduce results from the two degree Field Galaxy Redshift Survey (2dFGRS). We also measure the luminosity dependence of clustering in the SAM and find reasonable agreement with 2dFGRS results for bright galaxies. However, for fainter galaxies, L < L*, the SAM overpredicts the relative bias by ˜10-20 per cent. We use the SAM data to estimate the errors in the GLF estimates for a volume-limited survey of volume V ˜ 0.13 h-3 Gpc3. We find that different luminosity bins are highly correlated: for L < L* the correlation coefficient is r > 0.5. Our theory is in good agreement with these measurements. These strong correlations can be attributed to sample variance. For a flux-limited survey of similar volume, the estimates are only slightly less correlated. We explore the importance of these effects for GLF model parameter estimation. We show that neglecting to take into account the bin-to-bin covariances, induced by the large-scale structures in the survey, can lead to significant systematic errors in best-fitting parameters. For Schechter function fits, the most strongly affected parameter is the characteristic luminosity L*, which can be significantly underestimated.