  Structure formation in our Universe creates non-Gaussian random fields that
will soon be observed over almost the entire sky by the Euclid satellite, the
Vera-Rubin observatory, and the Square Kilometre Array. An unsolved problem is
how to analyze best such non-Gaussian fields, e.g. to infer the physical laws
that created them. This problem could be solved if a parametric non-Gaussian
sampling distribution for such fields were known, as this distribution could
serve as likelihood during inference. We therefore create a sampling
distribution for non-Gaussian random fields. Our approach is capable of
handling strong non-Gaussianity, while perturbative approaches such as the
Edgeworth expansion cannot. To imitate cosmological structure formation, we
enforce our fields to be (i) statistically isotropic, (ii) statistically
homogeneous, and (iii) statistically independent at large distances. We
generate such fields via a Monte Carlo Markov Chain technique and find that
even strong non-Gaussianity is not necessarily visible to the human eye. We
also find that sampled marginals for pixel pairs have an almost generic
Gauss-like appearance, even if the joint distribution of all pixels is markedly
non-Gaussian. This apparent Gaussianity is a consequence of the high
dimensionality of random fields. We conclude that vast amounts of non-Gaussian
information can be hidden in random fields that appear nearly Gaussian in
simple tests, and that it would be short-sighted not to try and extract it.
