{'fraction_non_alphanumeric': 0.035372108777702506, 'fraction_numerical': 0.015774702409755155, 'mean_word_length': 4.405963902694219, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Copula Cosmology: Constructing a Likelihood Function

To estimate cosmological parameters from a given dataset, we need to construct a likelihood function, which sometimes has a complicated functional form. We introduce the copula, a mathematical tool to construct an arbitrary multivariate distribution function from one-dimensional marginal distribution functions with any given dependence structure. It is shown that a likelihood function constructed by the so-called Gaussian copula can reproduce very well the n-dimensional probability distribution of the cosmic shear power spectrum obtained from a large number of ray-tracing simulations. This suggests that the Copula likelihood will be a powerful tool for future weak lensing analyses, instead of the conventional multivariate Gaussian likelihood.


I. INTRODUCTION
Weak gravitational lensing by intervening large scale cosmic structures provides an excellent tool to probe the nature of dark matter and dark energy. The so-called "cosmic shear" signal has been successfully measured in various groups since 2000 [e.g. [1][2][3][4][5][6][7][8][9][10]. If systematic errors are well under control, the weak lensing has the highest potential to constrain the physical parameters in the equation of state governing the dark energy among the cosmological observations, such as type Ia supernovae, baryon acoustic oscillations, and number count of galaxy clusters [11][12][13].
A number of wide-field weak lensing surveys have been planned for this purpose, such as the Subaru Hyper Suprime-Cam Survey [14], the Panoramic Survey Telescope & Rapid Response System (Pan-STARRS 1 ), the Dark Energy Survey (DES 2 ), the Large Synoptic Survey Telescope (LSST 3 ), the Joint Dark Energy Mission (JDEM 4 ) and Euclid [15]. It is expected that such wide-field weak lensing surveys will reduce statistical errors significantly compared to the past and ongoing surveys, because the number of observed galaxies increases proportional to the survey area. However, to make a maximal use of the full potential of planned weak lensing surveys for estimating cosmological parameters, it is of great importance to employ adequate statistical measures and methods for weak lensing. Particularly, one needs to take into account properly the correlations of the observables between different angular scales and redshifts, i.e. the covariances. Furthermore, we also need to use an appropriate likelihood function with given marginal distributions. If we do not use proper statistical measures and methods, or if we adopt an inaccurate covariance and/or a likelihood function, obtained results may be systematically biased [16,17].
For the cosmological parameter estimation, almost all previous authors used the χ 2 method in weak lensing analyses [e.g. [5][6][7][8][9][10]. However, it is found that the probability distribution function (PDF) of the weak lensing power spectrum is well approximated by the χ 2 distribution though it has a larger positive tail than expected from the χ 2 distribution [18]. The χ 2 distribution deviates from the Gaussian distribution on large scales because the number of modes corresponding to the degree of freedom is very small. Meanwhile, the χ 2 distribution approaches the Gaussian distribution at high ℓ due to the central limit theorem. We have to include this information accurately when we place constraints on the cosmological parameters.
If all of the marginal distributions are Gaussian distributed, it is straightforward to reconstruct the multivariate PDF, which is the so-called multivariate Gaussian PDF. However, it is not a trivial task to reconstruct the original multivariate PDF from general marginal distributions. There has been an infinite number of the degree of freedom in choosing the original PDF unless the dependence structure is specified. The copula provides us with a straightforward solution to this problem. The copula has been used in the field of mathematical finance, but not been widely used in the field of astronomy and cosmology, except only a few applications [19][20][21][22][23]. Hence, a copula may have a potential to open new fields of astronomy and cosmology.
In this paper, we construct a more plausible likelihood function using the Gaussian copula (hereafter "Copula likelihood") than the multivariate Gaussian likelihood for the cosmic shear power spectrum. We show that the Copula likelihood well reproduces the n-dimensional probability distribution of the cosmic shear power spectrum estimated from 1000 realizations which is obtained from ray-tracing simulations performed by [18]. Cosmological parameters employed for our ray-tracing simulations are consistent with the WMAP 3-year results [24]. The detail descriptions of our ray-tracing simulations are summarized in [18] (see also, [25]). In a companion paper [26], we estimate the cosmological parameters using both the copula likelihood constructed by this paper and the Gaussian likelihood in order to evaluate how the difference between two likelihoods affects the parameter estimation.

II. FORMULATION
The likelihood function is a central tool for any kind of parameter estimation. It is defined as a function of parameters in a given statistical model. The likelihood function L is related to the joint probability density function (JPDF) denoted by f (x 1 ,x 2 , . . . ,x n | θ) as wherex i (i = 1, 2, . . . , n) are independent and identically-distributed observed variables and θ denotes model parameters. We will suppress the argument θ hereafter. Recall that the n-point cumulative distribution function (CDF), denoted by F , is defined as and there is a following relation between JPDF and CDF: From the Sklar's theorem [27], we can obtain the following relation: where C denotes the function called copula and F i denotes one-point CDF defined by Therefore, the copula indicates how the one-point CDFs are jointed together to give the n-point CDF. A comprehensive proof of Sklar's theorem and rigorous definition of a copula are found in [23,28]. From the above relation, we can easily derive thatx i = F −1 i (u i ), and then derive the following relation, from Eq. (4). Differentiating Eq. (6) with Eq. (3) gives the density of copula c as where each f i is the marginal density function of the marginal CDF F i . The JPDF can then be expressed as If the variables x i are independent of each other, c = 1. In general, however, they are often correlated and c shows their correlations as a function of one-point CDF of each stochastic variable.

III. GAUSSIAN COPULA
In this section, we derive the Copula likelihood using a Gaussian copula which is more plausible than the multivariate Gaussian likelihood for the cosmic shear power spectrum. The multivariate Gaussian copula is the copula of the ndimensional random vector that is normally distributed. This copula is expressed by Here one-point Gaussian CDF Φ 1 is and Φ is n-point Gaussian CDF defined by where we consider the n-point Gaussian CDF with mean µ and n × n covariance matrix. Cov −1 shows the inverse covariance matrix. We define µ ≡ (µ 1 , µ 2 , . . . , µ n ) and x ≡ (x 1 , x 2 , . . . , x n ) and superscript 'T' stands for the transpose of vector.

A. Gaussian PDF
Let us consider a two-variables case for simplicity. First consider the case that one-point PDF is Gaussian with σ i standard deviation and µ i mean. In this case, the Gaussian copula is From Eq. (7), we derive the density of copula as By using we finally obtain the density of copula as where I stands for the identity matrix. Therefore, we can obtain the JPDF as from Eq. (8). Therefore, the JPDF using the Gaussian copula with a Gaussian one-point PDF results in the multivariate Gaussian distribution, as expected.

B. Beyond Gaussian PDF
Now, let us consider the case in which a one-point PDF is not a Gaussian distribution but a general probability distribution f i . In this case, the Gaussian copula is where one-point CDF F i is Defining , we obtain the copula density from Eq. (18) as where We use a formula of the differential of the inverse function at the second equality of the above equation. By using Eq. (13) and Eq. (21), we rewrite Eq. (20) as From Eq. (8), we can easily calculate the JPDF as We extend the above equation to the n-dimensional case and then finally obtain the n-dimensional JPDF as If f i is Gaussian, q i becomesx i . In this case, we can see that f (x 1 ,x 2 , . . . ,x n ) (Eq. 24) reduces to a multivariate Gaussian distribution. Therefore, when f i are not Gaussian PDFs, Eq. (24) carries the whole information on the correction to the Gaussian distribution. Both in practice and from theoretical point of view [29], it is more appropriate to work with the logarithm of the likelihood function, ln L, called log-likelihood. For example, the log-likelihood is used for the likelihood ratio test. The likelihood ratio test statistic is defined as twice the difference in these log-likelihoods with a minus sign. This quantity is also fundamental for the information statistics and information criterion theory (see, e.g. [29] and references therein). From Eq. (24), this test statistic is derived as for a general probability distribution. For a Gaussian case it reduces to the well-known form Here we abbreviate the irrelevant constant term in the above two equations. Finally, we give the relation betweenx i and q i because one has to calculate q i givenx i . Sincex i is related to u i through Eq. (19), we have only to derive the relation between u i and q i . From q i = Φ −1 1 (u i ), we get We change the variable x to y = (x − µ i )/σ i and then obtain At a first glance, we can recognize this equation as We can obtain the relation between u i and q i as

IV. RESULT AND DISCUSSION
In this section, we investigate whether the Copula likelihood (Eq. 25) reproduces the n-dimensional JPDF obtained from ray-tracing simulations performed by [18]. We then compare the shape of the Copula likelihood with that of the Gaussian likelihood and how the Copula likelihood is statistically better than the Gaussian likelihood. In our application, the observed variablesx i are the binned convergence power spectraP κ (ℓ i ), which are estimated from each realization. We estimate it for an assumed bin width ∆ ln ℓ. Throughout this paper, we employ the bin width ∆ ln ℓ = 0.3 and assume that the single source redshift distribution, i.e. all lensed galaxies lie at z s = 1.0 and we do not consider intrinsic ellipticity dispersion σ ǫ . In our ray-tracing simulations, the survey area is set as Ω s = 25 deg 2 . Therefore the fundamental mode of our ray-tracing simulations is at the multipole ℓ f = 72. We take 13 bins for our likelihood analysis. Therefore, we cover the multipole range from ℓ f = 72 to ℓ max = 2635. We call the convergence power spectrum estimated at bin i as "bin i power".
The left panel of Fig. 1 shows the two-dimensional JPDF between bin 1 power and bin 2 power from 1000 realizations of ray-tracing simulations. The one-point PDF is normalized so that the mean power spectrum is equal to unity. The red solid and blue solid contours show the 1σ and 2σ confidence level (CL) regions, respectively. In the right panel of Fig. 1, the blue and red contours show two-dimensional marginalized 1σ and 2σ confidence regions on the bin 1 power and bin 2 power plane, which are derived from the Gaussian likelihood (Eq. 26) and Copula likelihood (Eq. 25), respectively. One can clearly see that the Copula likelihood model reproduces the simulation data (left panel of Fig. 1) much better than the Gaussian likelihood model, as discussed below.
In our Copula likelihood model, we have chosen a χ 2 distribution for a general one-point probability distribution f i . This is a nice choice because the one-point PDF of convergence power spectrum is fairly well described by a χ 2 distribution where the mean and variance are P κ (ℓ i ) = P κ (ℓ i ) and σ 2 (ℓ i ) = P κ (ℓ i ) 2 − P κ (ℓ i ) 2 , respectively [see 18]. From [30], this χ 2 distribution is shown as forP κ (ℓ i ) > 0 and f χ 2 = 0 forP κ (ℓ i ) ≤ 0. Here, Γ(x) is the gamma function and we define Υ ≡ P κ (ℓ i ) 2 /σ 2 (ℓ i ) which corresponds to the number of independent modes. The covariance matrix is also estimated from the simulation. In order to obtain these likelihood contours, we employed the Markov Chain Monte Carlo method [31]. Assuming flat priors for bin i powers, we explored bin i power estimations in the multidimensional space (i.e. 13 dimensional space in this case). Eight parallel chains are computed and the convergence test is made based on the Gelman and Rubin statistics called "R − 1" statistics [32]. Each of our chains typically has 400,000 points and R − 1 < 0.05 in both models.
We can see that the results between our Copula likelihood function (red contours) and the Gaussian likelihood (blue contours) are very different. The difference mainly comes from the one-point PDF, which is taken as a χ 2 distribution or a Gaussian distribution. The χ 2 distribution denoted by Eq. (30) deviates from the Gaussian distribution on large scales, such as a case considered in Fig. 1, because the number of modes corresponding to the degrees of freedom is very small. We can also see that our likelihood function well reproduces the results from ray-tracing simulations than the Gaussian likelihood. In particular, the results of our Copula likelihood function capture the feature that the value which takes the maximum probability deviates from the mean value. The values at which the Copula likelihood takes its maximum are (0.619, 0.241) at bin 1 and bin 2. Figure 2 is same as Fig. 1, but results from bin 6 power and bin 7 power. The bin 6 and bin 7 correspond to the multipoles ℓ = 323 and ℓ = 436, respectively. The results from our likelihood are similar to the results from the Gaussian likelihood, because the χ 2 distribution approaches Gaussian distribution at these scales. However, we can see the small deviation from each mean value between the results obtained from the two likelihoods. Figure 3 is also the same as Fig. 1, but results from bin 12 power and bin 13 power. The bin 12 and bin 13 correspond to the multipoles ℓ = 1952 and ℓ = 2635, respectively. The contours from the copula likelihood are nearly identical to those from the Gaussian likelihood, because the χ 2 distribution approaches to Gaussian distribution at these small scales due to the central limit theorem. We can see a positive strong correlation which is mainly attributed to the fact that the nonlinear gravitational evolution causes the non-Gaussian contribution and correlations between different multipole bins (e.g., see in Fig. 7 in [18]).
To confirm a much better reliability of our Copula likelihood, we examine how the Copula likelihood is statistically better relative to the Gaussian likelihood based on the information criterion theory. We use the Akaike Information Criterion (AIC) [29,33,34] to evaluate which model is more preferable. The AIC is defined as where L max is the maximum likelihood achievable by the model under a certain dataset, and k is the number of free parameters of the model. In the present analysis L max is directly obtained by substitutingP κ (ℓ i ) into Eqs. (25) and (26) at each realization, because model parameters have been already fixed to give the maximum probability for the 1000 realized simulation data. The meaning of the AIC is clearly understood as an extension of the maximum likelihood method. An explanation for astronomers can be found in [29] and there are also useful applications of the AIC in their paper. Using the AIC enables us to compare the goodness of a certain model with that of another type directly. Figure 4 shows the PDF of the AIC difference, ∆AIC ≡ AIC Gaussian − AIC Copula among 1000 realizations. We see that the AIC is significantly reduced by the Copula likelihood (positive ∆AIC) for the vast majority of the realizations. This means that the Copula likelihood function represents the distribution of the convergence power spectra from Monte Carlo simulations better than the Gaussian likelihood one. We estimate that the mean value of the AIC difference is 2.35. Some rough rules of thumb are available and useful for estimation of the goodness of models. As shown in [35], the mean value of the AIC difference we have estimated is regarded as Considerably less which means the evidence to support the Gaussian likelihood is considerably less than the Copula likelihood. Therefore we can conclude that our Copula likelihood is strongly favored compared to the Gaussian likelihood for cosmic shear power spectrum. Finally, we discuss how the difference between the two likelihoods affects the cosmological parameter estimation, paying attention to the number of bins. Since weak lensing power spectra are expected to be very smooth in ℓ, it is very likely that little information is lost by binning as long as the bins are narrow compared to the width of any features. If we consider an observation over the sky coverage f sky (f sky = Ω s /4π) with useful signal at ℓ f < ∼ ℓ < ∼ ℓ max , the total number of modes one can obtain is estimated to be n t ∼ f sky (ℓ 2 max − ℓ 2 f ). If the number of bins is small enough, the distribution of the power at each bin can be approximated by a Gaussian because a sufficient number of modes are in that bin. As is shown in [36] if the number of bins satisfies a condition n b ≪ f 1/2 sky ℓ max assuming ℓ 2 max ≫ ℓ 2 f , the Gaussian approximation for the one-point PDF at each bin becomes very good. In that case one can use the multivariate Gaussian likelihood in order to estimate the cosmological parameters, which will simplify the parameter estimation analysis. Now let us apply this argument to our examples. In our example with f sky = 6.25 × 10 −4 and ℓ max ∼ 1000, the square root of the total number of the modes is n t ∼ f 1/2 sky ℓ max = 25. This number is not so large compared to the number of bins, i.e., n b = 13. Therefore, in this case an accurate likelihood function should be used in order to get unbiased cosmological parameter constraints, instead of the conventional multivariate Gaussian likelihood [see , 26]. Meanwhile, if we consider a future type survey which has f sky = 0.05 and ℓ max ∼ 1000 with the same number of bins, the square root of the total number of the modes is f 1/2 sky ℓ max = 224 which is much larger than the number of bins. Then, the Gaussian approximation could be fine in this case [26] 5 . Note, however, that Gaussian approximation may be violated by natural binning which has optimal ℓ resolution, regardless of how ℓ max value is. Therefore, we suggest the Copula likelihood should be used when an optimal binning is done to keep the cosmological information as much as possible. The impact of the difference between the two likelihoods on cosmological parameter estimations is illustrated and discussed carefully in a companion paper [26].

V. CONCLUSION
It has been becoming an important issue to obtain not only accurate statistics such as a power spectrum and covariance matrix but also an accurate likelihood function for the precision cosmology. In this work, we introduced a statistical tool called a copula into cosmology in a rather pedagogical way. The copula is a function to generate an n-point CDF from the given one-point CDFs and prescribed dependence structure of variables. We then applied the copula to the cosmic convergence power spectrum estimated from 1000 realizations obtained from ray-tracing simulations generated by [18].
By taking into account the fact that the one-point PDF of the convergence power spectrum is well approximated by χ 2 distribution, we showed that the Gaussian copula can reproduce the n-dimensional shape of the likelihood function of the convergence power spectrum better than the multivariate Gaussian likelihood, for the simulated 13-binned data expected from the assumed survey area of Ω s = 25 deg 2 . This is a main result of the present paper. The differences between the two likelihood models become significant at lowest multipole bins.
The deviation from the multivariate Gaussian will heavily depend on the width of the binning. We discussed that a sufficiently sparse binning will make it possible to use a multivariate Gaussian likelihood for the future cosmic shear survey. However, if one makes an optimal binning (∆ℓ ∼ 1) to keep the information as much as possible, non-Gaussian corrections will become important. In this case the copula will provide us an appropriate likelihood function in a convenient way by relating the one-point CDFs to the n-point CDF.