Correlation functions are an omnipresent tool in astrophysics, and they are routinely used to study phenomena as diverse as the large-scale structure of the Universe, time-dependent pulsar signals, and the cosmic microwave background. In many cases, measured correlation functions are analyzed in the framework of Bayesian statistics, which requires knowledge about the likelihood of the data. In the case of correlation functions, this probability distribution is usually approximated as a multivariate Gaussian, which is not necessarily good approximation -- hence, this work aims at finding a better description. To this end, we exploit fundamental mathematical constraints on correlation functions, which we use to construct a quasi-Gaussian likelihood. We explain how to compute the constraints, in particular for multi-dimensional random fields, where this can only be done numerically, check the quality of the quasi-Gaussian approximation, and compare it to alternative approaches -- most importantly, we test the new-found description of the likelihood in a toy-model Bayesian analysis. Finally, we compute correlation functions from the Millennium Simulation and show that they obey the constraints. By studying statistical properties of the measured correlation functions, we present further indications for the validity of the quasi-Gaussian approach
