Ergodic quantum many-body systems satisfy the eigenstate thermalization hypothesis (ETH). However, strong disorder can destroy ergodicity through many-body localization (MBL) -- at least in one dimensional systems -- leading to a clear signal of the MBL transition in the probability distributions of energy eigenstate expectation values of local operators. For a paradigmatic model of MBL, namely the random-field Heisenberg spin chain, we consider the full probability distribution of \emph{eigenstate correlation functions} across the entire phase diagram. We find gaussian distributions at weak disorder, as predicted by pure ETH. At intermediate disorder -- in the thermal phase -- we find further evidence for \emph{anomalous thermalization} in the form of heavy tails of the distributions. In the MBL phase, we observe peculiar features of the correlator distributions: a strong asymmetry in \(S_i^z S_{i+r}^z\) correlators skewed towards negative values; and a multimodal distribution for spin-flip correlators. A \emph{quantitative quasi-degenerate perturbation theory} calculation of these correlators yields a surprising \emph{agreement of the full distribution} with the exact results, revealing, in particular, the origin of the multiple peaks in the spin-flip correlator distribution as arising from the resonant and off-resonant admixture of spin configurations. The distribution of the \(S_i^zS_{i+r}^z\) correlator exhibits striking differences between the MBL and Anderson insulator cases.