You are an expert at summarizing long articles. Proceed to summarize the following text: 
the temperature fluctuations in the cosmic microwave background ( cmb ) are gaussian to a high degree of accuracy @xcite . non - gaussianity , if any , enters at a highly subdominant level . it could be either primordially generated along with gaussian fluctuations by exotic inflationary models , and/or it could arise from secondary anisotropies , such as gravitational lensing , sunyaev - zeldovich ( sz ) , or sachs - wolfe ( sw ) effects . quantifying the degree and nature of non - gaussianity in the cmb constrains specific inflationary models , as well as enhances our understanding of the secondary processes the cmb underwent beyond the surface of last scattering . interpretation of any such measurement is complicated by the fact that systematics and foreground contaminations might also produce non - gaussian signatures . given the nearly gaussian nature of the cmb , @xmath0-point correlation functions , and their harmonic counterparts , polyspectra , are the most natural tools for the perturbative understanding of non - gaussianity . if it were generated by inflationary models admitting a @xmath7 term , the leading order effect would be a @xmath8-point function . on the other hand some secondary anisotropies , such as lensing , are known to produce 4-point non - gaussianity at leading order @xcite . the skewness ( or integrated bispectrum ) was measured by @xcite and @xmath8-point correlation function by @xcite . many alternative statistics have been used to investigate non - gaussianity in cmb . a partial list includes wavelet coefficients @xcite , minkowski functionals @xcite , phase correlations between spherical harmonic coefficients @xcite , multipole alignment statistics @xcite , statistics of hot and cold spots @xcite , higher criticism statistic of pixel values directly @xcite . most of these measurements are consistent with gaussianity , although some claim detections of non - gaussianity up to 3-@xmath3 level . these alternative statistics , albeit often easier to measure , typically depend on @xmath0-point functions in a complex way , thus they can not pin - point as precisely the source of non - gaussianity .    among the three - point statistics , there is a perceived complementarity between harmonic and real space methods . the bispectrum can be relatively easily calculated for a full sky map @xcite , although the present methods have a somewhat slow @xmath9 scaling @xcite . methods put forward so far use the `` pseudo - bispectrum '' , ignoring the convolution with the complicated geometry induced by galactic cut and cut - out holes . in contrast with harmonic space , the corresponding pixel space edge effect corrections are trivial @xcite , since the window function is diagonal . unfortunately , simple methods to measure three - point clustering exhibit a prohibitive @xmath2 scaling if the full configuration space is scanned . to remedy the situation , most previous measurements of the @xmath8-point function only deal with an ad - hoc sub - set of triangular configurations @xcite . both of these papers covered the full configuration space on small scales ; the former paper also appears to have estimated most configurations on large scales , missing intermediate configurations with mixes scales . this work presents a novel method , which , at a given resolution , scans the full available configuration space for @xmath8-point level statistics using realistic computational resources . we find that the resulting configuration space itself is overwhelming to such a degree that interpretation of the results also requires novel methods . we introduce false discovery rate ( fdr ) technique as a tool to interpret three - point correlation function measurements . the next section introduces our algorithm to measure the @xmath8-point correlation function , 3 illustrates it with an application to the wmap first year data release , and  4 introduces the fdr method and applies it to our results . we summarize and discuss our results in  5 . the three point correlation function ( e.g. , * ? ? ? * ) is defined as a joint moment of three density fields @xmath10 at three spatial positions . for cmb studies @xmath11 denotes temperature fluctuations at position @xmath12 on the sky , and @xmath13 stands for ensemble average . if the underlying distribution is spatially isotropic , @xmath14 will only depend on the shape and size of a ( spherical ) triangle arising from the three positions . a number of characterizations of this triangle are possible and convenient . the most widely used are the sizes of its sides ( measured in radians ) , or two sizes and the angle between them . this latter angle is measured on the spherical surface of the sky . one can use the ergodic principle of replacing ensemble averages with spatial averages to construct a nearly optimal , edge corrected estimators with heuristic weights @xcite @xmath15 where we symbolically denoted a particular triangular configuration with @xmath16 ( any parametrization would suffice ) , and @xmath17 if pixels @xmath18 , and @xmath19 otherwise . we also defined a multiplicative weight @xmath20 for each pixel : this is @xmath19 if a pixel is masked out , and it could take various convenient values depending on our noise weighting scheme if the pixel is inside the survey ; e.g. , in the case of flat weights it is simply @xmath21 . this simple estimator has been widely used in large scale structure , and it is nearly optimal with appropriate weights . ( e.g. , * ? ? ? * ; * ? ? ? it is entirely analogous to the successful estimators used for the measurements of the @xmath22 s for the cmb ( up to harmonic transform , * ? ? ? * ; * ? ? ? the naive realization of equation  [ eq : estimator ] has a prohibitive @xmath2 scaling if one needs to scan through triplets of pixels and assign them to a particular bin . the summation can be restricted and thus made faster if one restricts the number of configurations and the resolution ( e.g. , * ? ? ? * ; * ? ? ? * ; * ? ? ? * ) , or it can be sped up by using tree - data structures @xcite . neither of these methods is able to scan through all possible configurations in megapixel maps with reasonable amount of computing resources . here we propose a new method which uses both hierarchical pixelization and fourier methods motivated by @xcite to scan through all the triangles simultaneously . note that @xcite comes closest to our aims , but their simple two - step approach is not systematic enough to cover all possible triangles at a given resolution , and it is not fast enough for massive monte carlo simulations .    in the following we will choose a parametrization of the triangle @xmath16 using two of its sides @xmath23 , and the angle @xmath24 between them . we define the configuration space as a set of ( logarithmic ) bins for the sides , and linear bins for the angle in their full possible range , i.e. , @xmath25 ( remember that the sides of the triangle on the sky are also measured in radians ) . the given resolution is determined by the number of bins for @xmath26 , and the number of bins for @xmath24 . note that a particular triangle might appear more than once in this scheme , albeit with different resolutions . different triangular bins of the three - point function are strongly correlated anyway , and the correlation from duplicating triangles can be taken into account in the general statistical framework over correlated bins .    given a triangular configuration , and a pixel @xmath12 , all other pixels which enter the summation in equation  [ eq : estimator ] are located on two concentric rings of size @xmath27 and @xmath28 . as a consequence , the summation over fixed @xmath24 can be thought of as an unnormalized ( raw ) two - point correlation function between two rings . to obtain three - point correlation function , one has to multiply this two - point correlation function with the value of the center pixel @xmath12 and finally sum over @xmath12 . calculating the two - point correlation function of rings can be fast if one repixellizes the map ( c.f . [ fig : rings ] ) into rings with sizes matching the binning scheme for @xmath29 , and uniform division in @xmath24 . such a repixellization , resulting in ring - pixels as shown in fig . [ fig : rings ] , would take only @xmath0 steps even in a naive way ; the healpix hierarchical scheme allows it to be done in @xmath30 time . we use the following algorithm : let us start a recursive tree walk at the coarsest map , @xmath31 in the healpix scheme . for each pixel in this map , we determine , using its center , which ring - pixel it would belong to . if the size of the pixel is much smaller than this ring - pixel ( how much smaller is a parameter of our algorithm : in this paper we used the condition that the pixel has to be smaller then @xmath32 the bin width which is also the approximate size of the ring - pixels ) , we record it . if not , algorithm splits the quad - tree , and calls itself recursively for each four sub - pixels . this procedure ends at the latest when the highest resolution ( i.e. the one of the underlying map ) is reached . if the bins are chosen appropriately such that large ring - pixels are set up for large triangles , for many pixels it will finish earlier . as noted above , the map has to be regridded around each pixel into rings of ring - pixels . in total , this takes @xmath33 time .    calculating the two - point correlation function between rings speeds up using fast fourier transform ( fft ) methods , such as those put forward in @xcite . the recipe is the following . first , fft every ring @xmath12 to obtain complex coefficients @xmath34 ; then calculate for every pair of rings @xmath35 the `` pseudo power spectrum '' @xmath36 , where @xmath37 means complex conjugate . due to the u(1 ) symmetry of the ring an inverse @xmath38 transform will give the ( raw ) two - point correlation function between the two rings ( c.f . * ? ? ? * ; * ? ? ? if we have @xmath39 rings , each of them @xmath40 ring - pixels , each fft can be done in @xmath41 time , and there is @xmath42 cross correlations to be calculated for a full scan of configurations . all the above needs to be performed for each pixel as a center point . the total scaling ( including the initial regridding ) takes @xmath43 , where we took into account that the two opposite pixels can be handled in one go if a symmetric set of bins around @xmath44 are used for @xmath29 . while the above procedure to calculate raw ( unnormalized ) correlation functions appears somewhat complex , we have checked with direct calculation that it gives numerically the same result as calculating correlations on the rings in a naive way . in order to obtain normalized correlation functions , the same procedure has to be followed for the rings associated with weights / masks . each configuration of the raw three - point function is divided with the mask / weight three - point function for the final result . for many realizations with the same mask , such as in the case of massive monte carlo simulations , the mask correlations need to be estimated only once , representing negligible cost . the above abstract scheme and calculation will be illustrated and further clarified with a practical application to wmap next . we demonstrate our method to calculate the three - point correlation function with an application to wmap . we downloaded first year foreground cleaned maps from lambda website@xcite . there are total 8 maps for 8 differencing assemblies ( da ) in q , v and w bands : q1 , q2 , v1 , v2 , w1 , w2 , w3 , and w4 , already in healpix format . following the two - point analysis of wmap @xcite , we only used cross correlations , i.e. three - point correlation functions calculated from three different das . we produced 100 ( gaussian ) simulations with synfast in healpix package . the input power spectrum , also from the lambda website , was taken from @xmath6cdm model using a scale - dependent ( running ) primordial spectral index which best fits the wmap , cbi and acbar cmb data , plus the 2df and lyman - alpha data . every simulation consists of 8 assembly maps as the data @xcite . these 8 maps were generated with a same random seed , representing the same primordial cmb , but 8 different beam transfer functions . then different simulated noise maps from lambda web @xcite were added to the synfast output maps .    since the non - gaussian signal is exceedingly small , and on the smallest scales the data are noise dominated , we degraded all maps ( simulations and data ) to @xmath45 after applying the kp2 mask . more precisely , we added up pixel and weight values for each map ; our two weighting schemes are presented in the next subsection .      at the heart of our algorithm is the regridding of figure  [ fig : rings ] which matches our binning of the triangles . we chose 19 rings for half of the sphere surface , and the same bins are repeated on the other half symmetrically around @xmath44 . the 19 bins are chosen to be uniformly distributed in logarithm between @xmath46 and @xmath47 . the number 19 was chosen such that @xmath48 , which gives a logarithmic resolution of @xmath49 . every ring was divided to 20 ring - pixels in @xmath24 . this number renders the resulting ring - pixels fairly compact , and it is also convenient for our chosen implementation of fft ( fftw , * ? ? ? weight maps were constructed using the kp2 mask and the noise profile of the maps . we used two weighting schemes : flat weighting where @xmath20 is 1 or 0 depending on the mask , and ( inverse ) noise weighting ; for the latter @xcite we used the effective number of observations of the pixel . the weights need to be determined only up to multiplicative factor , as their overall normalization cancels from the algorithm . the average noise level @xmath50 for each da is used when combining @xmath14 over different cross correlations . the total number of triangular configurations in @xmath51 rings with @xmath52 possible values of @xmath24 ( angles large than @xmath53 count to @xmath54 ) is @xmath55 for autocorrelations . the same number is valid for a cross correlations ( which we will be exclusively doing ) of 3 das . we introduce the notation ( da1 , da2 , da3 ) , for the central pixel , the first ring and the second ring sampled from the three da in this order . in addition we restrict the `` first ring '' has @xmath26 no larger than that of the `` second ring '' . then the total number of cross correlations between the 8 da s is @xmath56 . effectively , each triangle is calculated six times for a given triplet of da s due to the possible 6 permutations . however , each sample has a different resolution , therefore we opt to keep all possibilities . the resulting correlations are taken into account when dealing with correlated bins in general . in total , there are about @xmath57 triangular configurations for each data or simulation set . note that the total number of triplets of ring - pixels examined is @xmath58 more , or @xmath59 . this is still a lot smaller than checking @xmath60 triplets of pixels naively in an @xmath61 map . these numbers suggest that our algorithm even without fft should take order of days , while the naive algorithm would need over 200 years of cpu .    for a batch of @xmath62 simulations ( comprising of @xmath63 da maps ) , the calculation of the full three - point function in all the configuration takes about 90 hours on an intel xeon 2.4ghz cpu . this means each cross - correlation takes only about @xmath64 minutes on average ! about 10 hours are saved by batch processing 10 sets of simulations : to estimate @xmath14 in the data alone ( one set of 8 da files ) took 10 hours . clearly , the given resolution does not extract all the information from the data , as there are approximately @xmath65 distinct configurations of the bispectrum or three - point function . however , it surely must be redundant to extract more configurations than the amount of data . the ratio of data points vs. configurations is about @xmath66 for our chosen bins . it is unlikely that it were fruitful to push this number towards much smaller values , although the speed of our algorithm would allow higher resolution . figure  [ fig : zeta ] shows a typical set of @xmath14 measurements for the ( w2,w3,w4 ) da cross - correlation . the results from wmap lie comfortably in the 68% range of results from gaussian simulations . similar results are found for the other da combinations , or when all the 24 possible w - only da results are averaged . although different combinations have different effective beam , full averaging is meaningful on large scales . we checked that averaging all 336 possible combinations is also consistent with gaussian . finally , repeating all the measurement with noise weighting produced no obvious departure from gaussianity either . at the same time , the scatter in the simulations , i.e. , the probability density function ( pdf ) of @xmath14 from 100 simulation , shows a slightly non - gaussian signature . for ( w2,w3,w4 ) , figure  [ fig : hist ] shows a histogram derived from all simulated @xmath14 values normalized by their measured median and 68% levels . slight deviations from gaussian distribution are evident : student distribution of degree 3 fits better the overall distribution . this same distribution produces lower @xmath67 when applied to _ individual _ triangular configurations according to the inset , i.e. it is a ( marginally ) better fit than gaussian . we fully take this into account in our hypothesis testing which is described next . our goal is to test the null - hypothesis of gaussianity against our measurements by means of comparing the @xmath14 values measured from the data with the corresponding probability distribution function ( pdf ) determined from gaussian simulations . a crucial step in the traditional @xmath67 method appears to be computationally infeasible due to the large number of configurations : calculation of the ( pseudo ) inverse of an @xmath68 matrix for @xmath69 , the total number of our highly correlated configurations . moreover , as seen above , the underlying pdf marginally violates gaussian assumption , even for gaussian simulations . even if it were possible to calculate the inverse of the covariance matrix , and we were to accept the accuracy of the gaussianity in pdf of the individual bins , it is not possible to determine the underlying covariance matrix with sufficient accuracy . in fact , one would need ( e.g. , * ? ? ? * ) at least ( and likely much more than ) 2.6 million simulations for that purpose . @xcite have shown that using simulations with uncorrelated noise might result in spurious detection of non - gaussianity . therefore we chose to use only the wmap supplied correlated noise simulations , of which 110 is available at present . it is straightforward to test the null - hypothesis with a single configuration : we can calculate a @xmath70-value from the best fit student distribution from our simulations . the @xmath70-value is defined as the probability of obtaining a @xmath14 value that is at least as extreme as the one measured from wmap . for a threshold @xmath71 , the null - hypothesis is rejected at @xmath72 level if the @xmath70-value of the datum is smaller than @xmath71 . a problem arises from combining 2.6 million tests when all the data are used . for instance , even if the hypothesis were true , about 260 bins would still be rejected at the @xmath73 level ( ignoring the correlations in the data ) . fortunately a robust and simple method exists for massive hypothesis testing , which is insensitive to correlations between the tests , and makes no assumption on the gaussianity of the underlying error distribution : the method of false discovery rate ( fdr ) @xcite . in astronomy , it has been successfully applied in the context of image processing and finding outliers by @xcite , which can be consulted for a more detailed introduction . the fdr method combines the same @xmath70-value as defined above for individual tests using a threshold for rejecting the null hypothesis . this combination is insensitive to correlations and has more statistical power than naive combination . our goal is to adapt this powerful method for hypothesis testing of three - point correlation function measurements with overwhelming number of configurations . the fdr method gives a simple prescription for finding a threshold for rejection . in particular , the recipe suggests that we choose a threshold such that we control the rate of _ false rejections _ or fdr . the parameter , taking a similar role to the confidence interval in more traditional tests , is the maximum rate of fdr . if we fix an @xmath24 such that @xmath74 , the fdr procedure will guarantee @xmath75 in ensemble average . next we describe the recipe to control fdr ; more details can be found in @xcite . let @xmath76 denote the @xmath70-values calculated from the measurements of @xmath0 configurations , _ sorted _ from smallest to largest . let @xmath77 where @xmath78 is a constant depending on the level of correlations between different configurations . for uncorrelated data @xmath79 ; while @xmath80 can be used for correlated data @xcite . note that technically one would have to adjust @xmath78 to the degree of correlations in the data . the suggested value for correlated data is extremely conservative , and should be considered as a strong upper limit . even using this conservative adjustment decreases the statistical power of the technique only logarithmically ; the final results are expected to be robust regardless of the degree of correlations . if configurations with @xmath81 are rejected , equation  [ eq : fdr ] will hold @xcite , i.e. the fdr is controlled according to our preset parameter @xmath24 . the procedure is represented graphically on figure  [ fig : fdr ] : @xmath82 is plotted against @xmath83 superposed with the line through the origin of slope @xmath84 . all @xmath70-values reject the null hypotheses which are to the left from the last point at which @xmath82 falls below the line . these might include some false discoveries which are guaranteed to be a smaller fraction than @xmath24 in ensemble average . we have applied the fdr recipe to all of our individual cross - three point functions , as well as our full data set . since @xmath78 is a constant , initially we kept @xmath79 . for a fixed @xmath24 , the results can be subsequently reinterpreted in terms of any @xmath85 . for da combination ( w2 , w3 , w4 ) , there is no rejection for @xmath86 , i.e. , allowing as high as 81% false rejections , not a single configuration rejected our gaussian null hypothesis . correlations might increase @xmath78 , but it must be @xmath87 . the true @xmath24 , when correlations are taking into account , can only be larger than our effective @xmath24 for @xmath79 . in other words , the data are fully consistent with gaussianity .    as a sanity check , we repeated the fdr analysis in our simulations as well . by scanning through different @xmath24 values from 0 to 1 , we find that 50 out of 100 simulations have rejections with @xmath86 . this means that the wmap measurements are fully consistent with gaussianity at a level better than 1-@xmath3 in the traditional sense . in summary , at the three - point level , scanning all configurations , we did not find any significant non - gaussianty which would be localized in pixel space triangular configurations . we performed fdr analysis on all @xmath88 measurements individually , as well as on the combination of all those measurements with @xmath89 million configurations in total . none of these cases produced credible evidence for non - gaussianity and all of them were fully consistent with our null hypothesis at @xmath90 . we presented a new method to measure angular three - point correlation functions on spherical maps . we achieve an unprecedented @xmath1 scaling with a combination of hierarchical and fourier algorithms . the speed of our technique allows a systematical scan of the full available configuration space at a given resolution . such speed is especially useful for cross correlations and monte carlo simulations , where a vast number of configurations and measurements need to be performed . we have achieved a speed of about @xmath64 minutes per cross correlations , when 336 cross - correlations have been estimated simultaneously in @xmath45 healpix maps using a single intel xeon 2.4ghz cpu . this is to be contrasted with a naive approach , which would have taken about 200 years per cross - correlations ; a @xmath91 million fold speed up .    as a first application of our code we analyzed the first year wmap data along with 100 realistic simulations . we have calculated cross - correlations for about @xmath92 triangular configurations , or about @xmath59 triplets in total in maps of @xmath45 corresponding the @xmath93 da s . the ratio of pixels / configurations is about 50 for each measurement . comparing our measurements from 100 gaussian simulations with realistic correlated noise , we found wmap to be comfortably within the 68% percent range for most configurations . any significant departure from gaussianity at the three - point level , even if localized in particular triangular configurations , would have shown up clearly in our full scan of the available configuration space . our main result is that there is no credible evidence of non - gaussianity at the three - point level at any of the triangular configurations we examined . as a consequence , if the tentative detection of non - gaussianity claimed in previous works holds up , it should correspond to either 4-point or higher order correlations , or to spatially localized features which break rotational invariance ( e.g. , * ? ? ? * ; * ? ? ? in contrast with our measurements , all previous studies of higher order statistics used autocorrelations . comparison of our errorbars with that of @xcite appears to show that this increases the errors by a factor of two ( see the discussions below ) . in addition , many measurements used uncorrelated noise simulations . according to the findings of @xcite , this might increase the likelihood of finding spurious non - gaussianity . analysis of our gaussian simulations revealed that there is a slight non - gaussianity in the error distribution of individual configurations . this is not surprising , since three - point correlation function is a non - linear construction of the gaussian random variables ( c.f . the error distribution is well fit by a student distribution with 3 degrees of freedom . to quantify any possible departure of the overall data set from gaussianity , we introduced a new technique , fdr , to interpret three - point statistics . this corresponds to an optimized multiple hypothesis testing , and it is insensitive to the unavoidable correlations in the data . all of our fdr tests , whether applied to any of the 336 cross correlations , or the combined data set , were fully consistent with gaussianity with better then 1-@xmath3 . this quantifies our previous assertion based on examination of the individual configurations under the assumption of statistical isotropy . the above model independent tests showed that there is no credible evidence of any non - gaussianity in the data . next we illustrate how our measurements yield constraints on specific non - gaussian models . we choose a simple phenomenological model corresponding to the quadratic expansion of the density field in terms of one parameter , @xmath4 , as put forward by @xcite : @xmath94    to obtain constraints on this parameter , we construct an estimator for @xmath4 @xmath95 where @xmath14 is our measurement in data or simulation maps . we calculated the two - point correlation function @xmath96 analytically , to avoid any bias from the non - linear construction ( c.f . we used the same best fit power spectrum as for the simulations , as well as taking into account beam and pixel window functions . since previous measurements already established the weakness of non - gaussianity , our gaussian simulations should be accurate enough to calculate the variance @xcite . applying the same estimator to our 100 simulations , we obtained error bars for @xmath4 estimated from each particular configuration . the simplicity of the phenomenological model lies in the fact that a constant value of @xmath4 is assumed . we do not attempt to combine our estimates optimally , instead we use simple considerations . the signal increases towards small scales in this model , while noise dominates on the smallest scales . since we already discarded the smallest scales when using @xmath61 , it is intuitively clear that most signal pertaining to this model will be concentrated in the small fraction of the triangles corresponding to small scales , in particular the skewness . to confirm this we generated and analyzed a set of non - gaussian simulations according to equation  [ eq : nong ] , with @xmath4 equal to 1000 , 2000 , 3000 , 4000 , and 6000 . inspection of the configurations together with the errorbars from the gaussian simulations confirmed the above idea . therefore we decided to use the skewness , which corresponds to giving zero weight to all other configurations when combining our @xmath4 estimators . from these we obtain @xmath97 where the errorbar was estimated from the gaussian simulations . @xcite quotes similar constraints for a low quadrupole cdm model , but their errorbars are a factor of two larger for a @xmath6cdm model similar to the one we use . the fact that we obtained a factor of two tighter constraints than @xcite suggests that using cross - correlations is superior to auto - correlations for three - point statistics of wmap    as a sanity check , we calculated the mean value of the skewness estimator for 100 gaussian simulations ; it yields about @xmath98 . on the other hand , we also demonstrate that we can recover @xmath4 from the non - gaussian simulations . all the simulations have the same underlying gaussian signal and noise , the only difference is the value of @xmath4 . according to figure  [ fig : fnlt ] the errors might be underestimated when @xmath99 , and/or there might be a small low bias , but it is clear from the figure that we could detect non - gaussianity if it were present .    a suboptimal combination of @xmath4 estimates from all configurations weighted by their inverse variance yields about @xmath100 , a significantly weaker result , confirming the intuitive idea that most of the signal is concentrated on small scales .    some of the results in this paper have been derived using the healpix @xcite package . we acknowledge the use of the legacy archive for microwave background data analysis ( lambda ) . support for lambda is provided by the nasa office of space science . the authors were supported by nasa through aisr nag5 - 11996 , and atp nasa nag5 - 12101 as well as by nsf grants ast02 - 06243 , ast-0434413 and itr 1120201 - 128440 .
we present a new method to estimate three - point correlations in cosmic microwave background maps . 
 our fast fourier transform based implementation estimates three - point functions using all possible configurations ( triangles ) at a controlled resolution . 
 the speed of the technique depends both on the resolution and the total number of pixels @xmath0 . 
 the resulting @xmath1 scaling is substantially faster than naive methods with prohibitive @xmath2 scaling . as an initial application , 
 we measure three - point correlation functions in the first year data release of the wilkinson anisotropy probe . 
 we estimate 336 cross - correlations of any triplet of maps from the 8 differential assemblies , scanning altogether 2.6 million triangular configurations . 
 estimating covariances from gaussian signal plus realistic noise simulations , we perform a null - hypothesis testing with regards to the gaussianity of the cosmic microwave background . 
 our main result is that at the three - point level wmap is fully consistent with gaussianity . to quantify the level of possible deviations , we introduce false discovery rate analysis , a novel statistical technique to analyze for three - point measurements 
 this confirms that the data are consistent with gaussianity at better than 1-@xmath3 level when jointly considering all configurations . 
 we constrain a specific non - gaussian model using the quadratic approximation of weak non - gaussianities in terms of the @xmath4 parameter , for which we construct an estimator from the the three - point function . 
 we find that using the skewness alone is more constraining than a heuristic suboptimal combination of all our results ; our best estimate is @xmath5 assuming a @xmath6cdm concordance model .