MCMC generation of cosmological fields far beyond Gaussianity

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.


INTRODUCTION
Random fields are ubiquitous in cosmological research, as they permeate the Universe ever since its earliest phases. In fact, the theory of cosmological inflation reasons that the initial distribution of matter in our Universe traces back to quantum fluctuations, and our Universe is thus inherently random from the outset (Mukhanov et al. 1992;Durrer 1994).
At the cosmic time of recombination, a snapshot of the statistical state of the Universe's matter fields was generated, which has been observed as the cosmic microwave background (CMB) by a series of experiments (Chiang et al. 2010;Hinshaw et al. 2013;Sievers et al. 2013;Planck Collaboration 2018). These data reveal that the initial state of cosmic fields is compatible with Gaussian distributions whose power spectra are by now well measured.
Primordially generated non-Gaussianity is theoretically expected to be weak, but a detection of it would further constrain the physics of inflation. Accordingly, there exists a vast body of literature focusing on weakly non-Gaussian fields and their analysis, e.g. Komatsu (2010); Maldacena (2003); Planck Collaboration (2019).
In the post-CMB cosmic evolution, physical processes altered the statistical properties of the random cosmic fields, thereby creating structures which are richer than those of a Gaussian random field. Amongst these processes ranges gravity which -unlike primordial mechanisms -produces strong non-Gaussianity, thereby immediately hindering a direct transfer of CMB analysis sellentin@strw.leidenuniv.nl methodology to gravitationally evolved fields. In this paper, we hence focus on random fields which are strongly non-Gaussian.
Besides gravity, other mechanisms generating strong non-Gaussianity include also the radiation of the first stars which burn bubble-like patterns into the cosmic hydrogen distribution. Non-Gaussian fields therefore also occur in studies of cosmic reionization (Majumdar et al. 2018;Watkinson et al. 2019). Current observations of non-Gaussian fields include e.g. galaxy clustering and cosmic shear experiments (Abbott et al. 2019;Hildebrandt et al. 2020). In the near future, the Euclid satellite (Laureijs et al. 2011), the Vera-Rubin observatory (LSST Science Collaboration et al. 2009), and the Square Kilometre Array (SKA) (Mellema et al. 2013) will supersede contemporary observations both in sky-coverage, depth, and data quality.
With these new observational capabilities, the window on non-Gaussian fields is opened. It is thus timely to focus on the statistical challenge of analyzing strongly non-Gaussian fields. This study complements a significant body of literature handling the classical differential equations that encode the physical laws which generate the non-Gaussianity. Such deterministic models of physics include full N-body simulations (e.g. Springel et al. (2018)), and approximations including Zel'Dovich (1970) and the second order Lagrangian perturbation theory (e.g. described in Scoccimarro (1998)). These approaches model the deterministic growth of non-Gaussian structures from random initial conditions. Once these deterministic codes output a final non-Gaussian field, a statistical solution is required to compare the modelled random fields to observed fields. There currently exists no consensus on how this comparison is achieved best, and in this paper we study it from a Bayesian point of view.
Targeting the inference of parameters from a non-Gaussian field in the long run, we begin by generating non-Gaussian fields from a flexible non-Gaussian sampling distribution. The reason for doing so, is that the sampling distribution becomes the likelihood when analyzing an observed field. We clarify this setup in Sect. 2, and present our non-Gaussian ansatz, and fields generated from it, in sections 3-7.

