nips nips2013 nips2013-327 knowledge-graph by maker-knowledge-mining

327 nips-2013-The Randomized Dependence Coefficient

Source: pdf

Author: David Lopez-Paz, Philipp Hennig, Bernhard Schölkopf

Abstract: We introduce the Randomized Dependence Coefﬁcient (RDC), a measure of nonlinear dependence between random variables of arbitrary dimension based on the Hirschfeld-Gebelein-R´ nyi Maximum Correlation Coefﬁcient. RDC is deﬁned in e terms of correlation of random non-linear copula projections; it is invariant with respect to marginal distribution transformations, has low computational cost and is easy to implement: just ﬁve lines of R code, included at the end of the paper. 1

Reference: text

Summary: the most important sentenses genereted by tfidf model

sentIndex sentText sentNum sentScore

1 de Abstract We introduce the Randomized Dependence Coefﬁcient (RDC), a measure of nonlinear dependence between random variables of arbitrary dimension based on the Hirschfeld-Gebelein-R´ nyi Maximum Correlation Coefﬁcient. [sent-3, score-0.31]

2 RDC is deﬁned in e terms of correlation of random non-linear copula projections; it is invariant with respect to marginal distribution transformations, has low computational cost and is easy to implement: just ﬁve lines of R code, included at the end of the paper. [sent-4, score-0.342]

3 1 Introduction Measuring statistical dependence between random variables is a fundamental problem in statistics. [sent-5, score-0.153]

4 Commonly used measures of dependence, Pearson’s rho, Spearman’s rank or Kendall’s tau are computationally efﬁcient and theoretically well understood, but consider only a limited class of association patterns, like linear or monotonically increasing functions. [sent-6, score-0.119]

5 The development of non-linear dependence measures is challenging because of the radically larger amount of possible association patterns. [sent-7, score-0.226]

6 Despite these difﬁculties, many non-linear statistical dependence measures have been developed recently. [sent-8, score-0.183]

7 However, these methods exhibit high computational demands (at least quadratic costs in the number of samples for KCCA, HSIC, CHSIC, dCor or MIC), are limited to measuring dependencies between scalar random variables (ACE, MIC) or can be difﬁcult to implement (ACE, MIC). [sent-10, score-0.09]

8 This paper develops the Randomized Dependence Coefﬁcient (RDC), an estimator of the HirschfeldGebelein-R´ nyi Maximum Correlation Coefﬁcient (HGR) addressing the issues listed above. [sent-11, score-0.168]

9 RDC e deﬁnes dependence between two random variables as the largest canonical correlation between random non-linear projections of their respective empirical copula-transformations. [sent-12, score-0.446]

10 RDC is invariant to monotonically increasing transformations, operates on random variables of arbitrary dimension, and has computational cost of O(n log n) with respect to the sample size. [sent-13, score-0.057]

11 The following Section reviews the classic work of Alfr´ d R´ nyi [17], who proposed seven desirable e e fundamental properties of dependence measures, proved to be satisﬁed by the Hirschfeld-GebeleinR´ nyi’s Maximum Correlation Coefﬁcient (HGR). [sent-15, score-0.294]

12 Section 3 introduces the Randomized Depene dence Coefﬁcient as an estimator designed in the spirit of HGR, since HGR itself is computationally intractable. [sent-16, score-0.049]

13 Properties of RDC and its relationship to other non-linear dependence measures are analysed in Section 4. [sent-17, score-0.183]

14 1 2 Hirschfeld-Gebelein-R´ nyi’s Maximum Correlation Coefﬁcient e In 1959 [17], Alfr´ d R´ nyi argued that a measure of dependence ρ∗ : X × Y → [0, 1] between e e random variables X ∈ X and Y ∈ Y should satisfy seven fundamental properties: 1. [sent-19, score-0.312]

15 If (X, Y ) ∼ N (µ, Σ), then ρ∗ (X, Y ) = |ρ(X, Y )|, where ρ is the correlation coefﬁcient. [sent-30, score-0.097]

16 R´ nyi also showed the Hirschfeld-Gebelein-R´ nyi Maximum Correlation Coefﬁcient (HGR) [3, 17] e e to satisfy all these properties. [sent-31, score-0.282]

17 HGR was deﬁned by Gebelein in 1941 [3] as the supremum of Pearson’s correlation coefﬁcient ρ over all Borel-measurable functions f, g of ﬁnite variance: hgr(X, Y ) = sup ρ(f (X), g(Y )), (1) f,g Since the supremum in (1) is over an inﬁnite-dimensional space, HGR is not computable. [sent-32, score-0.182]

18 It is an abstract concept, not a practical dependence measure. [sent-33, score-0.135]

19 In the following we propose a scalable estimator with the same structure as HGR: the Randomized Dependence Coefﬁcient. [sent-34, score-0.027]

20 4 deﬁnes this concept formally, we describe the three necessary steps to construct the RDC statistic: copula-transformation of each of the two random samples (Section 3. [sent-37, score-0.041]

21 2) and computation of the largest canonical correlation between the two sets of non-linear random projections (Section 3. [sent-39, score-0.273]

22 αT Φ(P (x)) P (y) Figure 1: RDC computation for a simple set of samples {(xi , yi )}100 drawn from a noisy circular i=1 pattern: The samples are used to estimate the copula, then mapped with randomly drawn non-linear functions. [sent-42, score-0.046]

23 The RDC is the largest canonical correlation between these non-linear projections. [sent-43, score-0.191]

24 1 Estimation of Copula-Transformations To achieve invariance with respect to transformations on marginal distributions (such as shifts or rescalings), we operate on the empirical copula transformation of the data [14, 15]. [sent-45, score-0.307]

25 , Xd ) with continuous marginal cumulative distribution functions (cdfs) Pi , 1 ≤ i ≤ d. [sent-49, score-0.032]

26 , Pd (Xd )), known as the copula transformation, has uniform marginals: 2 Theorem 1. [sent-56, score-0.175]

27 (Probability Integral Transform [14]) For a random variable X with cdf P , the random variable U := P (X) is uniformly distributed on [0, 1]. [sent-57, score-0.063]

28 Crucially, U preserves the dependence structure of the original random vector X, but ignores each of its d marginal forms [14]. [sent-65, score-0.185]

29 The joint distribution of U is known as the copula of X: Theorem 2. [sent-66, score-0.175]

30 , Xd ) have continuous marginal cdfs Pi , 1 ≤ i ≤ d. [sent-70, score-0.069]

31 , Pd (Xd )), where the distribution C is known as the copula of X. [sent-77, score-0.175]

32 A practical estimator of the univariate cdfs P1 , . [sent-78, score-0.064]

