Dependency Clustering of Mixed Data with Gaussian Mixture Copulas

1 Introduction

The fundamental task of clustering has been studied extensively in machine learning [Jain, 2010]. Real world data often contains mixed - continuous and discrete - features as well as complex dependencies between features. Dependency clustering, that detects clusters corresponding to feature dependency patterns, for mixed data remains a challenging open problem [Plant, 2012; Wang et al., 2015].

Clustering methods for mixed data, that are much fewer than those for continuous data, often do not address dependencies satisfactorily. For example, K-prototypes [Huang, 1998] quantifies similarity between mixed observations but does not address dependency; K-Means-Mixed [Ahmad and Dey, 2007] models dependencies of discrete features only; mixture [Hunt and Jorgensen, 2011], assuming normally distributed data, captures dependence of continuous features only. A coupled representation (which can subsequently be used for clustering) for nominal data [Wang et al., 2011] and for continuous data [Wang et al., 2013] has recently been combined for mixed data [Wang et al., 2015]. The method, called CoupledMC, transforms mixed to continuous valued data through a series of steps that capture dependencies both within and between continuous and discrete features using non-parametric metrics like correlations and co-occurrence. The high computational complexity (upto cubic in dimension and quadratic in size of data, for CoupledMC) is a disadvantage of many of these distance-based methods.

Correlation clustering methods discover clusters in subspaces based on correlations revealed by low dimensional representations, e.g. through Principal Component Analysis. This approach has been extended to mixed data in INCONCO [Plant and B¨ohm, 2011], SCENIC [Plant, 2012] and SpectralCAT [David and Averbuch, 2012] that find low dimensional embeddings of the data, in different ways, to detect clusters. INCONCO models dependencies by distinct Gaussian distributions for each category of each discrete feature. While SCENIC is not as restrictive in the dependencies, it also assumes a Gaussian distribution to find the embedding space. SpectralCAT discretizes continuous features before spectral clustering using an adaptive Gaussian kernel. Normality assumptions limit the modeling capability of these techniques.

Model-based approaches for continuous data can leverage the flexible framework of copulas that provides a modular parameterization of multivariate distributions - arbitrary marginals independent of dependency models from copula families which can model a wide variety of linear and nonlinear dependencies. For example, with the Gaussian copula itself, using different marginals, many different joint distributions can be constructed, called meta-Gaussian distributions, that have been used in several applications (see [Elidan, 2013]) including multi-view clustering [Rey and Roth, 2012]. However, meta-Gaussian dependencies from the Gaussian copula cannot model many kinds of dependencies, notably asymmetric and tail dependencies, and the Gaussian mixture copula (GMCM) was proposed by [Tewari et al., 2011] to allow flexible dependency modeling going beyond metaGaussian distributions. These copula-based models can be used with continuous features only.

With discrete data, copula dependencies are not marginalfree but they can still be used effectively [Genest and Neslehova, 2007]. Clustering with a mixture of copulas was proposed by [Kosmidis and Karlis, 2015] but for discrete marginals it is computationally inefficient, the runtime for parameter estimation being exponential in data dimensions.

An efficient alternative not yet explored for dependency clustering, is the semiparametric extended rank likelihood (ERL) approach of [Hoff, 2007] that allows us to use the modeling flexibility of copulas with mixed data. We exploit this technique to enable the efficient use of GMCM on mixed data, to model complex dependencies, that in turn, improves clustering performance. We analyze, both theoretically and empirically, and find that the ERL approximation yields accurate parameter estimates, and good clustering performance, when all features are continuous or ordinal with not many levels.

To summarize, our contributions are:

1. We present a Gaussian mixture copula based clustering algorithm that:
2. canmodel complex dependencies beyond those captured by meta-Gaussian distributions,
3. can fit mixed - continuous, ordinal and binary valued - data, which previous GMCM-based methods cannot and
4. scales linearly with size and quadratically with dimensions of input, which is significantly faster than state-ofthe-art correlation clustering methods for mixed data.
5. We theoretically analyze the conditions necessary to obtain consistent and asymptotically unbiased estimates of GMCM parameters, which, to our knowledge, is the first such analysis of an ERL approach for mixed data clustering.
6. Our experimental results demonstrate the efficacy of our method, that outperforms state-of-the-art methods for correlation clustering on synthetic and real benchmark data sets with mixed features, thus illustrating the advantage of our copula-based approach for dependency clustering.