INFERENCE FROM NON-GAUSSIAN FIELDS
Increasingly accurate all-sky surveys confront contemporary cosmology with the challenge of inferring parameters from non-Guassian fields. In this section, we compare the problem to the more familiar case of inferring parameters from the Gaussian CMB, and we describe the challenges posed by the update to non-Gaussianity.
We imagine the initial data set is a pixelized random field, i.e. each pixel provides one data point and all pixels could be stored in a data vector. If ∆ i denotes the measured field value in the ith pixel, then the full data vector comprising all pixels will be ∆. In case of the CMB one typically works with the spherical harmonic transform of ∆, which are the usual modes a m . We collect all a m -modes in a data vector a m . For fixed , the variance of the a m -modes will be the spherical harmonic power spectrum C . The covariance matrix of a m will be diagonal for a full-sky observation, and we denote this covariance matrix by F = diag(C min , ..., C max ). The unconventional choice of naming the CMBs covariance matrix F provides a slightly smoother transition to the upcoming non-Gaussianity study.
Inferring physical parameters from the CMB implies the posterior probability of parameters θ in light of the data a m must be computed. We shall denote this posterior distribution as P(θ|a m ).
Inference begins by specifying a prior π(θ) on the physical parameters. In case of the CMB, drawing values for the parameters from the prior, the power spectra C (θ) can be computed, implying in Bayesian jargon that the matrix F is 'given' once the parameters θ are fixed.
Due to the CMB being a Gaussian random field, the a m -modes are drawn from a Gaussian sampling distribution of covariance matrix F. The likelihood must therefore be the Gaussian distribution G, and the posterior simplifies to the familiar P(θ|a m ) ∝ G a m |F(θ) π(θ). (1) We now extend to inference with a non-Gaussian field. In this paper, we shall work in real space: instead of using a m we use the pixel values ∆. For a non-Gaussian field, there will exist a matrix F as before, but we introduce two additional tensors S and Q that quantify the random field's non-Gaussianity, in a sense that will be described in due course. Inference with a non-Gaussian field will again need to begin by specifying a prior π(θ), but now all the tensors F(θ), S(θ), Q(θ) must be computed once θ is given. The sampling distribution for non-Gaussian fields will then be P(∆|F, S, Q) and the posterior is then P(θ|∆) ∝ P ∆|F(θ), S(θ), Q(θ) π(θ). (2) The probability distribution P(∆|F, S, Q) is the non-Gaussian sampling distribution of the fields, and the objective of this paper is hence to understand this probability distribution. The natural use of such a non-Gaussian distribution would be in a Monte-Carlo based cosmological inference, that ingests cosmic fields directly, see for example Wei & Murray (2016); Taylor et al. (2008); Alsing et al. (2017). Providing a good ansatz for this non-Gaussian distribution is an intellectual challenge in itself, and precedes any numerical challenges. We thus devote Sect. 3 to establishing a good ansatz, and in Sect. 4 we specialize our non-Gaussian distribution to statistically isotropic and homogeneous fields as in cosmology. In Sect. 4.1 we discuss the numerical challenges of non-Gaussian field generation arising from the many permutations in non-Gaussian calculations and the high dimensionality of the problem. Example fields are then presented in Sect. 6.

DISTRIBUTION
In this section we describe our ansatz for a probability distribution from which non-Gaussian random fields can be drawn.
Describing non-Gaussianity in a simple mathematical form is challenging: per definition it is anything else but Gaussian, but no further information is provided. The standard approach to handling non-Gaussian fields uses higher-order moments or cumulants, which in the case of studying random fields corresponds to computing 3-point and 4-point correlation functions, or even higher n-point correlation functions. Such cumulants can be used to construct the Edgeworth expansion, but this expansion neither yields positive-definite probability distributions, nor does it converge (Sellentin et al. 2017). The hierarchy of correlation functions is accordingly known to capture the information of non-Gaussian fields only incompletely and inefficiently (Carron & Szapudi 2017).
To overcome the problem of negative probabilities and the lack of convergence, Sellentin et al. (2014); Sellentin (2015) (S14 and S15 from now on) therefore developed a novel series expansion of non-Gaussian distributions, called DALI. In comparison to the Edgeworth expansion, DALI provides a guaranteed positive distribution at all expansion orders. It inherits its convergence properties from a Taylor series, and has been shown to converge extremely rapidly: in practice, DALI reproduces extremely non-Gaussian distributions already after one or two terms beyond Gaussianity.
DALI, as published in S14 and S15, was developed as an extension of a Fisher matrix formalism. As a Fisher matrix approach does not include 'real' data, we therefore begin by strongly modifying the original DALI setup: a data vector has to be introduced. The outcome of this reformulation will be that DALI is then able to analyze real data, or generate statistically compatible 'fake' data.

