Issue 
A&A
Volume 642, October 2020



Article Number  A156  
Number of page(s)  10  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202038585  
Published online  15 October 2020 
Modeling highdimensional dependence in astronomical data
^{1}
Chip Computers Consulting s.r.l., Viale Don L. Sturzo 82, S.Liberale di Marcon, 30020 Venice, Italy
email: robertovio@tin.it
^{2}
Mathematical Institute, Leiden University, Niels Bohrweg 1, 2333, CA Leiden, The Netherlands
email: t.w.nagler@math.leidenuniv.nl
^{3}
ESO, Karl Schwarzschild strasse 2, 85748 Garching, Germany
email: pandrean@eso.org
Received:
5
June
2020
Accepted:
4
August
2020
Fixing the relationship of a set of experimental quantities is a fundamental issue in many scientific disciplines. In the 2D case, the classical approach is to compute the linear correlation coefficient ρ from a scatterplot. This method, however, implicitly assumes a linear relationship between the variables. Such an assumption is not always correct. With the use of the partial correlation coefficients, an extension to the multidimensional case is possible. However, the problem of the assumed mutual linear relationship of the variables remains. A relatively recent approach that makes it possible to avoid this problem is the modeling of the joint probability density function of the data with copulas. These are functions that contain all the information on the relationship between two random variables. Although in principle this approach also can work with multidimensional data, theoretical as well computational difficulties often limit its use to the 2D case. In this paper, we consider an approach based on socalled vine copulas, which overcomes this limitation and at the same time is amenable to a theoretical treatment and feasible from the computational point of view. We applied this method to published data on the nearIR and farIR luminosities and atomic and molecular masses of the Herschel reference sample, a volumelimited sample in the nearby Universe. We determined the relationship of the luminosities and gas masses and show that the farIR luminosity can be considered as the key parameter relating the other three quantities. Once removed from the 4D relation, the residual relation among the latter is negligible. This may be interpreted as the correlation between the gas masses and nearIR luminosity being driven by the farIR luminosity, likely by the star formation activity of the galaxy.
Key words: methods: data analysis / methods: statistical
© ESO 2020
1. Introduction
Modeling the relationship of a set of experimental quantities is not straightforward. Often, no theoretical hints are available that would allow us to fix the dependence among the involved variables. Hence, the work has to be entirely based on the analysis of the data. In the 2D case, an example is represented by the scatterplots and the computation of the corresponding linear correlation coefficients ρ. Its extension to the multidimensional case is possible with the partial correlation coefficients. The main limitation of this approach is the implicit assumption of linear relationships among the variables under study. This is often an unrealistic condition. For this reason, a relatively recent alternative consists of modeling the joint probability distribution function (PDF) of the data. However, this task is not trivial even in the 2D case. Families of bivariate PDFs are available (Balakrishnan & Lai 2010), but are not very flexible and are difficult to use. Things worsen for the multidimensional case (Kotz et al. 2000). A relatively recent alternative is based on copulas. These are simply multivariate cumulative distribution functions (CDF) with standard uniform margins. They are used to describe the dependence between random variables, and their main role is to disentangle margins and the dependence structure (Nelsen 2006; Durante & Sempi 2016; Hofert et al. 2018). With copulas it is possible to decompose a joint probability distribution into their margins and a function that couples them. The copula is that coupling function.
In cosmology, 2D copulas have been used by Scherrer et al. (2010) for the determination of the PDF of the density field of the largescale structure of the Universe, by Lin & Kilbinger (2015) and Lin et al. (2016) to predict weaklensing peak counts, and by Sato et al. (2010, 2011) for the precise estimation of cosmological parameters. Other astronomical applications include the determination of the farUV (FUV) and farIR bivariate luminosity function of galaxies (Takeuchi 2010; Takeuchi et al. 2011), the determination of the Kband and the submillimeter luminosity function (Andreani et al. 2014), and the bivariate luminosity versus the mass functions of the of the local HRS galaxy sample (Andreani et al. 2018).
In principle, the copula approach can work with multidimensional data, but theoretical as well computational difficulties often limit its use to the 2D case. Recently, however, vine copulas have been proposed in the statistical literature as an approach that overcomes this limitation and at the same time is amenable to a theoretical treatment and feasible from the computational point of view. The strength of vine copulas is that they allow, in addition to the separation of margins and dependence by the copula approach, tail asymmetries and separate multivariate component modeling. This is accommodated by constructing multivariate copulas using only bivariate building blocks, which can be selected independently. These building blocks are glued together to form valid multivariate copulas by appropriate conditioning (Joe 2015; Czado 2019). This makes vine copulas a very flexible and reliable tool even in the case of very highdimensional data.
For this paper, we made use of multidimensional copulas, described in Sects. 2 and 3, and in particular of vine copulas, outlined in Sects. 4 and 5. We applied them to a data set related to a complete nearby sample of galaxies that has been observed at various wavelengths (Andreani et al. 2018, and references therein) and show its use to highlight the relation to the physical properties of the galaxies in Sect. 6.
2. What are copulas?
A ddimensional copula C_{1, …, d}(u), u = (u_{1}, …, u_{d})∈[0, 1]^{d} is simply a multivariate CDF with standard uniform univariate margins. Its importance is due to Sklar’s theorem: for any ddimensional CDF F(x), x = (x_{1}, …, x_{d})∈ℝ^{d}, with univariate margins F_{1}(x_{1}),…,F_{d}(x_{d}), a ddimensional copula C_{1, …, d}(u):[0, 1]^{d} → [0, 1] exists, such that
where u_{1} = F_{1}(x_{1}),…,u_{d} = F_{d}(x_{d}). The converse also holds, i.e. given a ddimensional copula C_{1, …, d}(u) and univariate CDFs F_{1}(x_{1}),…,F_{d}(x_{d}), the CDF F(x) defined by Eq. (1) is a ddimensional CDF with margins F_{1}(x_{1}),…,F_{d}(x_{d}). This means that copulas are those functions which combine the univariate margins F_{1}(x_{1}),…,F_{d}(x_{d}) to form the ddimensional CDF F(x). In other words, copulas link multivariate CDFs to their univariate margins. The importance of copula is more evident if the PDFs f(x) are considered. Indeed, it can be shown that
where
From Eq. (2), any joint PDF f(x) can be factorized into the product of two terms. One is the product of the marginal PDFs {f_{i}(x_{i})} and the other is the copula density c_{1, …, d}(u). The first term provides information on the statistical properties of the individual random variables {x_{i}} whereas the second term provides information on their mutual dependence. Therefore, the importance of c_{1, …, d}(u) lies in the fact that it describes the dependence structure among the random variables in separation of the associated marginal PDFs.
If a set of nddimensional random data {x_{k}}, k = 1, …, d, with x_{k} = {x_{p, k}}, p = 1, …, n, is available and the margins {F_{k}(x_{k})} and the corresponding PDFs {f_{k}(x_{k})} are known, the standard procedure to estimate f(x_{1}, …, x_{d}) is as follows: first, compute the standard uniform variates u_{k} = F_{k}(x_{k}), then fit their joint distribution by a copula C_{1, …, d}(uθ), which belongs to a continuous parametric family with characteristic parameters θ = {θ_{1}, …, θ_{np}}. After that, Eq. (2) provides the joint PDF. A common method for fitting the copula is based on an estimate of the parameters θ through a maximum likelihood method, but other techniques are also possible (Hofert et al. 2018). Often, however, the margins are not available. In this case, an alternative is to fit each set x_{k} with a PDF belonging to the Johnson, generalized Lambda, or any other family of parametric PDFs (Vio et al. 1994; Karian & Dudewicz 2011) (see also Appendix A), to compute the uniform random variates u_{k} = F_{k}(x_{k}) and then, as before, to fit a copula. When the margins are not estimable with sufficient accuracy (e.g., because of little available data), a useful nonparametric variant is the computation of the random variates u_{k} by means of the so called pseudoobservations u_{p, k} = R_{p, k}/(n + 1) with R_{p, k} the rank of x_{p, k} among (x_{1, k}, …x_{n, k}). Since in general one has no indication of which kind of copula is suited for the data of interest, the typical solution is to fit a set of copulas and to choose that which provides the best result.
In principle, the above procedures can be applied to any ddimensional data set. The point is that most of the parametric copula families available in literature are 2D (e.g., see Joe 2015), and the few available for a multidimensional analysis are not flexible enough. An alternative approach based on a nonparametric copula estimate has been also proposed (e.g., Nagler & Czado 2016; Nagler et al. 2017).
3. Preliminary considerations
Given that most of the available copula families are 2D, it is unclear how a ddimensional PDF f(x_{1}, …, x_{d}) can be computed. A possible solution is to express Eq. (2) in terms of 2D copulas. The starting point is that f(x_{1}, …, x_{d}) can be factorized into the form
with f(x_{k}y) being the conditional PDF of the random variable x_{k} given the vector of random variables y. Now, it can be proved (Czado 2019) that
where c_{xkyjy−j}(., .) is the conditional copula density,
C_{xkyjy−j}(., .) is the conditional copula, y_{j} is one arbitrarily chosen component of y, and y_{−j} denotes the yvector, excluding this component. The key point is that these conditional PDFs are expressed in terms of 2D copula densities. The same holds for the PDF f(x_{1}, …, x_{d}). For example, in the 3D case it is
In fact, the decomposition (4) is not unique since the indices of the variables {x_{k}} can be permuted. For instance, a decomposition equivalent to (7) is
Although the problem of estimating the PDF f(x_{1}, …, x_{d}) has been simplified by means of Eqs. (4)–(6), it is still hard to deal with. The conditional copulas C_{xkyjy−j} and corresponding densities c_{xkyjy−j} are difficult to estimate. For this reason, usually the conditional copula densities are simplified into the form
Something similarly occurs to the corresponding conditional copulas. This simplification does not only make the problem easier to deal with, but it permits the use of the large set of available continuous parametric families of 2D copulas. This makes the method quite flexible. For instance, in the 3D case, the decomposition can be written in the form
where the 2D copula densities c_{12}(., .; θ_{12}), c_{23}(., .; θ_{23}) and c_{132}(., ..; θ_{132}) can be chosen of different types.
4. Vine copulas: The theory
For highdimensional distributions, there is a huge number of possibilities for decompositions into 2D copulas, named paircopulas, like Eqs. (7) and (8). All these possibilities can be organized according to graphical models called "regular vines". Two special cases, called Dvine and Cvine (Aas et al. 2009), have been introduced as a simplification. Each model gives a specific way of decomposing a density.
Figure 1 shows the graphical structure of a Dvine for a fourdimensional problem. This structure is formed by three levels or trees. Each circle or ellipsis constitutes a node, and each pair of nodes is joined by an edge. The label of a node in a given tree is given by the label of the edges of the tree at its immediate left. The label of an edge is given by the indices contained in the joined nodes with the conditional index given by the common one. For example, in the central tree node (1, 2) is connected to node (2, 3). The common index is 2, hence the label of the joining edge is (1, 32). Each edge represents a paircopula density, and the edge label corresponds to the subscript of the paircopula density. The indices of the CDFs that appear as the argument of a specific paircopula density are given by the labels of the nodes connected by the corresponding edge. According to this rule, the first tree produces the terms c_{12}(F_{1}(x_{1}),F_{2}(x_{2})), c_{23}(F_{2}(x_{2}),F_{3}(x_{3})) and c_{34}(F_{3}(x_{3}),F_{4}(x_{4})). The second tree produces the terms c_{132}(F(x_{1}x_{2}),F(x_{3}x_{2})) and c_{243}(F(x_{2}x_{3}),F(x_{4}x_{3})). Finally, the last tree produces the term c_{1423}(F(x_{1}x_{2}, x_{3}),F(x_{4}x_{2}, x_{3})). The decomposition of f(x_{1}, x_{2}, x_{3}, x_{4}) is given by the product of these terms:
Fig. 1.
Example of tree structure of a 4D Dvine copula (see text). 