2 GaussianMixture Copulas

A p-dimensional copula is a multivariate distribution function C : [0↪ 1]p 7! [0↪ 1]. A theorem by [Sklar, 1959] proves that copulas can uniquely characterize continuous joint distributions: for every joint distribution with continuous marginals, F(Y1↪ ▷ ▷ ▷ ↪ Yp), there exists a unique copula function such that F(Y1↪ ▷ ▷ ▷ ↪ Yp) = C(F1(Y1)↪ ▷ ▷ ▷ ↪ Fp(Yp)) as well as the converse. In the discrete case, the copula is uniquely determined on Ran(F1) ⇥ ▷ ▷ ▷ ⇥ Ran(Fp), where Ran(Fj) is the range of marginal Fj. See [Joe, 2014] for a comprehensive treatment of copulas. Parametric copula families are typically defined on uniform random variables obtained through CDF transformations from the marginals. In a Gaussian Mixture Copula Model (GMCM) [Tewari et al., 2011], the dependence is obtained from a Gaussian Mixture (GMM). Note the difference from a mixture of copulas [Kosmidis and Karlis, 2015]. GMCM copula density is given by:

The generative process is:

zij ⇠ p-dimensional↪ G-component GMM, parameters: #

Notation: The observed data matrix is Y = ¶yij♦; indices i↪ j denote the observation (i = 1↪ ▷ ▷ ▷ ↪ n) and dimension (j = 1↪ ▷ ▷ ▷ ↪ p); Fj is the unknown marginal distribution for the j th dimension, and uniform random variables uij = Fj(yij). Latent variables Z 2 Rn⇥p = ¶zij♦, zi denotes p-dimensional row vectors and zj denotes ndimensional column vector along the j th dimension; j and j are the marginal CDF and density of the GMM along the j th dimension, φ is the multivariate Gaussian density and # = ¶⇡g↪θg↪⌃g♦ is the (unknown) parameter set representing mixing proportions (⇡g > 0, with PGg=1 ⇡g = 1), mean vectors (θg) and covariance matrices (⌃g; σgij : ij th entry of ⌃g) for GMM component g = 1↪ ▷ ▷ ▷ ↪ G. Vectors and matrices are in bold font.

Copula Parameter Estimation. Standard maximum likelihood inference requires specifying a parametric family for each marginal and inferring the parameters simultaneously with the copula parameters which is computationally prohibitive for even moderate dimensions. The two-step IFM procedure [Joe, 2014] requires fitting a marginal to each feature and then using the CDF transformation to obtain pseudo data for estimating copula parameters. Often marginal families are unknown and a semiparametric approach is to use rank-transformed scaled empirical marginals to obtain the pseudo data. The copula function for a joint distribution is invariant to monotone increasing transformations of marginals and the resulting pseudo-likelihood estimator is consistent for continuous marginals [Genest et al., 1995].

With discrete marginals, rank transformation leads to ties with non-zero probability. In general, technical hurdles have to be surmounted to use copulas with discrete marginals: the copula may be unidentifiable and may not be margin-free, inference algorithms based on Kendall's tau or Spearman's rho may be biased or inconsistent [Genest and Neslehova, 2007]. The Extended Rank Likelihood (ERL) approach is an approximation to the full likelihood that can accommodate both continuous and discrete (ordinal and binary) data [Hoff, 2007].

The key idea of ERL is that, since F -1 j is monotonic, although we don't observe variables zi directly, we still know the order of zi's induced by the data. Estimation of the copula parameters # using rank likelihood can be implemented by conditioning on the partial ordering induced by Y. The ERL function is given by P(Z 2 D♣#); D is the rank-induced set: D = ¶Z 2 Rn⇥p : zlj < zij < zuj♦ where zlj = max¶zkj : ykj < yij♦ and zuj = min¶zkj : yij < ykj♦, along the j th dimension. Since the ERL approach bypasses the estimation of marginals, it can be applied to both discrete and continuous data. Note that the ERL function is equivalent to the distribution of the ranks for continuous data. For discrete data, the distribution of the ranks depends on the marginals, which is ignored resulting in some loss of information. Bayesian inference can be achieved by constructing a Markov chain having a stationary distribution equal to P (#♣Z 2 D) / P (#) ⇥ P (Z 2 D♣#).