Review: DALI as a Fisher Matrix extension
Here we briefly repeat the main points in the derivation of DALI as presented in S14 and S15 and derive the sampling distribution used in this paper.
If d is the dimension of data space, and if θ = (θ 1 , ..., θ p ) is a p-dimensional parameter vector to be inferred, then the original DALI-formulation provides an analytic approximation for the expected non-Gaussian posterior of parameters, given Gaussian data with mean µ(θ) and data covariance matrix C. Hence, DALI approximates the posterior P[θ|µ(θ)] meaning as in any Fisher forecasting framework, also the non-Gaussian DALI extension replaces not-yet taken data by their expectation value µ(θ), evaluated at fiducial parameter valuesθ (Tegmark et al. 1997;Heavens et al. 2014;Amendola & Sellentin 2016).
Denoting partial derivatives with respect to θ α by commas, ∂ α µ = µ, α , the best fitting parameters byθ, and the offset of the αth parameter from its best fitting value as ∆ α = θ α −θ α , we can Taylor expand the mean as where Einstein's summation convention over repeated indices is implied. The DALI expansion of the posterior then follows to be Here the factors N ! were not cancelled against the numerical prefactors in angular brackets, as the former are a reminder of the Taylor expansion from which Eq. (4) originates, whereas the numerical prefactors arise from symmetries due to permuting derivatives and the symmetry of the scalar product. For details of this derivation, please refer to S14 and S15. In Eq. (4), we truncated the DALI expansion at second-order derivatives of the mean. For higher-order derivatives, further non-Gaussian distributions can be approximated, but already at second-order most singlepeaked distributions are well approximated (Sellentin & Schäfer 2016). In this paper, we thus focus on secondorder expansions. The components of the tensors F, S, Q then read We note that the tensors S and Q both contain second derivatives of µ as inputs. They hence appear simultaneously in the expansion, and satisfy relations with respect to each other. This guarantees the positive definiteness of the DALI-expansion.

DALI FIELD GENERATION
In this paper we aim to use DALI to generate random fields under the constraints imposed by cosmology. In this section, we first modify the DALI distribution to meet this aim. We then introduce our pixelization scheme of the random field, and then present the implementation of the constraints arising from statistical isotropy and homogeneity due to the cosmological principle. We further impose independence at large distances due to causality and the finite speed of light, and domination of Gaussianity at large distances, as demanded by cosmology's theory of inflation, and the CMB.
We begin by modifying the DALI-distribution such that it can serve as sampling distribution for non-Gaussian fields. In other words, we want to draw ∆ from the DALI distribution ∆ ∼ P(∆|F, S, Q).
In this case, the distribution Eq. (4) can be simplified: The matrix C −1 in Eq. (4) plays the role of a metric, and in order for Eq. (4) to be a positive definite function, the matrix C −1 has to be positive definite. It can therefore be Cholesky decomposed as where L is a unique lower triangular matrix, sometimes referred to as the 'root'. Likewise, L is upper triangular. This allows us to introduce the vectors Eq. (4) can thus be rewritten as We introduce three types of angles, ϑ, η, ζ such that where w i , W ij are the absolute values of w i , W ij . To satisfy identity constraints, we set The remaining values, w i , W ij , cos(θ ij ), cos(η ij,k ), cos(ζ ij,kl ) are free parameters that have to be set to satisfy physical constraints. For example, if all W ij were zero, a Gaussian random field would ensue, whose covariance matrix is defined by all w i and cos(θ ij ).
As the indices i, j, k, l run over the pixels in the field, and as physics depends on distance, we now first have to introduce the pixelization scheme, before the values for w i , W ij , cos(θ ij ), cos(η ij,k ), cos(ζ ij,kl ) can be further specified.

Pixelization scheme
The fields generated in this paper will be square and have N pixels. We enumerate the pixels with roman  indices i ∈ [1, N ]. If x(i) denotes the (Cartesian) xcoordinate of the ith pixel, and y(i) its y-coordinate, then the distance r between two pixels i and j is given by The ordering of the pixels is irrelevant, all physics will depend only on r. In extension, the distance between a triplet or quadruplet of pixels is given by We will denote the vector of all pixel values as ∆, the sidelength of the field f = √ N . This means we will have N random variables in ∆. Hence, field-generation quickly becomes a very high dimensional problem.