Open with DEXTER 
For a ddimensional density f(x_{1}, …, x_{d}), this procedure provides the decomposition formula
where index j identifies the trees, while index i runs over the edges in each tree.
Figure 2 shows the graphical structure of a Cvine again for a 4D problem. While in a Dvine no node in any tree is connected to more than two edges, in a Cvine each tree has a unique node, known as the root node, which is connected to all the other nodes. The rules for labeling the nodes and the edges are identical to those of the Dvine. For a Cvine, the decomposition formula is
Fig. 2.
Example of tree structure of a 4D Cvine copula (see text). 

Open with DEXTER 
Although in principle the decompositions provided by the Cvines and the Dvines should be equivalent, things are actually different because of the simplification (9). In general, Dvines are often useful when there is a natural ordering of the variables (e.g., by time), whereas Cvines might be advantageous when a particular variable is known to drive the interactions of the other variables. In such a situation, this variable can be located at the root node of the leftmost tree. In many practical applications, however, no a priori information is available allowing us to decide which kind of vine to use. As a consequence, the decision has to be based on which model better fits the data.
5. Vine copulas: Computational issues
The flexibility of vine copulas complicates the parameter estimation and the model selection. One needs to select the appropriate parametric families for each paircopula, estimate the parameters, and find a good structure for the vine trees. Thankfully, these problems can mostly be solved in separation per paircopula and per tree level.
We remind the reader that a copula C_{1, …, d}(u) is the distribution function of a random variate u = (u_{1}, …, u_{d}). In what follows, we assume that for all variables i = 1, …, d, n uniform variates u_{k} = {u_{p, k}}, p = 1, …, n, are available. As mentioned in Sect. 2, these are commonly obtained by transforming the original data {x_{k}} by means of u_{k} = F_{k}(x_{k}).
5.1. Model fitting in the 2D case
We first consider the simpler 2D case. Let us suppose that we have available the variates {(u_{p, 1}, u_{p, 2})}, p = 1, …, n, from a parametric copula model c_{12}(u_{1}, u_{2}; θ). Then, the parameters θ can be estimated by maximumlikelihood:
In practice, since the true copula is unknown, it is necessary to choose a parametric copula density from a set of families {ℱ_{κ}}, κ = 1, …, m, with n_{pk} parameters each. This is commonly done by estimating the parameters for each copula density and then by choosing the one with either the lowest Aikaike information criterion (AIC) or the lowest Bayesian information criterion (BIC; Appendix B; Czado 2019), where
5.2. Iterating through tree levels
The methods above work for a single paircopula. It is straightforward to apply them to all paircopulas in the first tree level, but the same is not true in later tree levels. The reason is that, as is shown by Eq. (5), the estimate of a ddimensional PDF requires the conditional uniform variates u_{k−j} = F(x_{k}y_{−j}) and u_{j−j} = F(y_{j}y_{−j}), which, however, are not directly available.
To solve this issue, for the moment we suppose that the tree structure is known and the paircopulas up to the (ℓ−1)th tree level have been fit. In the ℓth tree, paircopulas have the form c_{i, jD}, where D is a set of ℓ − 1 variable indices called "conditioning set". Then there are always edges with indices (i, rD ∖ k) and (j, sD ∖ s) in the (ℓ−1)th tree.^{1} With the help of the so called hfunctions,
and
it can be shown that c_{i, jD}(., .) is the joint copula density of the random variables
and
Here, the point is that the arguments of the hfunctions have the same form of the corresponding conditional uniform variables, but the conditioning set has one index fewer. Therefore, it is possible to iterate the above equation until D = ∅, which corresponds to the first tree, where data are available. Because the paircopulas c_{i, rD \ r} and c_{j, sD \ s} have already been estimated, we can substitute the estimated models in the expressions above. In this way, we can transform data from paircopulas in one tree into data required for the estimation in the next tree. For example, we can express u_{123} = F(x_{1}x_{2}, x_{3}) required in Eq. (11) as
where
and u_{1} = F_{1}(x_{1}),u_{2} = F_{2}(x_{2}),u_{3} = F_{3}(x_{3}). The analytical form of the hfunctions is available for the most common copulas (Joe 1997; Schepsmeier & Stöber 2013).
5.3. Finding the tree structure
The remaining issue is how to select the right tree structure. For a Cvine, we need to specify which variable serves as the root node in every tree. For a Dvine, it is sufficient to specify the order of variables in the first tree. If d > 4, there are also structures other than D and Cvines.
To select an appropriate structure, the heuristic proposed by Dissmann et al. (2013) can be used. Their idea is to capture the strongest dependencies as early as possible in the tree structure. Here, “strength” is defined as the absolute value of Kendall’s τ (for the definition of this quantity, see Appendix C). We start in the first tree and compute the (empirical) pairwise Kendall’s τ for all variable pairs. Then, we choose the tree that maximizes the sum of absolute pairwise Kendall’s τ. We fit paircopula models for the edges and compute data for the next tree. On these data, we again compute the Kendall’s τ for all possible pairs and select the maximum spanning tree^{2}. We continue this way, iterating between structure selection, model fitting, and transforming the data until the whole model is fit. A summary of the whole procedure is given in Algorithm ^{1} and implemented in the VineCopula Rpackage (Nagler et al. 2019).
6. Application to an experimental set of data
6.1. Data set
We made use of the data published in Andreani et al. (2018) and complemented the molecular mass values with additional CO(10) line data taken at the NRO 45m antenna at Nobeyama (Andreani et al. 2020, and in prep.). The data set consists of the Kband luminosity, L_{K}, the infrared luminosity, L_{FIR}, the atomic hydrogen mass, M_{HI}, and the molecular mass, M_{H2}, derived from the CO(10) line luminosity toward the volumelimited local galaxy sample, the Herschel reference survey (HRS; Boselli et al. 2010). The data set is extensively described in Andreani et al. (2018) and references therein. Being volume limited, the sample contains all the galaxies above a given threshold of Kband luminosity, and the analysis would not be largely affected by a flux selection effect.
These variables were chosen because they are related to the main overall physical properties of the sample and their relation to the star formation activity in the galaxies. We aimed to investigate the relationship of those properties and derive insights into the driving physical mechanism in their interstellar medium.
6.2. Data analysis and interpretation
As the first step of the analysis, the PDF of each of the quantities logLK = log_{10}L_{K}, logLIR = log_{10}L_{FIR}, logMHI = log_{10}M_{HI} and logMH2v = log_{10}M_{H2} were modeled by means of the generalized lambda distribution (GLD) family (Karian & Dudewicz 2011, and references therein). The members of this family are fourparameter PDFs, which are known for their high flexibility and the large range of shapes that they can reproduce. The starship method has been adopted to fix the parameters. The reason is that this method finds the parameters that transform the data closest to the uniform distribution, which is an attractive characteristic when working with copulas. The results of the fit are shown in Fig. 3. After this step, the procedure presented in the previous section was applied with the random variates u, computed by means of the estimated margins.
Fig. 3.
Histograms of logLK, logLIR, logMHI, and logMH2v data. The red lines provide the PDF obtained by the fit of these data with the generalized lambda distribution family. 