3 Clustering with GaussianMixture Copulas

Our algorithm, EGMCM (Extended GMCM), fits a Gaussian mixture copula to data in a semiparametric manner and infers clusters from the obtained parameters. Algorithm 1 outlines our iterative scheme for obtaining approximate samples to estimate the posterior P(#♣Z 2 D). Each iteration has 2 steps:

estimating GMM parameters # and resampling Z.

The algorithm is initialized with two steps. First, rank transformed data is obtained by transforming each data point yij of the feature vector Yj to R(yij) = 1n Pni=1 1(Yij  yij), the scaled rank of yij in all the observations of Yj, the j th feature (where 1(c) denotes the indicator function which is 1 if the condition c is true and 0 otherwise). To initialize the values of Z, we take the inverse CDF of a standard normal on the rank transformed data. We empirically find that this simple choice works well; note that the values of Z are resampled through the iterations.

We use Expectation Maximization (EM) to estimate #. Gibb's sampling could also be used. If a single Gibb's sweep is used to sample # instead of a full estimation, the clustering performance deteriorates - also explained by our theoretical analysis which shows that initial unbiased estimates of # are important to obtain final unbiased mean estimates.

In the Resample Z step, we sample from each Gaussian in the GMM, truncated by the bounds zuj and zlj which constrains the sampled z to be in the set D. The Truncated Normal (TN) random variates are then combined using the mixing proportions of the estimated GMM. The function unique, that returns the set of unique values in the input, is used to avoid repeated computation for ties.

Algorithm 1 EGMCM

Input: R(Y), scaled rank transformed data Y; G, number

of clusters Initialization Z = Φ-1 (R(Y)) loop Estimate, via EM, GMM parameters # = [⇡g↪ θg↪ ⌃g] Resample Z: for j: 1 to p do for all y 2 unique ¶y1j↪ ▷ ▷ ▷ ↪ ynj♦ do Compute zlj = max¶zij : yij < y♦ and zuj = min¶zij : y < yij♦ For each i such that yij = y: Sample rgij from TN(θgj↪ σgij↪ zlj↪ zuj) Set zij = PGg=1 ⇡grgij end for end for end loop Output: Cluster labels (latent variables of GMM (Z♣#))

The algorithm can easily handle data missing-at-random. If yij is missing, zij is obtained (in each Resample Z step) from the unconstrained GMM with current estimates of #.

Computational Complexity Comparison

The time complexity of EGMCM is O(tGnp2 ) where t is the number of iterations of the outer loop. The runtime is dominated by EM (O(Gnp2 )). Empirically we obtain good results with t < 50 and so, the runtime is comparable to that of GMM estimation algorithms.

In comparison, note that the runtime for SCENIC is O(t0n2p2v ) where t0 includes terms proportional to number of iterations in their algorithm and pv is the dimension of the low-dimensional embedding; the runtime for coupledMC is O(n2R3 ) where R is the maximal number of attribute values for all features, including continuous features that are discretized in their method, hence p  R  np.

4 Theoretical Analysis

Our iterative algorithm has two steps: estimation of # through EM and the second Resample Z step. EM estimates of GMM parameters are consistent and asymptotically unbiased when the component means are not too close to each other; and for large sample sizes (> 200) the bias in mean and variance can be ignored for practical purposes [Nityasuddhi and Bohning, 2003]. Thus, assuming that the initial estimates obtained via EM in the first iteration of the loop are consistent and asymptotically unbiased, for large number of iterations m, we analyze the deviation in the mean from the initial estimate in subsequent Resample Z steps for discrete (theorem 2) and continuous (theorem 3) marginals. We show that the deviation tends to zero under simple conditions that can be maintained in the algorithm.

Similar properties cannot be inferred for the variance parameters because we are assuming a full variance-covariance matrix; thus the behaviour of the variance components after truncating the GMM will be dependent on the data, the nature of truncation and position of the component means with respect to the mixture means, among other factors. To our knowledge, similar analysis of previous ERL based methods has not been done before.

Notation: Superscript (m) denotes the iteration; zl and zu depend on g and j, i.e. the corresponding cluster and the dimension but we do not use them for notational simplicity. Let σ2 gj denote the j th diagonal element of ⌃g, and δgj denote the contribution of the component mean of cluster g to the mixture mean along the j th dimension, i.e.

λ(m) gj is the deviance in mean due to truncation of the gth component of the normal distribution along the j th dimension in the mth iteration (which depends on the coefficient of variation, σθ ). Note that the mixture mean of the GMM is θ = PGg=1 ⇡gθg and the mixture variance is ⌃ = PGg=1 ⇡g h(θg -θ) (θg -θ)T +⌃g i. In addition to set D defined in section 2, we define the set D1 = ¶Z : zl1  zij  zu1♦ where zl1 = max¶zkj : ykj = yij♦, zu1 = min¶zkj : ykj = yij♦. The initial estimate of mean, θ(0) is obtained in the first iteration of the outer loop, through EM, and is considered asymptotically unbiased.

Lemma 1 After m iterations of the loop in algorithm EGMCM, the deviation of the g-th component of mixture mean θ(m) from the g-th component of the initial estimate θ(0) along the j th dimension is θ(0) gj hQmk=1 (1-λ(k-1) gj ) -1i.

Proof: According to the algorithm, for each j 2 ¶1↪ 2↪ ▷▷▷↪ p♦, for each i 2 ¶1↪ 2↪ ▷▷▷↪ n♦ and for each y 2 unique ¶y1j ↪ y2j ↪ ▷▷▷↪ ynj♦, at the mth iteration we sample rgij from the truncated normal distribution N(θ(m-1) gj ↪ σ(m-1) gj 2 ) with the truncation points (zl↪ zu), as described in the algorithm, for g = 1↪ 2↪ ▷▷▷↪G. Then zij is taken as PGg=1 ⇡grgij and θ is estimated based on zij. But since rgij is taken from a truncated normal distribution,

Repeating the iteration and applying the above equation, after m iterations we have θ(m) gj = θ(0) gj Qmk=1 (1-λ(k-1) gj )▷

Theorem 1 If the marginal distribution, Fj, is discrete and there exists some m0 such that (zl↪ zu) does not contain 0 at the mth 0 iteration, and if (z(m) u -z(m) l )◁σ(m) gj ! 0 asm!1 then λ(m) gj ! 0 as m!1.

Proof: For discrete Fj, if z(m) l and z(m) u are the two extremes at the mth iteration, then z(m) u - z(m) l  z(m-1) u - z(m-1) l for all m.

What if equality holds for all subsequent iterations? λ will not tend to 0 then. This is because of the discreteness of the marginals. Recall, for the extended rank likelihood, we consider the set D and not the set D1. For discrete Fj, the measure of D1 is non-zero. Thus (zD l ↪ zD u ) ⇢ (zD[D1 l ↪ zD[D1 u ) where zD l is the lowest order statistic taken based on the set D only, and zD[D1 l is the lowest order statistic taken based on the union of D and D1 (similarly for zu). Thus at each iteration some positive measure will be assigned to the set D1 and hence z(m) u - z(m) l < z(m-1) u - z(m-1) l , i.e. the relation will actually be a strict inequality.

Let us consider the sequence sm = z(m) u -z(m) l . From the discussion above, we can infer that sm is a positive sequence and sm < sm-1. Thus sm ! 0 as m ! 1. In other words, for each ¯m there exists an ✏ such that sm < ✏ for all m ≥ ¯m. Thus z(m) u - z(m) l ! 0.

Now, since θ(m) gj 2 ⇣z(m) l ↪ z(m) u ⌘, (z(m) u - z(m) l )◁σ(m) gj ! 0 implies that (z(m) u - θ(m) gj )◁σ(m) gj ! 0 and (z(m) l - θ(m) gj )◁σ(m) gj ! 0 as m ! 1. Thus, using the expression of λ(m) gj we get

where δ = (z(m) u -θ(m) gj )◁σ(m) gj . Using L'Hospital's rule,

Theorem 2 If the marginal distribution, Fj, is continuous and the initial values are sampled unrestrictedly (i.e. from the entire sample space), the deviation in the estimate of θ along the j th dimension after the mth iteration, λ(m) gj ! 0 as m ! 1.

Thus we show that zl and zu become closer to each other, and the change in difference of their respective CDFs also becomes smaller (since the Gaussian CDF is a uniformly continuous function of its arguments). Also, as zl and zu becomes closer and the subsequent iterations are based on samples taken from the shrinking set (zl↪ zu), σ(m) gj tends to 0. Also, since there exists m0 such that (zl↪ zu) does not contain 0 at the mth 0 iteration, it will not contain 0 for all the subsequent iterations since (zl↪ zu) is a shrinking set. Thus θ(m) gj 6! 0 for all m ≥ m0. Using these conditions in equation 2, we have λ(m) gj ! 0 as m!1.

Proof: For continuous Fj , we have ¶Z♦ = D [ D1 = D (since the measure of D1 is zero). Thus zl and zu will not change over the iterations, if the initial values are taken over the entire domain space of z along the j th dimension. Thus, after the mth iteration, as m!1,

and

Similarly σ(m) gj - σ(m-1) gj ! 0 and θ(m) gj - θ(m-1) gj ! 0 as m!1. Thus after iteration m, λ(m) gj ! 0 as m!1.

Summary of Results. The deviation (due to repeated Resample Z steps) from an initial asymptotically unbiased estimate and subsequent loss of information depends on the initial values of Z. For a continuous feature, the deviation along the dimension is less if the initial values range over the entire domain of Z. More interestingly, it turns out that for the discrete case the deviation is less if the highest and lowest order statistics of Z along the particular dimension are close to each other, that is, when the number of levels in the discrete feature is less. The analysis implies that EGMCMworks well when all features are continuous or when there is mixed data with ordinal features that do not have many levels. In the presence of discrete valued features, the number of iterations should be increased to obtain better estimates.

5 Experiments

We compare the performance of our algorithm, EGMCM, with that of SCENIC and CoupledMC. Previous experiments (in [Plant, 2012] and [Wang et al., 2015] respectively) have shown that SCENIC outperforms K-Means-Mixed and INCONCO; and CoupledMC outperforms K-Means-Mixed, KPrototype, mixture and SpectralCAT. In addition, we use GMM, with 'binarized' categorical features as a baseline.

PerformanceMetrics. To evaluate the performance of clustering algorithms we use Adjusted Rand Index (ARI) [Hubert and Arabie, 1985] (range [-1↪ 1]) and Adjusted Mutual Information (AMI) [Vinh et al., 2010] (range [0↪ 1]), higher values indicate better clustering in both metrics.

5.1 Simulations

We simulate both numerical, with only continuous valued features, and mixed datasets with dependencies.

Numerical Data. The number of clusters is chosen to be G = 2. The ratio of the number of the two cluster sizes is 1:1. Each data point in the first cluster is a product of a sample from a Multivariate Normal distribution (MVN(0↪ Ip◁2) where Ip is a p ⇥ p Identity matrix) and a sample from a F distribution (F(df1 = p◁2↪ df2 = p◁2)). Each data point in the second cluster is a product of a sample from a Multivariate Normal distribution (MVN(2↪∆) where ∆ is a matrix with (i↪ j)th element 0▷9♣i-j♣ ) and a sample from a Uniform distribution (Unif(0↪ 1)). We use 10 simulation settings and simulate 25 datasets for each setting as shown in table 1.

Mixed Data. We first generate two clusters, each with n◁2 p◁2-dimensional observations with discrete data. This data is concatenated with p◁2-dimensional numerical data (as above) to obtain an n ⇥ p dataset. Discrete data in the two clusters is sampled from 4 different 20-category multinomial distributions Pi (i = 1↪ 2↪ 3↪ 4) that are generated in such way that in each distribution one or two of the 20 classes are assigned a large probability, and the rest are assigned very small probability. The large probability class(es) vary from one distribution to another. Let xkij (yk ij ) denote the i-th observation of the j-th feature in the k-th cluster of the numerical (discrete) data. Dependencies between the numerical and discrete data are imposed by the following two rules. If x1ij < 0 then sample y1 ij fromP1, else fromP2. If x2ij < 0 then sample y2 ij from P3, else from P4. We use 10 simulation settings and 25 datasets are generated for each setting, as shown in table 1.

Results

Variations in n and p. Each value in the plots shown below is the average obtained over 25 runs. To avoid visual clutter, we do not plot standard error bars. We find that the variance of the algorithms across the runs were low and not significant.

Figure 1 shows the AMI of all four algorithms over 10 different simulation settings for numerical (left) and mixed

(right) data. As expected, GMM performs well on numerical data but not on mixed data while CoupledMC does reasonably well with mixed data. The performance of SCENIC is comparable to that of CoupledMC and GMM. In 2 out of the 10 settings for mixed data, it outperforms all other algorithms. Our algorithm EGMCM consistently outperforms all three algorithms for numerical data. For mixed data, it is the best in 6 out of 10 settings. Similar performance of all four algorithms is observed with respect to ARI in figure 1. Algorithm EGMCM has the highest ARI among the four methods in the case of numerical data and in 8 out of 10 settings for mixed data. In numerical data, the improvement is higher for simulations III-X (large n/p), lower for simulations I, II, indicating that the method seems to work better with large n/p (in accordance with our theoretical analysis). In mixed data, the improvement over GMM is high in all cases but with more improvement in settings with large n/p.

Variations in G. We study the performance of the algorithms when the number of clusters, G, is varied for a single simulation setting: n = 500↪ p = 20. For mixed data, we use 20 numerical and 20 discrete features. Each value in the plots is the average obtained over 25 runs. Figure 2 show the AMI of all four algorithms with increasing number of clusters. The performance of all the four methods drop with increase in G. But EGMCM is consistently better than the rest for all values of G. In figure 2 we also see that the performance of EGMCM is better than the other algorithms, however the increase in ARI is not much except for G = 2 in the numerical case. In general, for numerical data, all four algorithms do not perform well with increasing number of clusters. On mixed data, EGMCM remains better than the other three algorithms.

5.2 Experiments on UCI Benchmark Datasets

We compare the performance of our algorithm on 10 benchmark datasets obtained from the UCI repository [Bache and Lichman, 2013]. Results, in table 3, show that EGMCM performs significantly better (p-value 0.001) than the other algorithms (at 1% significance level; using Friedman-post-hocNemenyi test [Demˇsar, 2006]).

6 Discussion

We present a new clustering algorithm, EGMCM, that, using Gaussian mixture copulas (GMCM), can model a wide range of dependencies beyond those captured by meta-Gaussian distributions. We give the first semiparametric method to estimate the parameters of GMCM that can fit continuous, ordinal and binary data. We theoretically analyze the conditions required to obtain consistent and unbiased parameter estimates that, to our knowledge, is the first such analysis of an ERL approach. We empirically demonstrate that EGMCM outperforms state-of-the-art correlation clustering methods for mixed data. A limitation of EGMCM is that it is uninterpretable with nominal data.

EGMCM is more scalable than other dependency clustering methods for mixed data such as SCENIC and coupledMC. Processing of nominal data contributes to the computational expense of other correlation clustering methods but their runtime remains higher even when restricted to only ordinal and continuous data due to the expense of finding dependencies.

The time complexity of EGMCMmay be further improved using better sampling techniques such as those in [Kalaitzis and Silva, 2013]) and constrained covariance structures like Parsimonious Gaussian Mixture Model families [McNicholas and Murphy, 2008].

Our experiments have focussed on classification accuracy when the number of clusters is known. If unknown, the number of clusters can be inferred, for example, using the Bayesian Information Criterion that was found to be effective in the parameter estimation method for GMCM in [Bhattacharya and Rajan, 2014] for continuous data.

In our theoretical analysis we derive the conditions that ensure that the deviation in the mean, λ(m) gj in the mth iteration, from the initial unbiased EM estimate, tends to zero for large m. The value of λ(m) gj for small m depends on the choice of initial values and variance in the original dataset. If suitable conditions for λ(m) gj to be close to zero can be found for small m, then a low deviation from the original unbiased estimate of the mean can be maintained, which in turn ensures consistent and asymptotically unbiased final estimates. This follows from the fact that, if the maximum likelihood estimate of the mean lies in a small neighbourhood of the initial consistent and asymptotically unbiased estimate, then under certain regularity conditions, the estimate itself is also consistent and asymptotically unbiased [Huber, 1967]. Finding such conditions for small m remains an open problem.

The ERL function is equivalent to the distribution of the ranks for continuous data. For discrete data, the distribution of the ranks depends on the univariate marginals. By ignoring the marginal parameters, ERL loses information in the case of discrete marginals. It would be useful to have a precise characterization of this loss compared to the full likelihood involving marginal parameters. Preliminary work in this direction is in [Gu and Ghosal, 2009; Murray et al., 2013; Hoff et al., 2014], only for the Gaussian copula.