Statistical Isotropy and Homogeneity
On sufficiently large scales, standard cosmology assumes the Universe to be statistically homogeneous and isotropic. Statistical homogeneity is given by For future reference we note that the conditional distribution of density is no longer a simple function of distance on small scales. This implies that two point statistics are no longer sufficient to characterise a random field on such scales. In this paper we focus on the largest cosmic scales, but our framework could be applied to smaller scales by modifying the tensors F, S and Q to make them dependent on more than only distance on small scales. An emergence of filamentary structures is hence not expected for the results here shown, but can (in future) be provided for.

Independence and causality
To ensure that different pixels of the field are statistically independent at large distances we impose the limit which ensures that the covariance of widely separated pixels approaches zero. Furthermore, to ensure that pixels which are far apart are causally disconnected we impose the following three constraints

MCMC GENERATION OF NON-GAUSSIAN FIELDS
In this section we generate non-Gaussian random fields from DALI with a Monte Carlo Markov Chain (MCMC) approach.
As described in Sect. 4.1, we imagine the random field to be pixelated, as would occur when taking data. Each pixel ∆ i is then a random variable, and the collection of all pixels yields the random vector ∆. The joint distribution of all pixel values is given by D(∆|F, S, Q) and defines which random patterns the field will form.
A realization of one random field is then yielded by drawing from the distribution ∆ ∼ D(∆|F, S, Q).
We therefore must succeed in creating an algorithm that indeed generates samples drawn as defined by Eq. (17). The crux to succeeding will hereby be the high dimensionality of the distribution: for a square pixelized field of side length f , the number of random pixels to be jointly drawn will be f 2 , therefore very rapidly encountering the 'Curse of Dimensionality'. We therefore generate our random fields with a Hamiltonian Monte Carlo (HMC) sampler (Betancourt 2017). Intuitively, a Hamiltonian Monte Carlo sampler can be understood as a clever way of increasing the distance between samples of a more conventional Metropolis-Hastings sampler. This is achieved by integrating equations of motion to find a sample at similar probability. Readers not interested in the details of the HMC sampler may skip directly to the resulting fields presented in Sect. 6.
The Hamiltonian Monte Carlo algorithm begins by the observation that if a probability distribution is turned 'upside down', then the former maximum will become the lowest point of a valley-like potential. The sampler then receives random kicks at selected points in time, and is left to propagate through the potential between the kicks. It thereby explores the probability distribution much more efficiently than a sampler with a Gaussian proposal, which has no guidance, centres on the old point and is generally very narrow and hence slowly step-bystep explore its local neighbourhood.
Implementing the HMC algorithm, we define the (negative) potential L(∆) to be We introduce momenta p and a mass matrix M. The mass matrix is a tunable parameter that determines the efficiency of the sampler. The momenta p are auxiliary variables that are not of primary interest: they are introduced to upgrade the sampler to an evolving Hamiltonian system. The combination (∆, p) is accordingly referred to as 'phase space'. We define the Hamiltonian the augmented probability distribution is proportional to e −H , implying where we have ignored the normalisation coming from the mass matrix term as it is irrelevant for this work. As the original distribution D factors out, the probability distribution of ∆ is the marginal of e −H over the momenta. This underlines that the momenta are indeed only auxiliary variables.
The HMC algorithm begins by selecting a starting point ∆ 0 , and drawing an initial momentum vector The position and momentum are then updated by solving the Hamiltonian system 'Solving' here implies integrating these differential equations to track how the sampler propagates through phase space. Numerically solving this system is prone to numerical inaccuracies, which will affect the acceptance rate of the sampler and might also break its reversibility. We hence choose the Leapfrog algorithm to solve the equations of motion Eq. (22), which is a symplectic integrator conserving energy to second order yielding a trajectory with an energy that fluctuates around the true energy of the exact trajectory. The leapfrog algorithm is given by where is a small stepsize appearing due to having discretized time. The gradient and Hessian of a DALI distribution are given in App. A. The three steps of Eq. (23) are repeated N times for N updates of position. As i ∈ [1, f 2 ] denotes the value ∆ i of each pixel, the leapfrog loops runs in parallel (and synchronized) for all pixels.
As the leapfrog algorithm's energy oscillates around the true energy of the integrated trajectory, stopping the integration at any point in time is likely to yield an energy somewhat different from that of an exact integration. We thus decide upon whether the end point (∆ E , p E ) is accepted as a valid sample via the acceptance criterion known from the Metropolis-Hastings algorithm.
If the point (∆ E , p E ) is accepted, it takes the place of ∆ 0 , otherwise ∆ 0 remains unchanged. Having completed these steps, the algorithm reruns from Eq. (21) thereby building up a chain of samples.