Open with DEXTER 
The results are shown in Table 1 and Fig. 4. The original Kendall’s τ coefficients in Table 1 are related to the strengths of the relation between the quantities without being dependent on the derived margins. This shows that the strongest correlations occur between the farIR luminosity L_{FIR} and the gas masses, first molecular M_{H2} and then atomic M_{HI}, while L_{FIR} is weakly correlated with the nearIR Kband luminosity, L_{K}. On the other side, Fig. 4 indicates the type of vine structure selected, specifically a Cvine, and provides the list of paircopulas singled out for each edge. For each paircopula, the values of the corresponding coefficients and of the lower and upper tail dependence coefficients are also shown (for the meaning of last two quantities see Appendix D). The Kendall’s τ (for Tree 1) and partial Kendall’s τ (for Trees 2 and 3) associated to each edge are also shown. This last quantity measures the dependence between two variables after the effect of other variables (the common indices of two nodes) has been removed.^{3} In order to check the reliability of the obtained results, the procedure was repeated with the random variates u given by the pseudoobservations. As the comparison of Fig. 5 with Fig. 4 shows, the differences are not substantial. The fact that for the highest trees the copulas selected by the two methods are different is not significant. Indeed, one has to take present that when the Kendall’s τ between the random variates coming from two PDFs or two conditional PDFs is close to zero (i.e., they are almost uncorrelated), there are various kinds of copulas that can provide similar reconstructions of the corresponding bivariate data distribution. In other words, in the reconstruction of a multivariate distribution, the specific types of copula are only meaningful for values of the Kendall’s τ significantly different from zero.
Fig. 4.
Results concerning the application of the procedure described in Sects. 5 and 6.2 to the data corresponding to Fig. 3 with the random variates u computed by means of the estimated generalized Lambda PDFs (see text). Here, a copula is associated to each edge and Kendall’s τ for Tree 1, a partial Kendall’s τ for Trees 2 and 3 (column “tau”) and upper (column “utd”), respectively lower (column “ltd”) taildependence coefficients. These are theoretical quantities corresponding to the selected copulas of which the estimated parameters are given in the columns “par1” and “par2”. A description of these copulas can be found in Czado (2019). BB8_90 refers to copula BB8 rotated 90°. 