33 , Pd is the empirical cdf : Pn (x) := 1 n n i=1 (3) I(Xi ≤ x), which gives rise to the empirical copula transformations of a multivariate sample: (4) Pn (x) = [Pn,1 (x1 ), . [sent-81, score-0.317]

34 The Massart-Dvoretzky-Kiefer-Wolfowitz inequality [13] can be used to show that empirical copula transformations converge fast to the true transformation as the sample size increases: Theorem 3. [sent-85, score-0.294]

35 sample from a probability distribution over Rd with marginal cdf’s P1 , . [sent-92, score-0.051]

36 Let P (X) be the copula transformation and Pn (X) the empirical copula transformation. [sent-96, score-0.392]

37 Computing Pn (X) involves sorting the marginals of X ∈ Rd×n , thus O(dn log(n)) operations. [sent-98, score-0.038]

38 2 Generation of Random Non-Linear Projections The second step of the RDC computation is to augment the empirical copula transformations with non-linear projections, so that linear methods can subsequently be used to capture non-linear dependencies on the original data. [sent-100, score-0.253]

39 In an elegant result, Rahimi and Recht [16] proved that linear regression on random, non-linear projections of the original feature space can generate high-performance regressors: Theorem 4. [sent-102, score-0.107]

40 , αk for which fk (x) = i=1 αi φ(x; wi ) minimizes the empirical risk c(fk (x), y) has a distance from the c-optimal estimator in F bounded by EP [c(fk (x), y)] − min EP [c(f (x), y)] ≤ O f ∈F 1 1 √ +√ n k LC log 1 δ (6) with probability at least 1 − 2δ. [sent-113, score-0.105]

41 Intuitively, Theorem 4 states that randomly selecting wi in them causes only bounded error. [sent-114, score-0.03]

42 k i=1 αi φ(x; wi ) instead of optimising The choice of the non-linearities φ : R → R is the main and unavoidable assumption in RDC. [sent-115, score-0.03]

43 This choice is a well-known problem common to all non-linear regression methods and has been studied extensively in the theory of regression as the selection of reproducing kernel Hilbert space [19, §3. [sent-116, score-0.108]

44 3 We use random features instead of the Nystr¨ m method because of their smaller memory and como putation requirements [11]. [sent-119, score-0.048]

45 In our experiments, we will use sinusoidal projections, φ(wT x + b) := sin(wT x + b). [sent-120, score-0.044]

46 Arguments favouring this choice are that shift-invariant kernels are approximated with these features when using the appropriate random parameter sampling distribution [16], [4, p. [sent-121, score-0.048]

47 24], and that functions with absolutely integrable Fourier transforms are approxi√ mated with L2 error below O(1/ k) by k of these features [10]. [sent-123, score-0.046]

48 Let the random parameters wi ∼ N (0, sI), bi ∼ N (0, s). [sent-124, score-0.048]

49 Choosing wi to be Normal is analogous to the use of the Gaussian kernel for HSIC, CHSIC or KCCA [16]. [sent-125, score-0.08]

50 Tuning s is analogous to selecting the kernel width, that is, to regularize the non-linearity of the random projections. [sent-126, score-0.068]

51 T T φ(w1 xn + b1 ) · · · φ(wk xn + bk ) (7) the k−th order random non-linear projection from X ∈ Rd×n to Φ(X; k, s) ∈ Rk×n . [sent-139, score-0.034]

52 However, recent techniques using fast Walsh-Hadamard transforms [11] allow computing these feature expansions within a computational cost of O(k log(d)n) and O(k) storage. [sent-141, score-0.06]

53 3 Computation of Canonical Correlations The ﬁnal step of RDC is to compute the linear combinations of the augmented empirical copula transformations that have maximal correlation. [sent-143, score-0.253]

54 Canonical Correlation Analysis (CCA, [7]) is the calculation of pairs of basis vectors (α, β) such that the projections αT X and β T Y of two random samples X ∈ Rp×n and Y ∈ Rq×n are maximally correlated. [sent-144, score-0.105]

55 The correlations between the projected (or canonical) random samples are referred to as canonical correlations. [sent-145, score-0.143]

56 Canonical correlations ρ2 are the solutions to the eigenproblem: 0 −1 Cyy Cyx −1 Cxx Cxy 0 α β = ρ2 α β , (8) where Cxy = cov(X, Y ) and the matrices Cxx and Cyy are assumed to be invertible. [sent-147, score-0.028]

57 Therefore, the largest canonical correlation ρ1 between X and Y is the supremum of the correlation coefﬁcients over their linear projections, that is: ρ1 (X, Y ) = supα,β ρ(αT X, β T Y ). [sent-148, score-0.319]

58 4 Formal Deﬁnition or RDC Given the random samples X ∈ Rp×n and Y ∈ Rq×n and the parameters k ∈ N+ and s ∈ R+ , the Randomized Dependence Coefﬁcient between X and Y is deﬁned as: rdc(X, Y ; k, s) := sup ρ αT Φ(P (X); k, s), β T Φ(P (Y ); k, s) . [sent-151, score-0.064]

59 So there exist continuous (thus Borel measurable) functions f (X) and g(Y ) mapping both X and Y to the sorting ranks of Y , i. [sent-158, score-0.035]

60 Meaningful measures of dependence from ﬁnite samples thus must rely on some form of regularization. [sent-162, score-0.206]

61 Then, with probability greater than 1 − 2δ: hgr(X, Y ; F) − rdc(X, Y ; k) = O m F LC √ +√ n k log 1 δ , (10) where m := ααT + ββ T and n, k denote the sample √ and number of random features. [sent-164, score-0.037]

62 The error O(1/ n) is due to the convergence of CCA’s largest eigenvalue in the ﬁnite sample size regime. [sent-166, score-0.039]

63 This result [8, Theorem 6] is originally obtained by posing CCA as a least squares regression on the product space induced by the feature map √ ψ(x, y) = [φ(x)φ(x)T , φ(y)φ(y)T ,√ 2φ(x)φ(y)T ]T . [sent-167, score-0.043]

64 Because of approximating ψ with k random features, an additional error O(1/ k) is introduced in the least squares regression [16, Lemma 3]. [sent-168, score-0.039]

65 Therefore, an equivalence between RDC and KCCA is established if RDC uses an inﬁnite number of sinusoidal features, the random sampling distribution is set to the inverse Fourier transform of the shift-invariant kernel used by KCCA and the copula-transformations are discarded. [sent-169, score-0.112]

66 As parameters, we here count the kernel function for kernel methods, the basis function and number of random features for RDC, the stopping tolerance for ACE and the search-grid size for MIC, respectively. [sent-172, score-0.165]

67 When using random features φ linear for some neighbourhood around zero (like sinusoids or sigmoids), RDC converges to Spearman’s rank correlation coefﬁcient as s → 0, for any k. [sent-174, score-0.173]

68 2 n2 n2 n2 n log n Testing for independence with RDC: Consider the hypothesis “the two sets of non-linear projections are mutually uncorrelated”. [sent-181, score-0.064]

69 , ρk are the i k 2 5 canonical correlations between Φ(P (X); k, s) and Φ(P (Y ); k, s). [sent-185, score-0.102]

70 5 Experimental Results We performed experiments on both synthetic and real-world data to validate the empirical performance of RDC versus the non-linear dependence measures listed in Table 1. [sent-187, score-0.203]

71 Parameter selection: For RDC, the number of random features is set to k = 20 for both random samples, since no signiﬁcant improvements were observed for larger values. [sent-189, score-0.066]

72 The random feature sampling parameter s is more crucial, and set as follows: when the marginals of u are standard uniforms, w ∼ N (0, sI) and b ∼ N (0, s), then V[wT u + b] = s 1 + d ; therefore, we opt to set 3 1 s to a linear scaling of the input variable dimensionality. [sent-190, score-0.06]

73 1 Synthetic Data Resistance to additive noise: We deﬁne the power of a dependence measure as its ability to discern between dependent and independent samples that share equal marginal forms. [sent-199, score-0.207]

74 In the spirit of Simon and Tibshirani1 , we conducted experiments to estimate the power of RDC as a measure of non-linear dependence. [sent-200, score-0.039]

75 We chose 8 bivariate association patterns, depicted inside little boxes in Figure 3. [sent-201, score-0.09]

76 For each of the 8 association patterns, 500 repetitions of 500 samples were generated, in which the input sample was uniformly distributed on the unit interval. [sent-202, score-0.085]

77 Next, we regenerated the input sample randomly, to generate independent versions of each sample with equal marginals. [sent-203, score-0.038]

78 Figure 3 shows the power for the discussed non-linear dependence measures as the variance of some zero-mean Gaussian additive noise increases from 1/30 to 3. [sent-204, score-0.2]

79 RDC shows worse performance in the linear association pattern due to overﬁtting and in the step-function due to the smoothness prior induced by the sinusoidal features, but has good performance in non-functional association patterns. [sent-205, score-0.13]

80 Running times: Table 2 shows running times for the considered non-linear dependence measures on scalar, uniformly distributed, independent samples of sizes {103 , . [sent-206, score-0.206]

81 Table 2: Average running times (in seconds) for dependence measures on versus sample sizes. [sent-213, score-0.202]

82 0983 — — — Value of statistic in [0, 1]: Figure 4 shows RDC, ACE, dCor, MIC, Pearson’s ρ, Spearman’s rank and Kendall’s τ dependence estimates for 14 different associations of two scalar random samples. [sent-236, score-0.248]

83 When the associations are Gaussian (ﬁrst row), RDC scores values close to the Pearson’s correlation coefﬁcient (Section 2, 7th property). [sent-238, score-0.142]

84 2 Feature Selection in Real-World Data We performed greedy feature selection via dependence maximization [21] on eight real-world datasets. [sent-243, score-0.173]

85 More speciﬁcally, we attempted to construct the subset of features G ⊂ X that minimizes the normalized mean squared regression error (NMSE) of a Gaussian process regressor. [sent-244, score-0.051]

86 We do so by selecting the feature x(i) maximizing dependence between the feature set Gi = {Gi−1 , x(i) } and the target variable y at each iteration i ∈ {1, . [sent-245, score-0.179]

87 Given their quadratic computational demands, dCor, HSIC and CHSIC use up to 1, 000 points when measuring dependence; this constraint only applied on the sarcos and abalone datasets. [sent-252, score-0.086]

88 48 abalone 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 2 4 6 8 10 Number of features Figure 2: Feature selection experiments on real-world datasets. [sent-288, score-0.076]

89 Figure 2 summarizes the results for all datasets and algorithms as the number of selected features increases. [sent-289, score-0.03]

90 6 Conclusion We have presented the randomized dependence coefﬁcient, a lightweight non-linear measure of dependence between multivariate random samples. [sent-291, score-0.349]

91 Constructed as a ﬁnite-dimensional estimator in the spirit of the Hirschfeld-Gebelein-R´ nyi maximum correlation coefﬁcient, RDC performs well e empirically, is scalable to very large datasets, and is easy to adapt to concrete problems. [sent-292, score-0.287]

92 4 xvals cor dCor MIC ACE HSIC CHSIC RDC xvals power. [sent-314, score-0.295]

93 4 xvals 0 20 40 60 80 100 0 xvals 20 40 60 Noise Level 80 100 xvals Figure 3: Power of discussed measures on several bivariate association patterns as noise increases. [sent-321, score-0.543]

94 Insets show the noise-free form of each association pattern. [sent-322, score-0.043]

95 0 Figure 4: RDC, ACE, dCor, MIC, Pearson’s ρ, Spearman’s rank and Kendall’s τ estimates (numbers in tables above plots, in that order) for several bivariate association patterns. [sent-421, score-0.102]

96 A R Source Code rdc <- function(x,y,k=20,s=1/6,f=sin) { x <- cbind(apply(as. [sent-422, score-0.737]

97 Das statistische Problem der Korrelation als Variations- und Eigenwertproblem und sein Zusammenhang mit der Ausgleichsrechnung. [sent-439, score-0.068]

98 Zeitschrift f¨ r Angewandte Mathematik u und Mechanik, 21(6):364–379, 1941. [sent-440, score-0.034]

99 Convergence analysis of kernel canonical correlation analysis: theory and practice. [sent-475, score-0.221]

100 Measuring and testing dependence by correlae tion of distances. [sent-585, score-0.135]

similar papers computed by tfidf model

tfidf for this paper:

wordName wordTfidf (topN-words)

[('rdc', 0.737), ('hgr', 0.207), ('mic', 0.201), ('dcor', 0.184), ('copula', 0.175), ('chsic', 0.168), ('ace', 0.15), ('nyi', 0.141), ('dependence', 0.135), ('typ', 0.134), ('xvals', 0.134), ('hsic', 0.134), ('kcca', 0.122), ('correlation', 0.097), ('coef', 0.09), ('pearson', 0.079), ('canonical', 0.074), ('cca', 0.07), ('projections', 0.064), ('spearman', 0.061), ('xd', 0.059), ('transformations', 0.058), ('cbind', 0.05), ('cxx', 0.05), ('cxy', 0.05), ('cyy', 0.05), ('kernel', 0.05), ('measures', 0.048), ('associations', 0.045), ('kendall', 0.045), ('sinusoidal', 0.044), ('randomized', 0.044), ('association', 0.043), ('pd', 0.038), ('cdfs', 0.037), ('pn', 0.036), ('und', 0.034), ('alfr', 0.034), ('ncol', 0.034), ('rnorm', 0.034), ('marginal', 0.032), ('bivariate', 0.031), ('supremum', 0.031), ('features', 0.03), ('wi', 0.03), ('abalone', 0.03), ('reshef', 0.03), ('rq', 0.029), ('gretton', 0.029), ('measuring', 0.029), ('rank', 0.028), ('fk', 0.028), ('correlations', 0.028), ('sarcos', 0.027), ('cor', 0.027), ('cdf', 0.027), ('estimator', 0.027), ('kely', 0.026), ('nmse', 0.026), ('rahimi', 0.024), ('sch', 0.024), ('samples', 0.023), ('borel', 0.023), ('brownian', 0.023), ('sup', 0.023), ('lc', 0.022), ('expansions', 0.022), ('lkopf', 0.022), ('statistic', 0.022), ('spirit', 0.022), ('ud', 0.022), ('feature', 0.022), ('transformation', 0.022), ('gi', 0.021), ('wk', 0.021), ('regression', 0.021), ('hilbert', 0.021), ('demands', 0.02), ('marginals', 0.02), ('invariant', 0.02), ('largest', 0.02), ('empirical', 0.02), ('sample', 0.019), ('patterns', 0.019), ('sz', 0.018), ('sorting', 0.018), ('seven', 0.018), ('random', 0.018), ('wt', 0.018), ('tolerance', 0.017), ('ranks', 0.017), ('rp', 0.017), ('power', 0.017), ('multivariate', 0.017), ('transforms', 0.016), ('nonlinear', 0.016), ('jmlr', 0.016), ('bk', 0.016), ('selection', 0.016), ('depicted', 0.016)]

similar papers list:

simIndex simValue paperId paperTitle

same-paper 1 0.99999994 327 nips-2013-The Randomized Dependence Coefficient

Author: David Lopez-Paz, Philipp Hennig, Bernhard Schölkopf

2 0.15047716 126 nips-2013-Gaussian Process Conditional Copulas with Applications to Financial Time Series

Author: José Miguel Hernández-Lobato, James R. Lloyd, Daniel Hernández-Lobato

Abstract: The estimation of dependencies between multiple variables is a central problem in the analysis of ﬁnancial time series. A common approach is to express these dependencies in terms of a copula function. Typically the copula function is assumed to be constant but this may be inaccurate when there are covariates that could have a large inﬂuence on the dependence structure of the data. To account for this, a Bayesian framework for the estimation of conditional copulas is proposed. In this framework the parameters of a copula are non-linearly related to some arbitrary conditioning variables. We evaluate the ability of our method to predict time-varying dependencies on several equities and currencies and observe consistent performance gains compared to static copula models and other timevarying copula methods. 1

3 0.10850574 123 nips-2013-Flexible sampling of discrete data correlations without the marginal distributions

Author: Alfredo Kalaitzis, Ricardo Silva

Abstract: Learning the joint dependence of discrete variables is a fundamental problem in machine learning, with many applications including prediction, clustering and dimensionality reduction. More recently, the framework of copula modeling has gained popularity due to its modular parameterization of joint distributions. Among other properties, copulas provide a recipe for combining ﬂexible models for univariate marginal distributions with parametric families suitable for potentially high dimensional dependence structures. More radically, the extended rank likelihood approach of Hoff (2007) bypasses learning marginal models completely when such information is ancillary to the learning task at hand as in, e.g., standard dimensionality reduction problems or copula parameter estimation. The main idea is to represent data by their observable rank statistics, ignoring any other information from the marginals. Inference is typically done in a Bayesian framework with Gaussian copulas, and it is complicated by the fact this implies sampling within a space where the number of constraints increases quadratically with the number of data points. The result is slow mixing when using off-the-shelf Gibbs sampling. We present an efﬁcient algorithm based on recent advances on constrained Hamiltonian Markov chain Monte Carlo that is simple to implement and does not require paying for a quadratic cost in sample size. 1 Contribution There are many ways of constructing multivariate discrete distributions: from full contingency tables in the small dimensional case [1], to structured models given by sparsity constraints [11] and (hierarchies of) latent variable models [6]. More recently, the idea of copula modeling [16] has been combined with such standard building blocks. Our contribution is a novel algorithm for efﬁcient Markov chain Monte Carlo (MCMC) for the copula framework introduced by [7], extending algorithmic ideas introduced by [17]. A copula is a continuous cumulative distribution function (CDF) with uniformly distributed univariate marginals in the unit interval [0, 1]. It complements graphical models and other formalisms that provide a modular parameterization of joint distributions. The core idea is simple and given by the following observation: suppose we are given a (say) bivariate CDF F (y1 , y2 ) with marginals −1 −1 F1 (y1 ) and F2 (y2 ). This CDF can then be rewritten as F (F1 (F1 (y1 )), F2 (F2 (y2 ))). The func−1 −1 tion C(·, ·) given by F (F1 (·), F2 (·)) is a copula. For discrete distributions, this decomposition is not unique but still well-deﬁned [16]. Copulas have found numerous applications in statistics and machine learning since they provide a way of constructing ﬂexible multivariate distributions by mix-and-matching different copulas with different univariate marginals. For instance, one can combine ﬂexible univariate marginals Fi (·) with useful but more constrained high-dimensional copulas. We will not further motivate the use of copula models, which has been discussed at length in recent 1 machine learning publications and conference workshops, and for which comprehensive textbooks exist [e.g., 9]. For a recent discussion on the applications of copulas from a machine learning perspective, [4] provides an overview. [10] is an early reference in machine learning. The core idea dates back at least to the 1950s [16]. In the discrete case, copulas can be difﬁcult to apply: transforming a copula CDF into a probability mass function (PMF) is computationally intractable in general. For the continuous case, a common ˆ trick goes as follows: transform variables by deﬁning ai ≡ Fi (yi ) for an estimate of Fi (·) and then ﬁt a copula density c(·, . . . , ·) to the resulting ai [e.g. 9]. It is not hard to check this breaks down in the discrete case [7]. An alternative is to represent the CDF to PMF transformation for each data point by a continuous integral on a bounded space. Sampling methods can then be used. This trick has allowed many applications of the Gaussian copula to discrete domains. Readers familiar with probit models will recognize the similarities to models where an underlying latent Gaussian ﬁeld is discretized into observable integers as in Gaussian process classiﬁers and ordinal regression [18]. Such models can be indirectly interpreted as special cases of the Gaussian copula. In what follows, we describe in Section 2 the Gaussian copula and the general framework for constructing Bayesian estimators of Gaussian copulas by [7], the extended rank likelihood framework. This framework entails computational issues which are discussed. A recent general approach for MCMC in constrained Gaussian ﬁelds by [17] can in principle be directly applied to this problem as a blackbox, but at a cost that scales quadratically in sample size and as such it is not practical in general. Our key contribution is given in Section 4. An application experiment on the Bayesian Gaussian copula factor model is performed in Section 5. Conclusions are discussed in the ﬁnal section. 2 Gaussian copulas and the extended rank likelihood It is not hard to see that any multivariate Gaussian copula is fully deﬁned by a correlation matrix C, since marginal distributions have no free parameters. In practice, the following equivalent generative model is used to deﬁne a sample U according to a Gaussian copula GC(C): 1. Sample Z from a zero mean Gaussian with covariance matrix C 2. For each Zj , set Uj = Φ(zj ), where Φ(·) is the CDF of the standard Gaussian It is clear that each Uj follows a uniform distribution in [0, 1]. To obtain a model for variables {y1 , y2 , . . . , yp } with marginal distributions Fj (·) and copula GC(C), one can add the deterministic (n) (1) (1) (2) step yj = Fj−1 (uj ). Now, given n samples of observed data Y ≡ {y1 , . . . , yp , y1 , . . . , yp }, one is interested on inferring C via a Bayesian approach and the posterior distribution p(C, θF | Y) ∝ pGC (Y | C, θF )π(C, θF ) where π(·) is a prior distribution, θF are marginal parameters for each Fj (·), which in general might need to be marginalized since they will be unknown, and pGC (·) is the PMF of a (here discrete) distribution with a Gaussian copula and marginals given by θF . Let Z be the underlying latent Gaussians of the corresponding copula for dataset Y. Although Y is a deterministic function of Z, this mapping is not invertible due to the discreteness of the distribution: each marginal Fj (·) has jumps. Instead, the reverse mapping only enforces the constraints where (i ) (i ) (i ) (i ) yj 1 < yj 2 implies zj 1 < zj 2 . Based on this observation, [7] considers the event Z ∈ D(y), where D(y) is the set of values of Z in Rn×p obeying those constraints, that is (k) (k) D(y) ≡ Z ∈ Rn×p : max zj s.t. yj (i) < yj (i) (k) (i) (k) < zj < min zj s.t. yj < yj . Since {Y = y} ⇒ Z(y) ∈ D(y), we have pGC (Y | C, θF ) = pGC (Z ∈ D(y), Y | C, θF ) = pN (Z ∈ D(y) | C) × pGC (Y| Z ∈ D(y), C, θF ), (1) the ﬁrst factor of the last line being that of a zero-mean a Gaussian density function marginalized over D(y). 2 The extended rank likelihood is deﬁned by the ﬁrst factor of (1). With this likelihood, inference for C is given simply by marginalizing p(C, Z | Y) ∝ I(Z ∈ D(y)) pN (Z| C) π(C), (2) the ﬁrst factor of the right-hand side being the usual binary indicator function. Strictly speaking, this is not a fully Bayesian method since partial information on the marginals is ignored. Nevertheless, it is possible to show that under some mild conditions there is information in the extended rank likelihood to consistently estimate C [13]. It has two important properties: ﬁrst, in many applications where marginal distributions are nuisance parameters, this sidesteps any major assumptions about the shape of {Fi (·)} – applications include learning the degree of dependence among variables (e.g., to understand relationships between social indicators as in [7] and [13]) and copula-based dimensionality reduction (a generalization of correlation-based principal component analysis, e.g., [5]); second, MCMC inference in the extended rank likelihood is conceptually simpler than with the joint likelihood, since dropping marginal models will remove complicated entanglements between C and θF . Therefore, even if θF is necessary (when, for instance, predicting missing values of Y), an estimate of C can be computed separately and will not depend on the choice of estimator for {Fi (·)}. The standard model with a full correlation matrix C can be further reﬁned to take into account structure implied by sparse inverse correlation matrices [2] or low rank decompositions via higher-order latent variable models [13], among others. We explore the latter case in section 5. An off-the-shelf algorithm for sampling from (2) is full Gibbs sampling: ﬁrst, given Z, the (full or structured) correlation matrix C can be sampled by standard methods. More to the point, sampling (i) Z is straightforward if for each variable j and data point i we sample Zj conditioned on all other variables. The corresponding distribution is an univariate truncated Gaussian. This is the approach used originally by Hoff. However, mixing can be severely compromised by the sampling of Z, and that is where novel sampling methods can facilitate inference. 3 Exact HMC for truncated Gaussian distributions (i) Hoff’s algorithm modiﬁes the positions of all Zj associated with a particular discrete value of Yj , conditioned on the remaining points. As the number of data points increases, the spread of the hard (i) boundaries on Zj , given by data points of Zj associated with other levels of Yj , increases. This (i) reduces the space in which variables Zj can move at a time. To improve the mixing, we aim to sample from the joint Gaussian distribution of all latent variables (i) Zj , i = 1 . . . n , conditioned on other columns of the data, such that the constraints between them are satisﬁed and thus the ordering in the observation level is conserved. Standard Gibbs approaches for sampling from truncated Gaussians reduce the problem to sampling from univariate truncated Gaussians. Even though each step is computationally simple, mixing can be slow when strong correlations are induced by very tight truncation bounds. In the following, we brieﬂy describe the methodology recently introduced by [17] that deals with the problem of sampling from log p(x) ∝ − 1 x Mx + r x , where x, r ∈ Rn and M is positive 2 deﬁnite, with linear constraints of the form fj x ≤ gj , where fj ∈ Rn , j = 1 . . . m, is the normal vector to some linear boundary in the sample space. Later in this section we shall describe how this framework can be applied to the Gaussian copula extended rank likelihood model. More importantly, the observed rank statistics impose only linear constraints of the form xi1 ≤ xi2 . We shall describe how this special structure can be exploited to reduce the runtime complexity of the constrained sampler from O(n2 ) (in the number of observations) to O(n) in practice. 3.1 Hamiltonian Monte Carlo for the Gaussian distribution Hamiltonian Monte Carlo (HMC) [15] is a MCMC method that extends the sampling space with auxiliary variables so that (ideally) deterministic moves in the joint space brings the sampler to 3 potentially far places in the original variable space. Deterministic moves cannot in general be done, but this is possible in the Gaussian case. The form of the Hamiltonian for the general d-dimensional Gaussian case with mean µ and precision matrix M is: 1 1 H = x Mx − r x + s M−1 s , (3) 2 2 where M is also known in the present context as the mass matrix, r = Mµ and s is the velocity. Both x and s are Gaussian distributed so this Hamiltonian can be seen (up to a constant) as the negative log of the product of two independent Gaussian random variables. The physical interpretation is that of a sum of potential and kinetic energy terms, where the total energy of the system is conserved. In a system where this Hamiltonian function is constant, we can exactly compute its evolution through the pair of differential equations: ˙ x= sH = M−1 s , ˙ s=− xH = −Mx + r . (4) These are solved exactly by x(t) = µ + a sin(t) + b cos(t) , where a and b can be identiﬁed at initial conditions (t = 0) : ˙ a = x(0) = M−1 s , b = x(0) − µ . (5) Therefore, the exact HMC algorithm can be summarised as follows: • Initialise the allowed travel time T and some initial position x0 . • Repeat for HMC samples k = 1 . . . N 1. Sample sk ∼ N (0, M) 2. Use sk and xk to update a and b and store the new position at the end of the trajectory xk+1 = x(T ) as an HMC sample. It can be easily shown that the Markov chain of sampled positions has the desired equilibrium distribution N µ, M−1 [17]. 3.2 Sampling with linear constraints Sampling from multivariate Gaussians does not require any method as sophisticated as HMC, but the plot thickens when the target distribution is truncated by linear constraints of the form Fx ≤ g . Here, F ∈ Rm×n is a constraint matrix whose every row is the normal vector to a linear boundary in the sample space. This is equivalent to sampling from a Gaussian that is conﬁned in the (not necessarily bounded) convex polyhedron {x : Fx ≤ g}. In general, to remain within the boundaries of each wall, once a new velocity has been sampled one must compute all possible collision times with the walls. The smallest of all collision times signiﬁes the wall that the particle should bounce from at that collision time. Figure 1 illustrates the concept with two simple examples on 2 and 3 dimensions. The collision times can be computed analytically and their equations can be found in the supplementary material. We also point the reader to [17] for a more detailed discussion of this implementation. Once the wall to be hit has been found, then position and velocity at impact time are computed and the velocity is reﬂected about the boundary normal1 . The constrained HMC sampler is summarized follows: • Initialise the allowed travel time T and some initial position x0 . • Repeat for HMC samples k = 1 . . . N 1. Sample sk ∼ N (0, M) 2. Use sk and xk to update a and b . 1 Also equivalent to transforming the velocity with a Householder reﬂection matrix about the bounding hyperplane. 4 1 2 3 4 1 2 3 4 Figure 1: Left: Trajectories of the ﬁrst 40 iterations of the exact HMC sampler on a 2D truncated Gaussian. A reﬂection of the velocity can clearly be seen when the particle meets wall #2 . Here, the constraint matrix F is a 4 × 2 matrix. Center: The same example after 40000 samples. The coloring of each sample indicates its density value. Right: The anatomy of a 3D Gaussian. The walls are now planes and in this case F is a 2 × 3 matrix. Figure best seen in color. 3. Reset remaining travel time Tleft ← T . Until no travel time is left or no walls can be reached (no solutions exist), do: (a) Compute impact times with all walls and pick the smallest one, th (if a solution exists). (b) Compute v(th ) and reﬂect it about the hyperplane fh . This is the updated velocity after impact. The updated position is x(th ) . (c) Tleft ← Tleft − th 4. Store the new position at the end of the trajectory xk+1 as an HMC sample. In general, all walls are candidates for impact, so the runtime of the sampler is linear in m , the number of constraints. This means that the computational load is concentrated in step 3(a). Another consideration is that of the allocated travel time T . Depending on the shape of the bounding polyhedron and the number of walls, a very large travel time can induce many more bounces thus requiring more computations per sample. On the other hand, a very small travel time explores the distribution more locally so the mixing of the chain can suffer. What constitutes a given travel time “large” or “small” is relative to the dimensionality, the number of constraints and the structure of the constraints. Due to the nature of our problem, the number of constraints, when explicitly expressed as linear functions, is O(n2 ) . Clearly, this restricts any direct application of the HMC framework for Gaussian copula estimation to small-sample (n) datasets. More importantly, we show how to exploit the structure of the constraints to reduce the number of candidate walls (prior to each bounce) to O(n) . 4 HMC for the Gaussian Copula extended rank likelihood model Given some discrete data Y ∈ Rn×p , the task is to infer the correlation matrix of the underlying Gaussian copula. Hoff’s sampling algorithm proceeds by alternating between sampling the continu(i) (i) ous latent representation Zj of each Yj , for i = 1 . . . n, j = 1 . . . p , and sampling a covariance matrix from an inverse-Wishart distribution conditioned on the sampled matrix Z ∈ Rn×p , which is then renormalized as a correlation matrix. From here on, we use matrix notation for the samples, as opposed to the random variables – with (i) Zi,j replacing Zj , Z:,j being a column of Z, and Z:,\j being the submatrix of Z without the j-th column. In a similar vein to Hoff’s sampling algorithm, we replace the successive sampling of each Zi,j conditioned on Zi,\j (a conditional univariate truncated Gaussian) with the simultaneous sampling of Z:,j conditioned on Z:,\j . This is done through an HMC step from a conditional multivariate truncated Gaussian. The added beneﬁt of this HMC step over the standard Gibbs approach, is that of a handle for regulating the trade-off between exploration and runtime via the allocated travel time T . Larger travel times potentially allow for larger moves in the sample space, but it comes at a cost as explained in the sequel. 5 4.1 The Hough envelope algorithm The special structure of constraints. Recall that the number of constraints is quadratic in the dimension of the distribution. This is because every Z sample must satisfy the conditions of the event Z ∈ D(y) of the extended rank likelihood (see Section 2). In other words, for any column Z:,j , all entries are organised into a partition L(j) of |L(j) | levels, the number of unique values observed for the discrete or ordinal variable Y (j) . Thereby, for any two adjacent levels lk , lk+1 ∈ L(j) and any pair i1 ∈ lk , i2 ∈ lk+1 , it must be true that Zli ,j < Zli+1 ,j . Equivalently, a constraint f exists where fi1 = 1, fi2 = −1 and g = 0 . It is easy to see that O(n2 ) of such constraints are induced by the order statistics of the j-th variable. To deal with this boundary explosion, we developed the Hough Envelope algorithm to search efﬁciently, within all pairs in {Z:,j }, in practically linear time. Recall in HMC (section 3.2) that the trajectory of the particle, x(t), is decomposed as xi (t) = ai sin(t) + bi cos(t) + µi , (6) and there are n such functions, grouped into a partition of levels as described above. The Hough envelope2 is found for every pair of adjacent levels. We illustrate this with an example of 10 dimensions and two levels in Figure 2, without loss of generalization to any number of levels or dimensions. Assume we represent trajectories for points in level lk with blue curves, and points in lk+1 with red curves. Assuming we start with a valid state, at time t = 0 all red curves are above all blue curves. The goal is to ﬁnd the smallest t where a blue curve meets a red curve. This will be our collision time where a bounce will be necessary. 5 3 1 2 Figure 2: The trajectories xj (t) of each component are sinusoid functions. The right-most green dot signiﬁes the wall and the time th of the earliest bounce, where the ﬁrst inter-level pair (that is, any two components respectively from the blue and red level) becomes equal, in this case the constraint activated being xblue2 = xred2 . 4 4 5 1 2 3 0.2 0.4 0.6 t 0.8 1 1.2 1.4 1. First we ﬁnd the largest component bluemax of the blue level at t = 0. This takes O(n) time. Clearly, this will be the largest component until its sinusoid intersects that of any other component. 2. To ﬁnd the next largest component, compute the roots of xbluemax (t) − xi (t) = 0 for all components and pick the smallest (earliest) one (represented by a green dot). This also takes O(n) time. 3. Repeat this procedure until a red sinusoid crosses the highest running blue sinusoid. When this happens, the time of earliest bounce and its constraint are found. In the worst-case scenario, n such repetitions have to be made, but in practice we can safely assume an ﬁxed upper bound h on the number of blue crossings before a inter-level crossing occurs. In experiments, we found h << n, no more than 10 in simulations with hundreds of thousands of curves. Thus, this search strategy takes O(n) time in practice to complete, mirroring the analysis of other output-sensitive algorithms such as the gift wrapping algorithm for computing convex hulls [8]. Our HMC sampling approach is summarized in Algorithm 1. 2 The name is inspired from the fact that each xi (t) is the sinusoid representation, in angle-distance space, of all lines that pass from the (ai , bi ) point in a − b space. A representation known in image processing as the Hough transform [3]. 6 Algorithm 1 HMC for GCERL # Notation: T MN (µ, C, F) is a truncated multivariate normal with location vector µ, scale matrix C and constraints encoded by F and g = 0 . # IW(df, V0 ) is an inverse-Wishart prior with degrees of freedom df and scale matrix V0 . Input: Y ∈ Rn×p , allocated travel time T , a starting Z and variable covariance V ∈ Rp×p , df = p + 2, V0 = df Ip and chain size N . Generate constraints F(j) from Y:,j , for j = 1 . . . p . for samples k = 1 . . . N do # Resample Z as follows: for variables j = 1 . . . p do −1 −1 2 Compute parameters: σj = Vjj − Vj,\j V\j,\j V\j,j , µj = Z:,\j V\j,\j V\j,j . 2 Get one sample Z:,j ∼ T MN µj , σj I, F(j) efﬁciently by using the Hough Envelope algorithm, see section 4.1. end for Resample V ∼ IW(df + n, V0 + Z Z) . Compute correlation matrix C, s.t. Ci,j = Vi,j / Vi,i Vj,j and store sample, C(k) ← C . end for 5 An application on the Bayesian Gausian copula factor model In this section we describe an experiment that highlights the beneﬁts of our HMC treatment, compared to a state-of-the-art parameter expansion (PX) sampling scheme. During this experiment we ask the important question: “How do the two schemes compare when we exploit the full-advantage of the HMC machinery to jointly sample parameters and the augmented data Z, in a model of latent variables and structured correlations?” We argue that under such circumstances the superior convergence speed and mixing of HMC undeniably compensate for its computational overhead. Experimental setup In this section we provide results from an application on the Gaussian copula latent factor model of [13] (Hoff’s model [7] for low-rank structured correlation matrices). We modify the parameter expansion (PX) algorithm used by [13] by replacing two of its Gibbs steps with a single HMC step. We show a much faster convergence to the true mode with considerable support on its vicinity. We show that unlike the HMC, the PX algorithm falls short of properly exploring the posterior in any reasonable ﬁnite amount of time, even for small models, even for small samples. Worse, PX fails in ways one cannot easily detect. Namely, we sample each row of the factor loadings matrix Λ jointly with the corresponding column of the augmented data matrix Z, conditioning on the higher-order latent factors. This step is analogous to Pakman and Paninski’s [17, sec.3.1] use of HMC in the context of a binary probit model (the extension to many levels in the discrete marginal is straightforward with direct application of the constraint matrix F and the Hough envelope algorithm). The sampling of the higher level latent factors remains identical to [13]. Our scheme involves no parameter expansion. We do however interweave the Gibbs step for the Z matrix similarly to Hoff. This has the added beneﬁt of exploring the Z sample space within their current boundaries, complementing the joint (λ, z) sampling which moves the boundaries jointly. The value of such ”interweaving” schemes has been addressed in [19]. Results We perform simulations of 10000 iterations, n = 1000 observations (rows of Y), travel time π/2 for HMC with the setups listed in the following table, along with the elapsed times of each sampling scheme. These experiments were run on Intel COREi7 desktops with 4 cores and 8GB of RAM. Both methods were parallelized across the observed variables (p). Figure p (vars) k (latent factors) M (ordinal levels) elapsed (mins): HMC PX 3(a) : 20 5 2 115 8 3(b) : 10 3 2 80 6 10 3 5 203 16 3(c) : Many functionals of the loadings matrix Λ can be assessed. We focus on reconstructing the true (low-rank) correlation matrix of the Gaussian copula. In particular, we summarize the algorithm’s 7 outcome with the root mean squared error (RMSE) of the differences between entries of the ground-truth correlation matrix and the implied correlation matrix at each iteration of a MCMC scheme (so the following plots looks like a time-series of 10000 timepoints), see Figures 3(a), 3(b) and 3(c) . (a) (b) (c) Figure 3: Reconstruction (RMSE per iteration) of the low-rank structured correlation matrix of the Gaussian copula and its histogram (along the left side). (a) Simulation setup: 20 variables, 5 factors, 5 levels. HMC (blue) reaches a better mode faster (in iterations/CPU-time) than PX (red). Even more importantly the RMSE posterior samples of PX are concentrated in a much smaller region compared to HMC, even after 10000 iterations. This illustrates that PX poorly explores the true distribution. (b) Simulation setup: 10 vars, 3 factors, 2 levels. We observe behaviors similar to Figure 3(a). Note that the histogram counts RMSEs after the burn-in period of PX (iteration #500). (c) Simulation setup: 10 vars, 3 factors, 5 levels. We observe behaviors similar to Figures 3(a) and 3(b) but with a thinner tail for HMC. Note that the histogram counts RMSEs after the burn-in period of PX (iteration #2000). Main message HMC reaches a better mode faster (iterations/CPUtime). Even more importantly the RMSE posterior samples of PX are concentrated in a much smaller region compared to HMC, even after 10000 iterations. This illustrates that PX poorly explores the true distribution. As an analogous situation we refer to the top and bottom panels of Figure 14 of Radford Neal’s slice sampler paper [14]. If there was no comparison against HMC, there would be no evidence from the PX plot alone that the algorithm is performing poorly. This mirrors Radford Neal’s statement opening Section 8 of his paper: “a wrong answer is obtained without any obvious indication that something is amiss”. The concentration on the posterior mode of PX in these simulations is misleading of the truth. PX might seen a bit simpler to implement, but it seems one cannot avoid using complex algorithms for complex models. We urge practitioners to revisit their past work with this model to ﬁnd out by how much credible intervals of functionals of interest have been overconﬁdent. Whether trivially or severely, our algorithm offers the ﬁrst principled approach for checking this out. 6 Conclusion Sampling large random vectors simultaneously in order to improve mixing is in general a very hard problem, and this is why clever methods such as HMC or elliptical slice sampling [12] are necessary. We expect that the method here developed is useful not only for those with data analysis problems within the large family of Gaussian copula extended rank likelihood models, but the method itself and its behaviour might provide some new insights on MCMC sampling in constrained spaces in general. Another direction of future work consists of exploring methods for elliptical copulas, and related possible extensions of general HMC for non-Gaussian copula models. Acknowledgements The quality of this work has beneﬁted largely from comments by our anonymous reviewers and useful discussions with Simon Byrne and Vassilios Stathopoulos. Research was supported by EPSRC grant EP/J013293/1. 8 References [1] Y. Bishop, S. Fienberg, and P. Holland. Discrete Multivariate Analysis: Theory and Practice. MIT Press, 1975. [2] A. Dobra and A. Lenkoski. Copula Gaussian graphical models and their application to modeling functional disability data. Annals of Applied Statistics, 5:969–993, 2011. [3] R. O. Duda and P. E. Hart. Use of the Hough transformation to detect lines and curves in pictures. Communications of the ACM, 15(1):11–15, 1972. [4] G. Elidan. Copulas and machine learning. Proceedings of the Copulae in Mathematical and Quantitative Finance workshop, to appear, 2013. [5] F. Han and H. Liu. Semiparametric principal component analysis. Advances in Neural Information Processing Systems, 25:171–179, 2012. [6] G. Hinton and R. Salakhutdinov. Reducing the dimensionality of data with neural networks. Science, 313(5786):504–507, 2006. [7] P. Hoff. Extending the rank likelihood for semiparametric copula estimation. Annals of Applied Statistics, 1:265–283, 2007. [8] R. Jarvis. On the identiﬁcation of the convex hull of a ﬁnite set of points in the plane. Information Processing Letters, 2(1):18–21, 1973. [9] H. Joe. Multivariate Models and Dependence Concepts. Chapman-Hall, 1997. [10] S. Kirshner. Learning with tree-averaged densities and distributions. Neural Information Processing Systems, 2007. [11] S. Lauritzen. Graphical Models. Oxford University Press, 1996. [12] I. Murray, R. Adams, and D. MacKay. Elliptical slice sampling. JMLR Workshop and Conference Proceedings: AISTATS 2010, 9:541–548, 2010. [13] J. Murray, D. Dunson, L. Carin, and J. Lucas. Bayesian Gaussian copula factor models for mixed data. Journal of the American Statistical Association, to appear, 2013. [14] R. Neal. Slice sampling. The Annals of Statistics, 31:705–767, 2003. [15] R. Neal. MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, pages 113–162, 2010. [16] R. Nelsen. An Introduction to Copulas. Springer-Verlag, 2007. [17] A. Pakman and L. Paninski. Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians. arXiv:1208.4118, 2012. [18] C. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006. [19] Y. Yu and X. L. Meng. To center or not to center: That is not the question — An ancillaritysufﬁciency interweaving strategy (ASIS) for boosting MCMC efﬁciency. Journal of Computational and Graphical Statistics, 20(3):531–570, 2011. 9

4 0.08393386 173 nips-2013-Least Informative Dimensions

Author: Fabian Sinz, Anna Stockl, January Grewe, January Benda

Abstract: We present a novel non-parametric method for ﬁnding a subspace of stimulus features that contains all information about the response of a system. Our method generalizes similar approaches to this problem such as spike triggered average, spike triggered covariance, or maximally informative dimensions. Instead of maximizing the mutual information between features and responses directly, we use integral probability metrics in kernel Hilbert spaces to minimize the information between uninformative features and the combination of informative features and responses. Since estimators of these metrics access the data via kernels, are easy to compute, and exhibit good theoretical convergence properties, our method can easily be generalized to populations of neurons or spike patterns. By using a particular expansion of the mutual information, we can show that the informative features must contain all information if we can make the uninformative features independent of the rest. 1

5 0.068858713 310 nips-2013-Statistical analysis of coupled time series with Kernel Cross-Spectral Density operators.

Author: Michel Besserve, Nikos K. Logothetis, Bernhard Schölkopf

Abstract: Many applications require the analysis of complex interactions between time series. These interactions can be non-linear and involve vector valued as well as complex data structures such as graphs or strings. Here we provide a general framework for the statistical analysis of these dependencies when random variables are sampled from stationary time-series of arbitrary objects. To achieve this goal, we study the properties of the Kernel Cross-Spectral Density (KCSD) operator induced by positive deﬁnite kernels on arbitrary input domains. This framework enables us to develop an independence test between time series, as well as a similarity measure to compare different types of coupling. The performance of our test is compared to the HSIC test using i.i.d. assumptions, showing improvements in terms of detection errors, as well as the suitability of this approach for testing dependency in complex dynamical systems. This similarity measure enables us to identify different types of interactions in electrophysiological neural time series. 1

6 0.061252389 76 nips-2013-Correlated random features for fast semi-supervised learning

7 0.060511623 44 nips-2013-B-test: A Non-parametric, Low Variance Kernel Two-sample Test

8 0.053547669 116 nips-2013-Fantope Projection and Selection: A near-optimal convex relaxation of sparse PCA

9 0.052345306 9 nips-2013-A Kernel Test for Three-Variable Interactions

10 0.051752258 281 nips-2013-Robust Low Rank Kernel Embeddings of Multivariate Distributions

11 0.045410726 302 nips-2013-Sparse Inverse Covariance Estimation with Calibration

12 0.044004049 68 nips-2013-Confidence Intervals and Hypothesis Testing for High-Dimensional Statistical Models

13 0.043293159 153 nips-2013-Learning Feature Selection Dependencies in Multi-task Learning

14 0.04154858 55 nips-2013-Bellman Error Based Feature Generation using Random Projections on Sparse Spaces

15 0.037743922 252 nips-2013-Predictive PAC Learning and Process Decompositions

16 0.037533179 145 nips-2013-It is all in the noise: Efficient multi-task Gaussian process inference with structured residuals

17 0.036986452 170 nips-2013-Learning with Invariance via Linear Functionals on Reproducing Kernel Hilbert Space

18 0.036926363 144 nips-2013-Inverse Density as an Inverse Problem: the Fredholm Equation Approach

19 0.036053423 282 nips-2013-Robust Multimodal Graph Matching: Sparse Coding Meets Graph Matching

20 0.034007933 211 nips-2013-Non-Linear Domain Adaptation with Boosting

similar papers computed by lsi model

lsi for this paper:

topicId topicWeight

[(0, 0.104), (1, 0.043), (2, 0.026), (3, 0.014), (4, -0.004), (5, 0.045), (6, -0.002), (7, 0.008), (8, -0.028), (9, 0.019), (10, -0.01), (11, -0.015), (12, -0.076), (13, -0.032), (14, 0.081), (15, 0.074), (16, -0.038), (17, 0.053), (18, -0.109), (19, -0.0), (20, -0.077), (21, 0.015), (22, -0.001), (23, 0.05), (24, 0.107), (25, -0.007), (26, 0.102), (27, 0.028), (28, 0.047), (29, -0.127), (30, 0.027), (31, -0.049), (32, -0.043), (33, 0.003), (34, -0.042), (35, 0.006), (36, -0.021), (37, 0.042), (38, 0.023), (39, 0.076), (40, 0.024), (41, -0.063), (42, 0.071), (43, -0.036), (44, 0.032), (45, 0.085), (46, 0.019), (47, -0.042), (48, -0.011), (49, 0.048)]

similar papers list:

simIndex simValue paperId paperTitle

same-paper 1 0.88256615 327 nips-2013-The Randomized Dependence Coefficient

Author: David Lopez-Paz, Philipp Hennig, Bernhard Schölkopf

2 0.7944665 126 nips-2013-Gaussian Process Conditional Copulas with Applications to Financial Time Series

Author: José Miguel Hernández-Lobato, James R. Lloyd, Daniel Hernández-Lobato

3 0.640315 123 nips-2013-Flexible sampling of discrete data correlations without the marginal distributions

Author: Alfredo Kalaitzis, Ricardo Silva

4 0.54970449 9 nips-2013-A Kernel Test for Three-Variable Interactions

Author: Dino Sejdinovic, Arthur Gretton, Wicher Bergsma

Abstract: We introduce kernel nonparametric tests for Lancaster three-variable interaction and for total independence, using embeddings of signed measures into a reproducing kernel Hilbert space. The resulting test statistics are straightforward to compute, and are used in powerful interaction tests, which are consistent against all alternatives for a large family of reproducing kernels. We show the Lancaster test to be sensitive to cases where two independent causes individually have weak inﬂuence on a third dependent variable, but their combined effect has a strong inﬂuence. This makes the Lancaster test especially suited to ﬁnding structure in directed graphical models, where it outperforms competing nonparametric tests in detecting such V-structures.

5 0.52741534 44 nips-2013-B-test: A Non-parametric, Low Variance Kernel Two-sample Test

Author: Wojciech Zaremba, Arthur Gretton, Matthew Blaschko

Abstract: A family of maximum mean discrepancy (MMD) kernel two-sample tests is introduced. Members of the test family are called Block-tests or B-tests, since the test statistic is an average over MMDs computed on subsets of the samples. The choice of block size allows control over the tradeoff between test power and computation time. In this respect, the B-test family combines favorable properties of previously proposed MMD two-sample tests: B-tests are more powerful than a linear time test where blocks are just pairs of samples, yet they are more computationally efﬁcient than a quadratic time test where a single large block incorporating all the samples is used to compute a U-statistic. A further important advantage of the B-tests is their asymptotically Normal null distribution: this is by contrast with the U-statistic, which is degenerate under the null hypothesis, and for which estimates of the null distribution are computationally demanding. Recent results on kernel selection for hypothesis testing transfer seamlessly to the B-tests, yielding a means to optimize test power via kernel choice. 1

6 0.48868257 173 nips-2013-Least Informative Dimensions

7 0.48136449 145 nips-2013-It is all in the noise: Efficient multi-task Gaussian process inference with structured residuals

8 0.47754788 310 nips-2013-Statistical analysis of coupled time series with Kernel Cross-Spectral Density operators.

9 0.46484363 76 nips-2013-Correlated random features for fast semi-supervised learning

10 0.45231321 153 nips-2013-Learning Feature Selection Dependencies in Multi-task Learning

11 0.43565392 43 nips-2013-Auxiliary-variable Exact Hamiltonian Monte Carlo Samplers for Binary Distributions

12 0.43191898 130 nips-2013-Generalizing Analytic Shrinkage for Arbitrary Covariance Structures

13 0.43150085 302 nips-2013-Sparse Inverse Covariance Estimation with Calibration

14 0.41632053 117 nips-2013-Fast Algorithms for Gaussian Noise Invariant Independent Component Analysis

15 0.40291691 281 nips-2013-Robust Low Rank Kernel Embeddings of Multivariate Distributions

16 0.39807072 156 nips-2013-Learning Kernels Using Local Rademacher Complexity

17 0.38643253 306 nips-2013-Speeding up Permutation Testing in Neuroimaging

18 0.38323593 110 nips-2013-Estimating the Unseen: Improved Estimators for Entropy and other Properties

19 0.38155639 144 nips-2013-Inverse Density as an Inverse Problem: the Fredholm Equation Approach

20 0.37934643 225 nips-2013-One-shot learning and big data with n=2

similar papers computed by lda model

lda for this paper:

topicId topicWeight

[(2, 0.063), (16, 0.042), (19, 0.013), (33, 0.109), (34, 0.079), (36, 0.011), (39, 0.011), (41, 0.04), (49, 0.034), (56, 0.074), (70, 0.022), (71, 0.305), (85, 0.023), (89, 0.028), (93, 0.048), (95, 0.016)]

similar papers list:

simIndex simValue paperId paperTitle

same-paper 1 0.7042616 327 nips-2013-The Randomized Dependence Coefficient

Author: David Lopez-Paz, Philipp Hennig, Bernhard Schölkopf

2 0.60046422 226 nips-2013-One-shot learning by inverting a compositional causal process

Author: Brenden M. Lake, Ruslan Salakhutdinov, Josh Tenenbaum

Abstract: People can learn a new visual class from just one example, yet machine learning algorithms typically require hundreds or thousands of examples to tackle the same problems. Here we present a Hierarchical Bayesian model based on compositionality and causality that can learn a wide range of natural (although simple) visual concepts, generalizing in human-like ways from just one image. We evaluated performance on a challenging one-shot classiﬁcation task, where our model achieved a human-level error rate while substantially outperforming two deep learning models. We also tested the model on another conceptual task, generating new examples, by using a “visual Turing test” to show that our model produces human-like performance. 1

3 0.57213414 153 nips-2013-Learning Feature Selection Dependencies in Multi-task Learning

Author: Daniel Hernández-Lobato, José Miguel Hernández-Lobato

Abstract: A probabilistic model based on the horseshoe prior is proposed for learning dependencies in the process of identifying relevant features for prediction. Exact inference is intractable in this model. However, expectation propagation offers an approximate alternative. Because the process of estimating feature selection dependencies may suffer from over-ﬁtting in the model proposed, additional data from a multi-task learning scenario are considered for induction. The same model can be used in this setting with few modiﬁcations. Furthermore, the assumptions made are less restrictive than in other multi-task methods: The different tasks must share feature selection dependencies, but can have different relevant features and model coefﬁcients. Experiments with real and synthetic data show that this model performs better than other multi-task alternatives from the literature. The experiments also show that the model is able to induce suitable feature selection dependencies for the problems considered, only from the training data. 1

4 0.5086695 123 nips-2013-Flexible sampling of discrete data correlations without the marginal distributions

Author: Alfredo Kalaitzis, Ricardo Silva

5 0.50630212 51 nips-2013-Bayesian entropy estimation for binary spike train data using parametric prior knowledge

Author: Evan W. Archer, Il M. Park, Jonathan W. Pillow

Abstract: Shannon’s entropy is a basic quantity in information theory, and a fundamental building block for the analysis of neural codes. Estimating the entropy of a discrete distribution from samples is an important and difﬁcult problem that has received considerable attention in statistics and theoretical neuroscience. However, neural responses have characteristic statistical structure that generic entropy estimators fail to exploit. For example, existing Bayesian entropy estimators make the naive assumption that all spike words are equally likely a priori, which makes for an inefﬁcient allocation of prior probability mass in cases where spikes are sparse. Here we develop Bayesian estimators for the entropy of binary spike trains using priors designed to ﬂexibly exploit the statistical structure of simultaneouslyrecorded spike responses. We deﬁne two prior distributions over spike words using mixtures of Dirichlet distributions centered on simple parametric models. The parametric model captures high-level statistical features of the data, such as the average spike count in a spike word, which allows the posterior over entropy to concentrate more rapidly than with standard estimators (e.g., in cases where the probability of spiking differs strongly from 0.5). Conversely, the Dirichlet distributions assign prior mass to distributions far from the parametric model, ensuring consistent estimates for arbitrary distributions. We devise a compact representation of the data and prior that allow for computationally efﬁcient implementations of Bayesian least squares and empirical Bayes entropy estimators with large numbers of neurons. We apply these estimators to simulated and real neural data and show that they substantially outperform traditional methods.

6 0.5005331 182 nips-2013-Manifold-based Similarity Adaptation for Label Propagation

7 0.49952883 248 nips-2013-Point Based Value Iteration with Optimal Belief Compression for Dec-POMDPs

8 0.49722725 310 nips-2013-Statistical analysis of coupled time series with Kernel Cross-Spectral Density operators.

9 0.49313292 201 nips-2013-Multi-Task Bayesian Optimization

10 0.48860085 173 nips-2013-Least Informative Dimensions

11 0.48710793 171 nips-2013-Learning with Noisy Labels

12 0.48602405 350 nips-2013-Wavelets on Graphs via Deep Learning

13 0.48368731 45 nips-2013-BIG & QUIC: Sparse Inverse Covariance Estimation for a Million Variables

14 0.48328224 341 nips-2013-Universal models for binary spike patterns using centered Dirichlet processes

15 0.48302457 49 nips-2013-Bayesian Inference and Online Experimental Design for Mapping Neural Microcircuits

16 0.48207116 191 nips-2013-Minimax Optimal Algorithms for Unconstrained Linear Optimization

17 0.4820396 318 nips-2013-Structured Learning via Logistic Regression

18 0.48184574 262 nips-2013-Real-Time Inference for a Gamma Process Model of Neural Spiking

19 0.4806689 308 nips-2013-Spike train entropy-rate estimation using hierarchical Dirichlet process priors

20 0.47973722 99 nips-2013-Dropout Training as Adaptive Regularization