RESULTS
In this section we present results from a non-Gaussian field generated using the DALI distribution and an HMC sampler. We will use this non-Gaussian field to show the feasibility of generating such fields using DALI and an HMC sampler, and to demonstrate the subtleties involved in analysing them.
The non-Gaussian field used in this section only has a Q tensor (see equation 4) using w i = 0, cos(η ij,k ) = 0, W ij = 1.5 r ij , R = r 2 13 + r 2 14 + r 2 23 + r 2 24 , We will therefore call it a Q-field. Figure 3 shows the distance dependence of these functions. Their decay towards zero at large pixel separations implies pixel far apart will be statistically independent. Homogeneity and isotropy are ensured by using functions of distances only. We generate fields with side length f = 70, which gives a total dimensionality of 4900. Sampling such non-Gaussian field is extremely expensive, and we therefore switch off non-Gaussianity beyond 1/8 th of the field's side length. This approximation is motivated by physical processes acting locally, with e.g. the electromagnetic and gravitational force diminishing as a function of distance and becoming causally disconnected at infitite distance.
Though our sampling is run using C++ and makes use of parallelization on 128 CPUs as well as 128GB of memory, we are limited to a 70x70 pixel field. This is largely due to the many dot products in Eq. 9 which require retrieving stored pixel values a large number of times. We have extensively optimized the algorithm to make use of the isotropy and homogeneity of the field, both in storing the W ij W kl values, the ordering of the for-loops and the combination of identical pixel combinations.  One aim is to compare this non-Gaussian Q-field with a Gaussian field of equal covariance matrix. To this end, after generating the non-Gaussian field, we compute the sample covariance matrix from the MCMC chain and use it in the HMC sampler to generate a Gaussian field. The thus created Gaussian field can be interpreted as the closest Gaussian approximation to our non-Gaussian field.
As sampling in 4900 dimensions is difficult, we iteratively improved the quality of our chains by plotting the HMC phase space trajectory in all dimensions. We equally plotted the trajectory in pixel-space only, to judge which step sizes are adequate. With suitable step sizes found, we computed the correlation lengths (Goodman & Weare 2010) and the Gelman-Rubin convergence diagnostics (Gelman & Rubin 1992). These are listed in Table 1 with NG abbreviating 'non-Gaussian', 'G' abbreviating 'Gaussian', and followed by the number of chains from which the convergence diagnostics were computed. We conclude our chains decorrelate after 3 to 30 steps, and indeed converge to the target distribution. Figure 4 shows seven consecutive samples from an arbitrary point in the chain for both the Gaussian and non-Gaussian field. From these images, it would be near impossible to tell them apart by eye. It illustrates that non-Gaussian fields can look deceivingly similar to Gaussian fields, even when their sampling distributions differ strongly.
The extreme non-Gaussianity of the Q-field's sampling distribution is illustrated by Figure 5. The top row shows the conditionals of the non-Gaussian field, where the bottom row gives the marginals. The conditionals are clearly non-Gaussian. Figure 5 also shows that for small distances (left in the figure) between two pixels the non-Gaussian conditional cannot be factorized into to independent distributions for the two pixels. This implies close pixels are statistically co-dependent. Towards the right of Figure 5 the distance between the studied pixels increases, and their joint probability distribution approaches a factorizable distribution. This is a visual impression of pixels at large distances becoming increasingly independent of each other. Figure 5 also illustrates that the strong non-Gaussianity of the sampling distribution is essentially invisible when plotting marginals of pixel pairs: the marginals in the lower row are all similar to an ellipsoidally contoured distribution. The washing out of the non-Gaussian shape of the sampling distribution is due to the high dimensionality of the sampling problem: the plotted marginals are integrals over 4898 dimensions and evidently the projection over many thousand dimensions will cause a loss of structure in the remaining two dimensions. Likewise, the two-dimensional conditionals in the top row only show such clear non-Gaussianity because all other dimensions are held fixed at zero when generating the conditional. When sampling, the probability of having a high dimensional system in such a state that all dimensions except two are exceedingly close to zero is negligible.
The appearance of shapeless marginals is particularly interesting in the context of how non-Gaussian fields are usually approached: often, histograms of pixel values are computed, and if these do not show a strong deviation from a Gaussian, then the field is judged as 'close to Gaussian'. The Q-field here studied shows however that significant non-Gaussian information can be hidden in a field whose marginals look almost Gaussian. Accordingly, non-Gaussian analyses of cosmological fields should be attempted even when 1-point and 2-point analyses do not evidence any striking departure from Gaussianity.
For example, we found a better discriminator between Gaussianity and non-Gaussianity to be the calculation of cumulants. We present measurements of the skewness and kurtosis in Figure 6. The skewness is defined as Due to the use of the third power, the skewness measures asymmetry in a probability distribution. The Q-field has a symmetric distribution, due to the lack of a symmetry breaking S tensor. Accordingly, the skewness must be compatible with zero in our case, and non-Gaussianity is detected in the excess kurtosis, defined as The skew and kurtosis are computed for each sampled field, with N the number of pixels in the field. The resulting histograms for the skewness and curtosis are shown in Figure 6. As expected, neither the Gaussian field, nor the Q-field show statistically significant departures of the Fig. 4.-Seven consecutive samples generated for a non-Gaussian Q-field (top) and its closest approximating Gaussian field which has the same covariance (bottom). This figure illustrates that the human eye struggles tell apart the Gaussian approximation and the original non-Gaussian field. Nonetheless, the non-Gaussian information is detected at high statistical significance in an analysis. Going from left to right, the conditional joint distribution morphs from a star to a box, which is caused by it beginning to factorize for increasing distance. This transition was implemented to enforce causality, such that close pixels are co-dependent, while pixels with increasing separation become ever more independent (see Sect. 4.3). The upper row (conditionals) shows clear non-Gaussian co-dependency between two pixels. The lower row imitates an an analysis of observed pixelized fields, where histograms of the pixel values are taken. It evidences that taking histograms of pixel values (marginalizing) easily triggers the false impression that the field were close to Gaussian. skewness from zero. The kurtosis measurements show a clear difference between the Gaussian and non-Gaussian fields. While the kurtosis of the Gaussian fields is consistent with zero, the non-Gaussian fields show a positive offset with an average kurtosis of about 0.2.
We further note that all measurements of n-point functions where n < 4 show no disagreement between the Gaussian and non-Gaussian field. We verify isotropy and homogeneity of the field in App. B.
Finally, we study how an initially Gaussian field evolves randomly into a non-Gaussian field. To this end, we initialize the HMC sampler with a Gaussian field of mean zero and the same covariance matrix as the Q-field.
We then let the sampler run in the potential well of the Q-field, such that the burn-in phase of the sampler causes a transition from Gaussianity to non-Gaussianity. Keeping the colour bar of the plotting range fixed then results in Figure 7: We here show a 60x60 field, different from the 70x70 field shown before, that starts from a Gaussian and remains a Gaussian in the bottom row. Accordingly, the bottom row simply depicts consecutive samples of a Gaussian field from left to right. In contrast, the non-Gaussian Q-field is plotted in the top row, where the order left-to-right again shows consecutive samples: here we see that the initially Gaussian field transforms into a non-Gaussian field as a function of time. The emergence of more pronounced structures than in the Gaussian field can be seen by comparing the top and bottom row. Both fields start from the exact same initial configuration, but the non-Gaussian field can be seen to quickly deviate from this and form much stronger fluctuations in pixel values.