Open with DEXTER 
Fig. 5.
As in Fig. 4 but with random variates u computed by means of the pseudoobservations (see text). Clayton_90 and BB8_180 mean copula Clayton rotated 90° and copula BB8 rotated 180°, respectively. 

Open with DEXTER 
Sample Kendall’s τ.
These results can be more clearly interpreted by looking at Fig. 6. As explained in Sect. 5, the structure selection algorithm tries to capture the strongest dependencies first. Figure 6 shows a plot of the tree structure labeled with the Kendall’s τ (for Tree 1) and partial Kendall’s τ (for Trees 2 and 3). Here, the logLIR quantity as been selected as root node. This means that it is strongly correlated to all other variables and that it drives part of the dependence between the other variables. In the second tree, the effect of logLIR on the dependence between the others has been removed. There is only some weak negative dependence left.
Fig. 6.
Cvine structure selected by the procedure described in Sects. 5 and 6.2 to the data corresponding to Fig. 3. Each edge in Tree 1 is associated to a Kendall’s τ, whereas for Trees 2 and 3 they are associated to a partial Kendall’s τ. Here, 1→ logLK, 2→ logLIR, 3→ logMHI, 4→ logMH2v. 

Open with DEXTER 
All this can be interpreted with the fact that although from Table 1 the quantities M_{HI} and M_{H2} appear positively dependent, such dependence appears to be driven entirely by their dependence on the quantity L_{FIR}. Once the dependence of L_{FIR} is removed from the relation with the other quantities the residual relations M_{HI} with M_{H2} and L_{K} with M_{HI} are negatively dependent (albeit this dependence is quite weak). This means that the dependence shown in Table 1 is driven entirely by their dependence on L_{FIR}.
Since L_{FIR} is dominated by the thermal dust emission heated by FUV photons by massive stars and residing in molecular clouds, which are cocoons of star formation processes, this result confirms that the physical properties of the galaxies are driven by their star formation.
For completeness, in Fig. 7 we show the original data versus the data simulated from the estimated 4D joint PDF of which the 2D slices are shown in Fig. 8. Figure 7 shows a good agreement between original and simulated data, while Fig. 8 demonstrates the onetoone relation between the couple of variables.
Fig. 7.
Original logLK, logLIR, logMHI, and logMH2v data (blue circles) versus the corresponding simulated data obtained from the estimated 4D joint PDF (red circles). 

