The need for a method to construct multidimensional distribution function is increasing recently, in the era of huge multiwavelength surveys. We have proposed a systematic method to build a bivariate luminosity or mass function of galaxies by using a copula. It allows us to construct a distribution function when only its marginal distributions are known, and we have to estimate the dependence structure from data. A typical example is the situation that we have univariate luminosity functions at some wavelengths for a survey, but the joint distribution is unknown. Main limitation of the copula method is that it is not easy to extend a joint function to higher dimensions ($d > 2$), except some special cases like multidimensional Gaussian. Even if we find such a multivariate analytic function in some fortunate case, it would often be inflexible and impractical. In this work, we show a systematic method to extend the copula method to unlimitedly higher dimensions by a vine copula. This is based on the pair-copula decomposition of a general multivariate distribution. We show how the vine copula construction is flexible and extendable. We also present an example of the construction of an stellar mass--atomic gas--molecular gas 3-dimensional mass function. We demonstrate the maximum likelihood estimation of the best functional form for this function, as well as a proper model selection via vine copula.

(abridged) The scale of the acoustic oscillation of baryons at the baryon-photon decoupling is imprinted on the spatial distribution of galaxies in the Universe, known as the baryon acoustic oscillation (BAO). The correlation functions and power spectrum are used as a central tool for the studies on the BAO analysis. In this work, we analyzed the spatial distribution of galaxies with a method from the topological data analysis (TDA), in order to detect and examine the BAO signal in the galaxy distribution. The TDA provides a method to treat various types of "holes" in point set data, by constructing the persistent homology (PH) group from the geometric structure of data points and handling the topological information of the dataset. We can obtain the information on the size, position, and statistical significance of the holes in the data. A particularly strong point of the persistent homology is that it can classify the holes by their spatial dimension, i.e., a 0-dim separation, 1-dim loop, 2-dim shell, etc. We first analyzed the simulation datasets with and without the baryon physics to examine the performance of the PH method. We found that the PH is indeed able to detect the BAO signal: simulation data with baryon physics present a prominent signal from the BAO, while data without baryon physics does not show this signal. Then, we applied the PH to a quasar sample at $z <1.0$ from extended Baryon Oscillation Spectroscopic Survey in Sloan Digital Sky Survey Data Release 14. We discovered a characteristic hole (a hollow shell) at a scaler $\sim150\ [{\rm Mpc}]$. This exactly corresponds to the BAO signature imprinted in the galaxy/quasar distribution. We performed this analysis on a small subsample of 2000 quasars. This clearly demonstrates that the PH analysis is very efficient in finding this type of topological structures even if the sampling is very sparse.