CONCLUSIONS
Non-Gaussian random fields are ubiquitous in observational cosmology. A commonly found wish is to compare an observed non-Gaussian field to a simulated one, in order to infer parameters via this statistical comparison. The ideal case for such a comparison would be if one knew the probability distribution that generated these fields: given this distribution, one could calculate how likely one field is, given the other.
In this paper we therefore developed a highly adaptable probability distribution from which non-Gaussian fields  Fig. 7.-A non-Gaussian, and a Gaussian 60x60 field showing successive samples from the burn-in phase (left to right). Both fields have the same covariance matrix. As the sampler proceeds, the non-Gaussian field (top) builds up non-Gaussianities, and can be seen to quickly deviate from the best Gaussian approximation (bottom).
can be generated (Sect. 4). In Sect. 4.2 we specialized to cosmology, by enforcing that our fields are statistically isotropic and homogeneous. We additionally enforced in Sect. 4.3 that pixels at increasing distance become statistically ever more independent of each other.
Our field generation proceeds by MCMC sampling, and once the statistical properties of a field are defined by this distribution, our technique easily generates many thousands of fields with the prescribed statistics. However, the maximum resolution we could achieve with a university high performance cluster is 70 by 70 pixels.
In fact, we found the numerical demands in generating the non-Gaussian fields to be unconventional and therefore (for now) limiting: drawing from a joint probability distribution of 4900 pixels implicitly implies all conditional distributions of one pixel's dependency on all others are to be evaluated. We therefore found the ideal numerical setup for this calculation to be a multithread access to an unconventionally large shared memory, whereas more conventional distributed computing almost immediately produced long waiting times for accessing the required joint information.
As presented, our setup has the advantage of handling strong rather than perturbative non-Gaussianity. It is in this respect the first one to our knowledge. Although an extension to larger pixel numbers is highly desirable, we already conclude from our 70 by 70 fields that the human eye struggles to recognize non-Gaussianity, and also histograms of pixel values are an inaccurate predictor of how non-Gaussian a field is. The extraction of non-Gaussian information should therefore also be attempted on fields which appear almost Gaussian.