Open with DEXTER 
Fig. 8.
Twodimensional slices of the estimated 4D joint PDF for the data corresponding to Fig. 3. Along the diagonal are the fit PDFs of Fig. 3. The right panels show the slices in which colors correspond to the intensity of the relation, while the left panels report, on the same slices, the data points and the isocontours. 

Open with DEXTER 
7. Conclusions
In this work, a flexible and effective approach to modeling the relationship of a set of experimental multidimensional quantities is presented. This approach consists of modeling the joint PDF of the data by means of a special type of copula called a vine copula. Classical copulas are functions that contain all of the information on the relationship between two random quantities. Their major limitation is that they are unable to model multidimensional data. Vine copulas overcome this limitation by expressing the joint PDFs as the product of a set of 2D copula densities and the 1D PDFs corresponding to each quantity. In particular, two types of vine copulas have been considered: the Cvine and the Dvine. This approach makes the estimation of the joint PDFs amenable to a theoretical treatment and feasible from the computational point of view.
We applied this method to published data on the nearIR and farIR luminosities and atomic and molecular masses of the HRS. We find that the farIR luminosity, L_{FIR}, is the key player in driving the galaxy properties in this sample. Despite its original selection in the Kband, the HRS sample shows that it is L_{FIR} that plays a fundamental role. Removing its dependence from the other variables, the Kband luminosity, and the atomic and molecular masses, makes it clear that the established relation among these quantities does not show up any more.
The L_{FIR} in this sample is dominated by the thermal dust emission heated by FUV photons produced by massive stars in molecular clouds. Our analysis therefore highlightsthat the star formation activity of these galaxies is the key parameter driving the galaxy evolution.
References
 Aas, K., Czado, C., Frigessi, A., & Bakken, H. 2009, Insurance: Math. Econ., 44, 182 [CrossRef] [Google Scholar]
 Andreani, P., Spinoglio, L., Boselli, A., et al. 2014, A&A, 566, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Andreani, P., Boselli, A., Ciesla, L., et al. 2018, A&A, 617, A33 [CrossRef] [EDP Sciences] [Google Scholar]
 Andreani, P., Miyamoto, Y., Kaneko, H., et al. 2020, A&A, submitted [arXiv:2008.11134] [Google Scholar]
 Balakrishnan, N., & Lai, C. D. 2010, Continuous Bivariate Distributions (New York: Springer) [Google Scholar]
 Boselli, A., Eales, S., Cortese, L., et al. 2010, PASP, 122, 261 [NASA ADS] [CrossRef] [Google Scholar]
 Burnham, K. P., & Anderson, D. R. 2002, Model Selection and Multimodel Inference (New York: Springer) [Google Scholar]
 Czado, C. 2019, Analyzing Dependent Data with Vine Copulas (New York: Springer) [CrossRef] [Google Scholar]
 Dissmann, J., Brechmann, E. C., Czado, C., & Kurowicka, D. 2013, Comput. Stat. Data Anal., 59, 52 [CrossRef] [Google Scholar]
 Durante, F., & Sempi, C. 2016, Principles of Copula Theory (New York: CRC Press) [Google Scholar]
 Hofert, M., Kojadinovic, I., Mächler, M., & Yan, J. 2018, Elements of Copula Modeling with R (New York: Springer) [CrossRef] [Google Scholar]
 Joe, H. 1997, Multivariate Models and Dependence Concepts (Dordrecht: SpringerScience+Business Media) [CrossRef] [Google Scholar]
 Joe, H. 2015, Dependence Modeling with Copulas (New York: CRC Press) [Google Scholar]
 Karian, Z. A., & Dudewicz, E. J. 2011, Handbook of Fitting Statistical Distributions with R (New York: CRC Press) [Google Scholar]
 Kotz, S., Balakrishnan, N., & Johnson, N. L. 2000, Continuous Multivariate Distributions (New York: John Wiley & Sons, Inc.), 1 [CrossRef] [Google Scholar]
 Lin, C. A., & Kilbinger, M. 2015, A&A, 583, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lin, C. A., Kilbinger, M., & Pires, S. 2016, A&A, 593, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nagler, T., & Czado, C. 2016, J. Multivariate Anal., 151, 69 [CrossRef] [Google Scholar]
 Nagler, T., Schellhase, C., & Czado, C. 2017, Depend. Model., 5, 99 [CrossRef] [Google Scholar]
 Nagler, T., Schepsmeier, U., Stoeber, J., et al. 2019, VineCopula: Statistical Inference of Vine Copulas. R package version 2.3.0, https://CRAN.Rproject.org/package=VineCopula [Google Scholar]
 Nelsen, R. B. 2006, An Introduction to Copulas (New York: Springer Science + BusinessMedia, Inc.) [Google Scholar]
 Sato, M., Ichiki, K., & Takeuchi, T. 2010, Phys. Rev. Lett., 105, 251301 [NASA ADS] [CrossRef] [Google Scholar]
 Sato, M., Ichiki, K., & Takeuchi, T. 2011, Phys. Rev. D, 83, 023501 [NASA ADS] [CrossRef] [Google Scholar]
 Schepsmeier, U., & Stöber, J. 2013, Stat. Pap., 55, 525 [CrossRef] [Google Scholar]
 Scherrer, R. J., Berlind, A. A., Mao, Q., & McBride, C. K. 2010, AJ, 708, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuchi, T. T. 2010, MNRAS, 406, 1830 [NASA ADS] [Google Scholar]
 Takeuchi, T. T., Sakurai, A., Yuan, F. T., & Burgarella, D. 2011, Earth Planet Space, 65, 281 [CrossRef] [Google Scholar]
 Vio, R., Fasano, G., Lazzarin, M., & Lessi, O. 1994, A&A, 289, 640 [NASA ADS] [Google Scholar]
