  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 $S_i^z S_{i+r}^z$ 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
$S_i^zS_{i+r}^z$ correlator exhibits striking differences between the MBL and
Anderson insulator cases.