THE GRADIENT AND HESSIAN
The gradient and Hessian of a DALI distribution appear as intermediate steps in the HMC sampler here used. Components-wise, the gradient is given by ∂(−2 log(D)) ∂∆ e = 2 w e w i ∆ i + w e W ij ∆ i ∆ j + 2 W ei w j ∆ i ∆ j + W ei W jk ∆ i ∆ j ∆ k .
The components of the Hessian matrix then follow to be Evaluated at the peak, where ∆ ≡ 0, the Hessian is thus identical to the usual Fisher matrix. In cases of weak non-Gaussianity, it can provide a reliable estimator for the mass matrix of the sampler.

VERIFICATION OF STATISTICAL ISOTROPY AND HOMOGENEITY
Standard cosmology demands cosmic fields to be statistically homogeneous and isotropic at the largest scales. We implemented these constraints in Sect. 4.2, and verified isotropy and homogeneity as follows.
We compute the 4th multivariate central moment for irregular 4-pointed shapes of pixels. An example of such a 4-pixel shape is depicted in Figure 8. The multivariate 4th central moment is estimated as: where α, β, γ, δ denote the four corners of the shape and are held fixed. The overbar denotes the sample mean. -Illustration of statistical isotropy and homogeneity: plotted are histograms of the fourth generalized (since multivariate) central moment of four pixels. This fourth moment is measured for the non-Gaussian Q-field. The four pixels are arranged in a geometric shape as depicted in the rightmost panel. The shape is shifted and rotated around the field and the fourth moment is histogrammed for each instance. The grey histograms arises from shifting the shape as shown, the red, green and yellow histograms first rotate the shape by 90, 180 or 270 degrees and then shift all over the field. All four histograms agree, which supports statistical isotropy and homogeneity.
The fixed shape is then shifted across the field both horizontally and vertically, and the positions are enumerated by the index i. The histogram of the 4th moment at each point is shown in Fig. 8. The different thin lined histograms in colour show the same measurement, but with the shape rotated by either 90 deg, 180 deg or 270 deg. All these different rotations yielding the same result proves isotropy of the field. The 6 different panels show different separations d between the points in our 4-point shape to indicate that the field is homogeneous and isotropic on all scales. The noisiness of the histograms at larger separations is simply because a smaller number of unique shape positions fit into the field.