Appendix A: The Johnson’s distribution and the generalized Lambda distribution families
The Johnson’s distributions and the generalized Lambda distributions (GLD) are both fourparameter families that are used for fitting distributions to a wide variety of data sets. In particular, the Johnson’s system is based on three different PDFs, f_{U}(x), f_{B}(x), and f_{L}(x) according to the fact that the random variable x is unbounded, bounded both above and below, and bounded only below:
for −∞< x < ∞,
for ϵ ≤ x ≤ λ + ϵ, and
for x ≥ ϵ. In literature, various methods are available for the selection of the appropriate type of PDF as well for the estimate of the parameter (e.g., Vio et al. 1994; Karian & Dudewicz 2011).
The PDFs corresponding to the GLD family are given by
where x = Q(y; λ_{1}, λ_{2}, λ_{3}, λ_{4}) with
and 0 ≤ y ≤ 1. Concerning this family, various methods are also available for the estimate of the parameters (e.g., Karian & Dudewicz 2011).
Appendix B: AIC and BIC criteria
The Akaike information criterion (AIC) and the Bayesian information criterion (BIC) are two criteria for model selection from a finite set of models (Burnham & Anderson 2002). They are based on the maximum value of the likelihood function for the model as well on the number n_{p} of free parameters it contains. The idea is that, when fitting models, it is possible to increase the likelihood by adding parameters, but doing so may result in overfitting. Both the BIC and AIC attempt to resolve this problem by introducing a penalty term for the number of parameters in the model. In particular,
whereas
with n being the number of data. In practical applications, a set of models is chosen, the corresponding quantity evaluated, Eqs. (B.1) or (B.2) used, and finally the model with the lowest AIC or BIC selected. The difference between AIC and BIC is how much model complexity (i.e., the number of parameters) is penalized. For n ≥ 8, the BIC penalty is stronger. Both criteria give a mathematical guarantee to find the “best” model as the sample size increases. The BIC assumes that the true model is among the set of candidates, but the AIC does not. These criteria are useful in the context of the vine copulas since the different types of bivariate copulas considered for their construction contain a different number of free parameters.
Appendix C: The Kendall’s τ
When working with copulas the relationship between two random quantities is typically measured by means of the Kendall’s τ. The reason can be understood by looking at Fig. C.1, which shows the realization of 1000 independent copies of a bivariate random vector (x_{1}, x_{2}) from the Gaussian, exponential and Cauchy PDFs and of the same number of a bivariate random vector (u_{1}, u_{2}) from the uniform PDF. These realizations appear quite different from one another, as well as the corresponding linear correlation coefficients ρ. Here, the point is that the first three sets of random numbers {(x_{1, i}, x_{2, i})} were obtained from the set of uniform random pairs {(u_{1, i}, u_{2, i})} by means of the transformations:
Fig. C.1.
Numerical realization of 1000 independent copies of a bivariate random vector (x_{1}, x_{2}) from the Gaussian, exponential, and Cauchy PDFs obtained from the set of uniform random pairs {(u_{1, i}, u_{2, i})}, shown in the bottomright panel, by means of the transformations (x_{1}, x_{2}) = (F^{−1}(u_{1}),F^{−1}(u_{2})) where F^{−1}(u) is the inverse CDF corresponding to the various PDFs. 

