  We calculate the two-, three- and (for the first time) four-point correlation
functions of the COBE-DMR 4-year sky maps, and search for evidence of
non-Gaussianity by comparing the data to Monte Carlo-simulations of the
functions. The analysis is performed for the 53 and 90 GHz channels, and five
linear combinations thereof. For each map, we simulate an ensemble of 10,000
Gaussian realisations based on an a priori best-fit scale-invariant
cosmological power spectrum, the DMR beam pattern and instrument-specific noise
properties. Each observed COBE-DMR map is compared to the ensemble using a
simple chi^2 statistic, itself calibrated by simulations. In addition, under
the assumption of Gaussian fluctuations, we find explicit expressions for the
expected values of the four-point functions in terms of combinations of
products of the two-point functions, then compare the observed four-point
statistics to those predicted by the observed two-point function, using a
redefined chi^2 statistic. Both tests accept the hypothesis that the DMR maps
are consistent with Gaussian initial perturbations.
