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 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 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 SziSzi+r correlators skewed towards negative values; and a multimodal distribution for spin-flip correlators. A quantitative quasi-degenerate perturbation theory calculation of these correlators yields a surprising 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 SziSzi+r correlator exhibits striking differences between the MBL and Anderson insulator cases.