Open with DEXTER 
where F^{−1}(u) is the inverse CDF corresponding to the various PDFs. This is a common method to simulate random numbers from a given PDF. What this figure indicates is that the different appearance of the realizations is not due to the intrinsic relationship between the random quantities, but rather to their margins. Since with copulas one wants to disentangle margins from the dependence structure, the latter should be measured in a way that does not depend on the marginal distributions. This is what the Kendall’s τ does.
If is and independent copy of (x_{1}, x_{2}), τ is defined as
that is, it is the probability of concordance minus the probability of discordance of the random pairs (x_{1}, x_{2}) and . The rationale behind this definition is that if there is positive dependence between the variable x_{1} and x_{2}, then when x_{1} increases or decreases, a similar behavior has to be expected for x_{2}. It can be demonstrated (Hofert et al. 2018) that
meaning that τ effectively depends only on the underlying copula.
The sample version of τ is given by
where n is the number of observations and sign[x] = 1 if x > 0, sign[x] = 0 if x = 0 and sign[x] = − 1 if x < 0. As expected, is the same for all the realizations in Fig. C.1.
Appendix D: Taildependence coefficients
There are situations where in the 2D scatterplot of a set of data, the points appear concentrated in one or both the tails of their joint distribution. For instance, this is the case for the scatterplots in Fig. C.1, where a concentration of points in the lowerleft tail of the joint distribution is evident. Joint distributions characterized by welldeveloped tails indicate a high probability of joint occurrence of extremely small and/or large values. In some practical applications, it is useful to have an estimate of this probability. Given the margins F_{1}(x_{1}) and F_{2}(x_{2}) and the copula C(u_{1}, u_{2}), the coefficients of lower and upper tail dependence provide such an estimate and are defined as
respectively,
These coefficients are conditional probabilities that measure the tendency of the random variable x_{2} to behave as the random variable x_{1}. When their value is close to one, it means tail dependence (i.e., high probability of joint extreme values), when close to zero it means tail independence (i.e., low probability of joint extreme values). The analytical expression of λ_{l} and λ_{u} is available for various parametric copulas. For instance, the random points in the bottomleft panel of Fig. C.1 has been generated through a Clayton copula with coefficient θ = 5, for which λ_{l} = 0.87 and λ_{u} = 0. These values also hold for the other distributions in the same figure.
All Tables
All Figures
Fig. 1.
Example of tree structure of a 4D Dvine copula (see text). 

Open with DEXTER  
In the text 
Fig. 2.
Example of tree structure of a 4D Cvine copula (see text). 

Open with DEXTER  
In the text 
Fig. 3.
Histograms of logLK, logLIR, logMHI, and logMH2v data. The red lines provide the PDF obtained by the fit of these data with the generalized lambda distribution family. 

Open with DEXTER  
In the text 
Fig. 4.
Results concerning the application of the procedure described in Sects. 5 and 6.2 to the data corresponding to Fig. 3 with the random variates u computed by means of the estimated generalized Lambda PDFs (see text). Here, a copula is associated to each edge and Kendall’s τ for Tree 1, a partial Kendall’s τ for Trees 2 and 3 (column “tau”) and upper (column “utd”), respectively lower (column “ltd”) taildependence coefficients. These are theoretical quantities corresponding to the selected copulas of which the estimated parameters are given in the columns “par1” and “par2”. A description of these copulas can be found in Czado (2019). BB8_90 refers to copula BB8 rotated 90°. 

Open with DEXTER  
In the text 
Fig. 5.
As in Fig. 4 but with random variates u computed by means of the pseudoobservations (see text). Clayton_90 and BB8_180 mean copula Clayton rotated 90° and copula BB8 rotated 180°, respectively. 

Open with DEXTER  
In the text 
Fig. 6.
Cvine structure selected by the procedure described in Sects. 5 and 6.2 to the data corresponding to Fig. 3. Each edge in Tree 1 is associated to a Kendall’s τ, whereas for Trees 2 and 3 they are associated to a partial Kendall’s τ. Here, 1→ logLK, 2→ logLIR, 3→ logMHI, 4→ logMH2v. 

Open with DEXTER  
In the text 
Fig. 7.
Original logLK, logLIR, logMHI, and logMH2v data (blue circles) versus the corresponding simulated data obtained from the estimated 4D joint PDF (red circles). 

Open with DEXTER  
In the text 
Fig. 8.
Twodimensional slices of the estimated 4D joint PDF for the data corresponding to Fig. 3. Along the diagonal are the fit PDFs of Fig. 3. The right panels show the slices in which colors correspond to the intensity of the relation, while the left panels report, on the same slices, the data points and the isocontours. 

Open with DEXTER  
In the text 
Fig. C.1.
Numerical realization of 1000 independent copies of a bivariate random vector (x_{1}, x_{2}) from the Gaussian, exponential, and Cauchy PDFs obtained from the set of uniform random pairs {(u_{1, i}, u_{2, i})}, shown in the bottomright panel, by means of the transformations (x_{1}, x_{2}) = (F^{−1}(u_{1}),F^{−1}(u_{2})) where F^{−1}(u) is the inverse CDF corresponding to the various PDFs. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.