Parametric Copula-GP model for analyzing multidimensional neuronal and behavioral relationships
Authors: Nina Kudryashova
Theoklitos Amvrosiadis
Nathalie Dupuy
Nathalie Rochefort
Arno Onken
Venue:  (2021)
Keywords: hierarchical vine copula model
copula
gaussian process
mutual information
neuroscience
neuronal activity
calcium imaging
visual cortex

Abstract: One of the main challenges in current systems neuroscience is the analysis of high-dimensional neuronal and behavioral data that are characterized by different statistics and timescales of the recorded variables. We propose a parametric copula model which separates the statistics of the individual variables from their dependence structure, and escapes the curse of dimensionality by using vine copula constructions. We use a Bayesian framework with Gaussian Process (GP) priors over copula parameters, conditioned on a continuous task-related variable. We improve the flexibility of this method by 1) using non-parametric conditional (rather than unconditional) marginals; 2) linearly mixing copula elements with qualitatively different tail dependencies. We validate the model on synthetic data and compare its performance in estimating mutual information against the commonly used non-parametric algorithms. Our model provides accurate information estimates when the dependencies in the data match the parametric copulas used in our framework. Moreover, even when the exact density estimation with a parametric model is not possible, our Copula-GP model is still able to provide reasonable information estimates, close to the ground truth and comparable to those obtained with a neural network estimator. Finally, we apply our framework to real neuronal and behavioral recordings obtained in awake mice. We demonstrate the ability of our framework to 1) produce accurate and interpretable bivariate models for the analysis of inter-neuronal noise correlations or behavioral modulations; 2) expand to more than 100 dimensions and measure information content in the whole-population statistics. These results demonstrate that the Copula-GP framework is particularly useful for the analysis of complex multidimensional relationships between neuronal, sensory and behavioral data.

Under review as a conference paper at ICLR 2021
PARAMETRIC COPULA-GP MODEL FOR ANALYZING MULTIDIMENSIONAL NEURONAL AND BEHAVIORAL RELATIONSHIPS
Anonymous authors Paper under double-bling review
ABSTRACT
One of the main challenges in current systems neuroscience is the analysis of high-dimensional neuronal and behavioral data that are characterized by differ- ent statistics and timescales of the recorded variables. We propose parametric copula model which separates the statistics of the individual variables from their dependence structure, and escapes the curse of dimensionality by using vine cop- ula constructions. We use a Bayesian framework with Gaussian Process (GP) priors over copula parameters, conditioned on continuous task-related variable. We improve the flexibility of this method by 1) using non-parametric conditional (rather than unconditional) marginals; 2) linearly mixing copula elements with qualitatively different tail dependencies. We validate the model synthetic data and compare its performance in estimating mutual information against the com- monly used non- conaamametric algorithms. Our model provides accurate information estimates when the dependencies the data match the parametric copulas used in our framework. Moreover, when the exact density estimation with parametric model is not possible, our Copula-GP GP model is still able to provide reasonable information estimates, close to the ground truth and comparable to those obtained with a neural network estimator. Finally, we apply our framework to real neuronal and behavioral recordings obtained in awake mice. We demonstrate the ability of our framework to 1) produce accurate and interpretable bivariate models for the analysis of nnter-neuronal noise correlations or behavioral modulations; ex- pand to more than 100 dimensions and measure information content in the whole- population statistics. These results demonstrate that the Copula- G GP framework is particularly useful for the analysis of complex multidimensional relationships between neuronal, sensory and behavioral data.
INTRODUCTION
Recent advances in imaging and recording techniques have enabled monitoring the activity of hun- dreds several thousands of neurons simultaneously (Jun et al. 2017 Helmchen 2009 Dombeck jt al. 2007). These recordings can be made in awake animals engaged in specifically designed tasks or natural behavior (Stringer et al., 2019 Pakan al. [0118 b). which further augments these al- ready large datasets with a variety of behavioral variables. These complex high dimensional datasets necessitate the development of novel analytical approaches (Saxena & Cunningham, 2019 Steven- son & Kording 2011 Staude et al., 2010) to address two central questions of systems and behavioral neuroscience: how do populations of neurons encode information? And how does this neuronal ac- tivity correspond to the observed behavior? In machine learning terms, both of these questions translate into understanding the high-dimensional multivariate dependencies between the recorded variables (Kohn al. 2016 Shimazaki et al. 2012 Ince et al., 2010 Shamir & Sompolinsky 2004).
There are two major methods suitable for recording the activity of large populations of neurons from behaving animals: the multi-electrode probes (Jun et al. 2017). and calcium imaging meth- ods (Grienberger et al. 2015 Helmchen 2009 Dombeck et al. 2007) that use changes in intracel- lular calcium concentration
proxy for neuronal spiking activity at lower temporal precision. While neuronal spiking occurs on a temporal scale milliseconds, the behavior spans the timescales from milliseconds to hours and even days (Mathis al. 2018). As result, the recorded neuronal

---

Under review as a conference paper at ICLR 2021
and behavioral variables may operate different timescales and exhibit different statistics, which further complicates the statistical analysis of these datasets.
The natural approach to modeling statistical dependencies between the variables with drastically different statistics is based on copulas, which separate marginal (i.c. single variable) statistics from the dependence structure (Joe) 2014). For this reason, copula models are particularly effective for mutual information estimation (Jenison & Reale 2004 Calsaverini & Vicente (209b), which quan-
tifies how much knowing one variable reduces the uncertainty about another variable (Quiroga & Panzeri 2009). Copula models can also escape the 'curse of dimensionality' by factorising the
multi-dimensional dependence into pair-copula constructions called vines (Aas al.
2009
Czado
2010
Copula models have been successfully applied spiking activity (Onken et al. 2009: Hu et al.
2015 Shahbaba et al. 2014; Berkess et al. 2009). 2-photon calcium recordings (Safaai 2019) and multi-modal neuronal datasets (Onken Panzeri 2016). However, these models assumed that the dependence between variables was static, whereas in neuronal recordings it may be dynamic or modulated by behavioral context (Doiron et al., 2016 Shimazaki et al., 2012). Therefore, it might be helpful to explicitly model the continuous time- or context-dependent changes in the relationships between variables, which reflect changes in an underlying computation.
Here, we extend copula-basee approach by adding explicit conditional dependence to the param- eters of the copula model, approximating these latent dependencies with Gaussian Processes (GP). It was previously shown that such combination of parametric copula models with GP priors out- performs static copula models (Lopez-Paz et al.] 2013) and even dynamic copula models on many real-world datasets, including weather forecasts, geological data or stock market data (Heránndez- Lobato et al. 2013). Yet, this method has never been applied to neuronal recordings before.
In this work, increase the complexity of both marginal and copula models in order adequately describe the complex dependencies commonly observed in neuronal data. In particular, we use conditional marginal models to account for changes of the single neuron statisties and mixtures of parametric copula models to account for changes in tail dependencies. We also improve the scalabil- ity of the method by using stochastic variational inference. We develop model selection algorithms, based on the uully-Bayesian Watanabe -akiik information criterion (WAIC). Finally and most im- portantly, we demonstrate that our model suitable for estimating mutual information. It performs especially well when the parametric model can closely approximate the target distribution. When it is not the case, our copula mixture model demonstrates sufficient flexibility and provides close in- formation estimates, comparable to the best statc-of-the-arr non- -onpaaametiic information estimators.
The goal of this paper is to propose and validate the statistical Copula-GP method, and illustrate that it combines multiple desirable properties for neuroscience applications: interpretability of paramet- ric copula models, accuracy in density and information estimation and scalability to large datasets. We first introduce the copula mixture models and propose model selection algorithms (Sec. 2). We then validate our model on synthetic data and compare its performance against other commonly used information estimators (Sec. B) Next, we demonstrate the utility of the method on real neuronal and behavioral data (Sec. 4 We show that our Copula-GP method can produce bivariate models that emphasize the qualitative changes in tail dependencies and estimate mutual information that ex- poses the structure of the task without providing any explicit cues to the model. Finally, we measure information content in the whole dataset with behavioral variables and more than 100 neurons.
PARAMETRIC COPULA MIXTURES WITH GAUSSIAN PROCESS PRIORS
Our model is based on copulas: multivariate distributions with uniform marginals. Sklar's theo- rem (Sklar 1959) states that any multivariate joint distribution can be written in terms of univariate marginal distribution functions p(y.) and unique copula which characterizes the dependence struc- ture: p(y1..... . yN) = c(Fi (y1)... F' (yN)) x II p(u.). Here, E.(-) are the marginal cumulative distribution functions (CDF). Thus, for each i, ((())) is uniformly distributed on [0.1].
For high dimensional datasets (dim y), maximum likelihood estimation for copula parameters may become computationally challenging. The two-stage inference for margins (IFM) training scheme is typically used in this case (Joe) 2005). First, univariate marginals are estimated and used to map the data onto multidimensional unit cube. Second, the parameters of the copula model are inferred.
2

---

Under review as a conference paper at ICLR 2021
Goussian BE 1.1
Frank A EER
Clayton 9 [0, w)
Gumbel 0 [1. m)
Clayton + 90" Gumbel
1.0
27%
1.0
10-
0.8
0.8
0.8
08
0.6
0.6
8
0.4
0.4
0.2
02-
0.2
0.2
0.2
8 5.0
5,0
0.0
0.0
0.5
10
0.
0.5 1.0 0.0
0.5
1.0
0.0
0.5
1.0
0.0
0.5
10
ys
Mixture copulas
Rotation-equivarian copulas
x 4 rotations
= copula
Figure 1: Copula families used in the mixture models in our framework. The percentage in the upper-left corner shows how often each of the families was selected to be used a copula mixture for pairwise relationships in the real neuronal data from Pakan et al. (2018a) (see Sec. 4).
Conditional copulas Following the approach by Hernández-Lobato et al. (013),, we are using Gaussian Processes (GP) model the conditional dependencies of copula parameters:
NV II P[
ply|=) = (unl=)), - FN(UN())(Z) x
In the most general case, the marginal PDFs p(us(z) and CDFs (B (r) and the copula c(... [2) it- self can conditioned on z. our framework, æ is assumed to be one-dimensional. Gaussian Process is ideally suited for copula parametrization, as it provides an estimate of the uncertainty in model parameters, which we utilize our model selection process (Sec. 2.T
Conditional marginals Previous works conditional copulas (Hernández-Lobato
2013)
relied on the assumption that marginal distributions remain constant. This assumption might not hold for some real- life datasets, including neuronal recordings. Thus, we propose to use conditional marginal distributions instead. Note, that the interpretation of the copula model itself would depend on the assumptions made for the marginal distributions (see Appx. C.T for further discussion).
In order to estimate marginal CDFs F (():(z), we use the non-paraettric fastKDE (O'Brien et al. 2016) algorithm, which allows for direct estimation of the conditional distributions. These CDFs are then used to map the data onto unit hypercube using the probability integral transform: F F[y(z) - W; ~ 1]. such that uN is uniformly distributed for any x.
Bivariate copula families We use copula families as the building blocks for our copula models: Gaussian, Frank, Clayton and Gumbel copulas (Fig. IT. All of these families have a single parameter, corresponding to the rank correlation (Table Dp. We also use rotated variants (90', 180". 270") of Clayton and Gumbel copulas in order to express upper tail dependencies and negative correlation.
Table 1: Bivariate copula familics and their GPLink functions
Copula
Domain
GPLink(f): R dom(())
Independence
Gaussian
[-1,1]
Erf(f/1.4)
Frank
(-o0, oo)
0.3 f +iign(f) -)
Clayton
[0,oc)
Exp(0.3 f)
Gumbel
(1,00) 1 Exp
Since we are primarily focused on the analysis of neuronal data, we have first visualized the de- pendencies in calcium signal recordings after a probability integral transform, yielding empirical conditional copulas. As a distinct feature in neuronal datasets, we observed changes in tail depen- dencies with regard to the conditioning variable. Since none of the aforementioned families alone could describe such conditional dependency, we combined multiple copulas into linear mixture model: e c(u(r) = -i= Φy (I)); (w; 0g (x)), where M is the number of elements, ()(x) is the con- centration of the jthh copula in mixture, Gj is the pdf of the jth copula, and 03 is its parameter.
3

---

Under review as a conference paper at ICLR 2021
Each of the copula families includes the Independence copula a special case. To resolve this over- completeness, we add the Independence copula as separate model with zero parameters (Table [ For independent variables Yind, the Independence model will be preferred over the other models in our model selection algorithm (Sec. 2.1 l. since it has the smallest number of parameters.
Gaussian Process priors We parametrize each copula in the mixture model with an independent latent GP: f ~ Np x 1, (x æ)). For each copula family, we constructed GPLink functions (Table that map the GP variable onto the copula parameter domain: 8 GPlinke, (fz), R - dom Next, we also use GP to parametrize concentrations Φ, (z), which are defined on a simplex 1):
-- II f. im 1m= (Tm (fm+
M (M M m
, = (1 (1
ta = O,
where & is a CDF of a standard normal distribution and fm (imm x 1, 1,K,, (2, We use the RBF kernel Kx (zz, æ) with bandwidth parameter A. Therefore, the whole mixture model with M copula elements is parameterized by [2M independent GPs and requires [2M - 1] hyperparam- eters: {A} M for 0 and MA}M-1 for Φ
Approximate Inference Since our model has latent variables with GP priors and intractable poste- rior distribution, direct maximum likelihood Type-Il estimation is not possible and an approximate inference needed. Such an inference problem with copula models has previously been solved with the expectation propagation algorithm (Hermández-Lobato et al. (213), see comparison in (B4). which was not suitable for large scale data. Recently, a number of scalable approximate inference methods were developed. including stochastic variational inference (SVI) (Titsias 2009)
Cheng & Boots 2017), scalable expectation propagation (SEP) (Hernández-Loban
& Hernández- al., 2015). well as scalable exact GP (Wang et al. 2019). We chose to use SVT due to availability of the well-established GPU-
Lobato (2016), and MCMC based algorithms (Hensman
accelerated libraries: PyTorch (Paszke al. 2017) and GPy GPyTorch (Gardner al. 2018).
2. BAYESIAN MODEL SELECTION
We use the Watanabe Akaike information criterion (WAIC, Watanabe (2013) for model selection. WAIC a fully Bayesian approach estimating the Akaike information criterion (AIC) (see Eq. 31 in the original paper by Watanabe 22033).. The main advantage of the method is that it avoids the empirical estimation of the effective number of parameters, which is often used for approximation of the out-of-sample bias. It starts with the estimation of the log pointwise posterior predictive density ((ppd) (Gelman et al. 2014):
S Σ p(u: 11w.))))
1ppd Σ log
PWAIC Σ VsS [log (logp(v:)))
where {0*}s is draw from posterior distribution, which must be large enough represent the posterior. Next, the PWAIC approximates the bias correction, where V represents sample vari- ance. Therefore, the bias-corrected estimate of the log pointwise posterior predictive density is given by: elppdware lppd PWAIC = -N WAICoriginal-
In the model selection process, we aim to choose the model with the lowest WAIC. Since our copula probability densities are continuous, their values can exceed I and the resulting WAIC is typically negative. Zero WAIC corresponds to the Independence model (pdf 1 on the whole unit square).
Since the total number of combinations of 10 copula elements (Fig. considering rotations) large, exhaustive search for the optimal model is not feasible. In our framework, we propose two model algorithms for constructing close-to-optimal copula mixtures: greedy and heuristic (see Appx. (..x for details). The greedy algorithm is universal and can be used with any other copula families without adjustment, while the heuristic algorithm fine-tuned to the specific copula families used in this paper (Fig. []. Both model selection algorithms were able to select the correct 1- and 2- component model on simulated data and at least find close approximation (within WAIC = 0.005) for complex models (see validation of model selection in Appx.
4

---

Under review as a conference paper at ICLR 2021
2.2 ENTROPY AND MUTUAL INFORMATION
Our framework provides tools for efficient sampling from the conditional distribution and for cal- culating the probability density p(y|z). Therefore, for each x=t the entropy H x=t) can be estimated using Monte Carlo (MC) integration: Hy - (p( log (=))) The prob- ability p(=)x= factorizes into the conditional copula density and marginal densities (1, hence the entropy also factorizes (Jenison Reale 2004) as H (y|==()) = > H (wi(z=1) + (u (x-t), where u = Fz) The conditional entropy can be integrated as Hy]z) = Σi= ())))) + HH (t)d,, separating the entropy of the marginals [ye]x from the copula entropy.
Now, y) [ar, G(y)) if G(y) is 1) homeomorphism, 2) independent of (Kraskov et al. 2004). If marginal statistics are independent of æ, then the probability integral transform u = Fly satisfies both requirements, and 7 (xr, [(x, u). Then, in order to calculate the mutual infor- mation [ar, u) = (u)) Hu)r), we must also rewrite it using only the conditional distribution p(u|z), which is modelled with our Copula-GP model. This can done as follows:
T(a u) = H(u)) - H( ()((t)dt =
E logp(ul>) - E log E p(u|z).
(2)
p(u,a)
p(u) p(a))
The last term in 2 involves nested integration, which is computationally difficult and does not scale well with N dim Therefore, we propose an alternative way of estimating I(x, y), which avoids double integration and allows us to use the marginals conditioned on æ (u² = F(y x)), providing better estimate of Hylzz). We can two separate copula models, for estimating p(y) and calculating H Hy), and another one for estimating p(y (ar) and calculating H(y (r):
N
((r,y) = Σ Mi, vi) + (ai,........... UN) / Ho(wg..... uN 1s = ()(())dt, 1=1
(3)
where both entropy terms are estimated with MC integration. Here only integrate over the unit cube [0, 1] and then dom x. whereas 2 required integration over [0, 1| x dom x.
The performance of both 2 and critically depends on the approximation of the dependence structure, i.c. how well the parametric copula approximates the true copula probability density. If the joint distribution p(y1 ... yNv has a complex dependence structure, as we will see in synthetic examples, then the mixture of parametric copulas may provide poor approximation of p(y) and overestimate H.(wi,.... u.), thereby overestimating (xx,y), The direct integration (2) on the other hand, typically underestimates the [xx, y) due to imperfect approximation of p (z), and under assumption that the marginals can be considered independent of x..
We further refer to the direct integration approach 2 as "Copula-GP integrated" and to the alterna- tive approach 3 as "Copula-GP estimated" and assess both of them on synthetic and real data.
2.3 COPULA VINE CONSTRUCTIONS
High-dimensional copulas can be constructed from bivariate copulas by organizing them into hierar- chical structures called copula vines (Aas et al. 2009). In this paper, we focus on the canonical vine or C-vine, which factorizes the high dimenssina copula probability density function as follows:
N
N
c(u) II ca (11, Uz)
x
II
II 1=11
ey (K)xe (F(uel F(u,((m))<<i))
(4)
(=2
where flitkei = 1,...,, and F(.l.) is a conditional CDF. Note, that all of the copulas in (4 also be conditioned via Copula- GP model. We choose the first variable in to be the one with the highest rank correlation with the rest (sum of absolute values of pairwise Kendall's 7), and condition all variables on the first one. We repeat the procedure until variable is left (see Appx. It was shown by Czado et al. (202)) that this ordering facilitates C-vine modeling.
Code availability Code will be made available on GitHub upon paper acceptance.
5

---

Under review as a conference paper at ICLR 2021
3 VALIDATION ON ARTIFICIAL DATA
We compare our method with the other commonly used non-parametric algorithms for mutual infor- mation estimation: KSG by Gao et al. (BI-KSG Gao et al. (2016)) and the Mutual Information Neural Estimator relatively low dimensional data (< for which we can still directly calculate the true mutual information.
(KSG, Kraskov al. (2004)) Bias-Imaroo Improved--
(MINE, Belghazi et al. (2018). ) In this section,
A Gaussian with carr= Al.
B StudenT with corr=0.7, df=2/150
Gauss carr= parr==0.1.1.0, ý1 = ( 1
p(u)
p(
p(u)
p(ujk=0)
plujx=1)
p(u)
p(ujx=0)
[x=xx
0-
of
150
Capula-GP integr.
Copula estim.
MINE100
MINE200
100
0.2
Hly
1[x-2
dip
y
0.0
8
8 #variables (dim y)
6 #variables (dim y)
10
#varlables (dim y)
Figure 2: Conditional entropy H (y[zr) and mutual information (xx, measured by different meth- ods on synthetic data. Upper row shows the dependency structures p(u) and conditional dependency structures at the beginning and the end of the dom = [0. 1]. A Multivariate Gaussian samples. Multivariate Student T samples. C Multivariate Gaussian samples y (same as A), morphed into another distribution p(ý) with a tail dependence, while (xx, y) = Tx, 8). Gray intervals show either standard error of mean (SE, 5 repetitions), or ((SE)² + (M (tot)) for integrated variables. Note, that MINE estimates depend on the choice of hyper parameters (e.g. number of hidden units).
First, if there is no model mismatch, we expect our information estimates to be unbiased due to the MC estimator being unbiased (Luengo et al., 2020p. We confirm this fact using the dataset sampled from a multivariate Gaussian distribution, with cov 31,) Kronecker's delta and p = -0.1 + 1.1z, x E 1]. Our algorithm selects a Gaussian copula on these data, which perfectly matches the true distribution. The same applies to any linear mixture of copulas from Tab. as well (see Appx. B.3). and is covered by the automated tests on simulated data (see Appx. Ap.).[B.Ip. Therefore, Copula- GP measures both entropy and mutual information exactly (within integration tolerance, see Fig. 2 The performance of the non-parametric methods on this dataset is lower. It was shown before that KSG and MINE both severely underestimate the MT for high-ding sech-imennsioall Gaussians with high correlation (e.g. see Fig. 1 in Belghazi et al. (2018)). The Copula-GP model (integrated) provides accurate estimates for highly correlated (up p 0.999, at least up to 20D) Gaussian distributions (see Appx. B.2).
(1 p) dzj. where 6,3
Next, we test the Copula-GP performance on the Student T distribution, which can only be approx- imated by our copula mixtures, but would not exactly match any of the parametric copula families from Fig. []] We keep the correlation coefficient p fixed 0.7, and only change the number of de- grees of freedom exponentially: df = exp(5x:) 1. E [0. 1]. This makes the dataset particularly challenging, as all of the mutual information T(x, encoded in tail dependencies of ((y])r). The (2) of the Student T distribution calculated analytically (see Eq. A.12 in Calsaverini & Vicente (2009a)) and I [x, y) was integrated numerically according to 2 given the true p(y|z).
true H [y(z)
Fig. 3 shows that most of the methods underestimate y(z,y). Copula-GP (integrated) and MINE (with 100 hidden units) provide the closest estimates. The training curve for MINE with more hidden units (200,500) showed signs of overfitting (abrupt changes in loss at certain permutations) and the resulting estimate was higher than the true [xx,y) y) at higher dimensions. was shown before that
MINE provides inaccurate and inconsistent results datasets with low 7 yx, (Song & Ermon

---

Under review as a conference paper at ICLR 2021
We also demonstrate [2x y) estimation with combination of two copula models for (ly) and H [y[z): "Copula-GP estimated" (see Eq. 3). In lower dimensions, it captures less information than "Copula-Gi integrated"" but starts overestimating the true MI at higher dimensions, when the inaceuracy of the density estimation for p(y) builds up. This shows the limitation the "estimated" method, which can either underestimate overestimate the correct value due parametric model mismatch, whereas "integrated" method consistently underestimates the correct value. We conclude that Copula-GP and MINE demonstrate similar performance in this example, while KSG-based methods significantly underestimate [zz, y) in higher dimensions.
2019).
Finally, we created another artificial dataset that is not related to any of the copula models used in our framework (Table D. We achieved that by applying a homeomorphic transformation F(y) to a multivariate Gaussian distribution. Since the transformation
independent of the conditioning variable, it does not change the y) [xx F(y)) (Kraskov et al. 2004). Therefore, we possess the true T (1r, y), which is the same as for the first example Fig. Note, however, that there is no ground truth for the conditional entropy in this example, since 17 (F())). We transform
1/)//N and again transform the marginals using the empirical probability integral transform u = Fyy. Both conditional p(u zz) and unconditional p(u) densities here do not match any of the parametric copulas from Table As result, "Copula-GP estimated" overestimated the correct value, while "Copula- GP integrated" underestimated it similarly to the MINE estimator with 100 hidden units.
the Gaussian copula samples y E "fo.il from the first example as 31 3/s + [[]) N
Fig. demonstrates that the performance of the parametric Copula-GP model critically depends on the match between the true probability density and the best mixture of parametric copula ele- ments. When the parametric distribution matches the true distribution (e.g. Fig. 2 or Fig. A5 our Copula-GP framework provides unbiased estimates and predictably outperforms all non-parametric methods. Nonetheless, even when the exact reconstruction of the density is not possible (Figs. 2 C), the mixtures of the copula models are still able to model the changes in tail dependencies, least qualitatively. In those challenging examples, our method performs similarly to the neural-network based method (MINE) and still outperforms KSG-like methods.
VALIDATION ON REAL DATA
X
60
120
1440
160 cm
Neuron
y
Nouron
Reward
Neuron Licking rate
time
O-
60-120
120-140
140-160 160
E
$9.%%
0.2
HING, N631
20-
H(u [x]
0.1
15-
0.0
Neuran
Neuron
Neuren
Ncuron
F
Position,
Position
Nouron 60
90%
58%
0.2
/LLicks, N60)
Copula-GP int,jest.
MINE
0.1
0.0
Licks
Licks
Licks
Licks
Posttio,
variables (subset u*)
Figure 3: Applications of the Copula-GP framework to neuronal and behavioral data from the visual cortex of awake mice. Schematic of the experimental task (Pakan et al.) 2018a Henschke et al. 2020) in virtual reality (VR); B Example traces from ten example trials: is a vector of neuronal recordings (blue) and behavioral variables (red); C-D Density plots for: the noise correlation (C) and the behavioral modulation (D) examples; E-G Conditional entropy for the bivariate examples (E-F) and the population statistics (G); H Comparison of Copula-GP methods ("estimated" and "integrated") vs. non-parametric MI estimators on subsets of variables.
position in VR, y
We investigate the dependencies observed in neuronal and behavioral data and showcase possible applications of the Copula-GP framework. We used two-photon calcium imaging data of neuronal
7

---

Under review as a conference paper at ICLR 2021
population activity in the primary visual cortex of mice engaged in a visuospatial navigation task in virtual reality (data from Henschke et al. (2020). ). Briefly, the mice learned to run through a virtual corridor with vertical gratings the walls (Fig. E3. 0-120 cm) until they reached reward zone (Fig. A 120-144 cm), where they could get a reward by licking a reward spout. We condition our Copula-GP model the position in the virtual environment x and studied the joint distribution of the behavioral (ii... .....)).)))))))) Wis and neuronal (26... (109) variables (dim y=19). Fig. [ shows part of the dataset (trials 25-35 out of 130). The traces here demonstrate changes in the position æ of the mouse as well as the activity of 3 selected neurons and the licking rate. These variables have different patterns of activity depending on and different signal-to to-lltto-oise ratios. Both differences are reflected in marginal statistics, which are shown on the right with the density plots of equal area.
Constructing interpretable bivariate models We first studied bivariate relationships between neurons. In order to do this, we transformed the raw signals (shown in Fig. (3)) with a probability integral transform c(u(z)) subject to the position in the virtual reality x and related visual information (Fig. 3C) Such stimulus-related changes in the joint variability of two neuronal signals are commonly described as noise correlations. The Copula-GP model provides a more detailed description of the joint proba- bility that goes beyond linear correlation analysis. In this example, the dependence structure is best characterized by combination of Gaussian and Clayton copula (rotated by 90%). The density plots Fig. BC demonstrate the match between the true density (outlines) and the copula model density (blue shades) for each part of the task. We measure the accuracy of the density estimation with the proportion of variance explained R?", which shows how much of the variance of the variable y/2 can be predicted given the variable y1 (see Eq. in Appx. A.1). The average R² for all y is provided in the upper right corner of the density plots.
F(y). We observed strong non-trivial changes in the dependence structure
Next, we show that our model can be applied not only to the neuronal data, but also to any of the behavioral variables. Fig. BD shows the dependence structure between one of the neurons and the licking rate. The best selected mixture model here is Frank + Clayton 02 + Gumbel 270", which again provides accurate estimate of the conditional dependence between the variables. Such copula mixture models allow us to analyze the contribution of the tail dependencies to the mutual information, and backtrack the outliers in order to relate them to certain trials or behaviors (see Appx. C.2). Therefore, Figs. 3DD demonstrate that our Copula-GP model provides both an accurate fit for the probability distribution and an interpretable visualization of the dependence structure.
Figs. 3 F show the absolute value of the conditional entropy |H (u"|==)) which is equivalent to the mutual information between two variables [(u uz)) For both examples, the MI peaks in the reward zone. The bivariate Copula-GP models were agnostic of the reward mechanism in this task, yet they revealed the position of the reward zone as an anomaly in the mutual information.
Measuring information content in a large neuronal population Finally, we constructed a C- vine describing the distribution between all neuronal and behavioral variables (dim u² 109) and measured the conditional entropy H Huu (x) for all variables in the dataset {uf_..uffo9l The conditional entropy in Fig. 3G peaks in the reward zone (similarly Figs. 3--F) and also at the beginning of the trial, where the velocity of the animal varied the most on the trial-to-trial basis.
While constructing the C-vine, we ordered the variables according to their pairwise rank correlations (see Sec. (23)). We considered subsets of the first N variables and measured the MI with the position for each subset. We compared the performance of our Copula-GP method on these subsets of u² vs. KSG and MINE. Fig. 3 shows that all 3 methods provide similar results on subsets of up to 10 variables, yet in higher dimensions both MINE and KSG show smaller T(x (uā²N)) compared to our Copula-GP method, which agrees with the results obtained on the synthetic data (Fig. ( The true values of (x, {) are unknown, yet we expect the integrated Copula-GP (solid line) to underestimate the true value due to parametric model mismatch. The Copula-GP "estimated" (dashed line) can either under- or over over-estimate the information (see eq.3), but here it almost per- fectly matches the "integrated" result, which suggests that the model was able to tightly approximate both p(u" (x) and p(u"), and, as a result (cq.3). [(x, x H). These results demonstrate superior performance of our Copula-CH GP model on high-dimensional neuronal data.
8

---

Under review as a conference paper at ICLR 2021
5 DISCUSSION
We have developed a Copula-GP framework for modeling conditional multivariate joint distribu- tions. The method based on linear mixtures of parametric copulas, which provide flexibility for estimating complex dependencies but still have interpretable parameters. We approximate condi- tional dependencies of the model parameters with Gaussian Processes which allow us to implement Bayesian model selection procedure. The selected models combine the accuracy in density es- timation with the interpretability of parametric copula models. We demonstrated the performance of our framework in mutual information estimation synthetically generated data. Despite the general limitations of the parametric models, our framework performs either significantly better (when the true dependence matches parametric model), at least as good as the state-of-the-an non-parametric information estimators. Unlike black-box machine-learning methods for entropy es- timation, our Copula-GP model uses interpretable elements (pair copulas) and latent GP parameters, which can be isolated, visualized, and analyzed (Appx. C The framework is also well suited for describing the dependencies observed in neuronal and behavioral data. We demonstrated that the model scales well at least to 109 variables and 21k samples, while theoretically. the parameter inference scales as O(n m²), where n a number of samples and m is the (effective) number of variables (see Appx. Appx.A.6). Alternative popular non-copula methods that could be applied to our real data, include GPFA or GLMs on the deconvolved spikes. Contrary to these approaches, copula methods allow to model the dependencies between elements with utterly different statistics (e.g. licks vs. velocity). In addition, Copula-Gr GP explicitly represents the dependencies as a function of position, revealing insightful information about the task structure. To the best of our knowledge, there are currently no other methods with this combination of features. The possible applications of Copula-GP include, but not limited to, studying the neuronal population statistics, noise cor- relations, behavioral or contextual modulations. In summary, we demonstrated that the Copula-GP approach can make stochastic relationships explicit and generate accurate models of dependencies between neuronal responses, sensory stimuli, and behavior. Future work will focus on implement- ing model selection for the vine structure and improving the scalability of the mutual information estimation algorithm.
REFERENCES
Kjersti Aas, Claudia Czado, Amoldo Frigessi, and Henrik Bakken. Pair-copula constructions of multiple de- pendence. Insurance: Mathematics and economics, (44(2)182-199) 2009.
Mohamed Ishmael Belghazi, Aristide Baratin, Sai Rajeswar, Sherjil Ozair, Yoshua Bengio, Aaron Courville, and Devon Hjelm. MINE: Mutual Information Neural Estimation, 2018
Pietro Berkes, Frank Wood, and Jonathan W Pillow. Characterizing neural dependencies with copula models. In Advances in neural information processing systems, pp. 129-136, 2009.
R. S. Calsaverini and R. Vicente. An information-theoorticc approach to statistical dependence: Copula infor- mation. EPL (Europhysics Letters), 886))68003, Dee 2009a. ISSN 1286-4854. doi: 10. 209902995-507 88/68003. URL (tRt:t/t::///at http://............./////////.//
Rafael Calsaverini and Renato Vicente. An information- approach to statistical dependence: Copula information. EPL (Europhysics Letters), 88(6):68003, 2009b.
Ching-An An Cheng and Byron Boots. Variational inference for gaussian process models with lincar complexity. In Advances in Neural Information Processing Systems, pp. 5184 5194, 2017.
Claudia Czado. Pair-copula constructions of multivariate copulas. In Copula theory and its applications, pp. 93-109. Springer, 2010.
Claudia Czado, Ulf Schepsmeier, and Aleksey Min. Maximum likelihood estimation of mixed c-vines with application exchange rates. Statistical Modelling, 12(3):229-255 2012.
Brent Doiron, Ashok Litwin-Ke Kumin, Robert Rosenbaum, Gabriel K Ocker, and Krešimir Josié, The mechanics of state-dependent neural correlations. Nature neuroscience, 19(3):383, 2016.
Daniel A Dombeck, Anton N Khabbaz, Forrest Collman, Thomas I. Adelman, and David W Tank. Imaging large-scale neural activity with cellular resolution in awake, mobile mice. Newron, 56(1):43557 2007.

---

Under review as a conference paper at ICLR 2021
Weihao Gao, Sewoong Oh, and Pramod Viswanath. Demystifying Fixed k-Nearest Neighbor Information Estimators, 2016.
Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. GPyTorch: Black- box matrix-matrix Gaussian process inference with GPU acceleration. Advances in Neural Information Processing Systems, pp. 7576-7586, 2018.
Andrew Gelman, Jessica Hwang. and Aki Vehtari. Understanding predictive information criteria for bayesian models. Statistics and computing, 24(6):997-1016, 2014.
Christine Grienberger, Xiaowei Chen, and Arthur Konnerth. Dendritic function in vivo. Trends in xciences, 38(1):45-54, 2015.
Fritjof Helmehen. Two-photon functional imaging of neuronal activity. CRC Press, 2009.
Julia U Henschke, Evelyn Dylda, Danai Katsanevaki, Nathalie Dupuy, Stephen P Currie, Theoklitos Amvrosiadis, Janelle MP Pakan, and Nathalie L Rochefort. Reward association enhances siimulus-ppccine representations in primary visual cortex. Current Biology, 2020.
James Hensman, Alexander G Matthews, Maurizio Filippone, and Zoubin Ghahramani. MCMC for varia- tionally sparse gaussian processes. In Advances in Neural Information Processing Systems, pp. 1648-1656, 2015.
Daniel Hermández-Lobenn and José Miguel Hernández-Lobato. Scalable gaussian process classification via expectation propagation. Artificial Intelligence and Statistics, pp. 168-176, 2016.
José Miguel Hernández-Lobater James R Lloyd, and Daniel Hernnndez ditional copulas with applications financial time series. C. C. Burges, Bottou, Welling,
Gaussian process
Z. Ghahramani, and K. Q. Weinberger (eds.), Advances in Neural Information Processing Systems
26, pp. 1736-1744. Curran Associates, Inc., 2013. URL. httt:://paperae http:/pppprraa.. wtwpppppccpppprr.///pp wit th apph www.thhaapalaaaiinsesseeeee
5084
t ticnal coll-coopllee as aa-wwa--eee
-time-se series.
pdf
Meng Hu, Kelsey L Clark, Xiajing Gong, Behrad Noudoost, Mingyao Li, Tirin Moore, and Hualou Liang. Copula regression analysis of simultaneously recorded frontal eye field and inferotemporal spiking activity during object-based working memory. Journal of Newroscience, 35(23):8745 8757, 2015.
Robin AA Ince. Riccardo Senatore, Ehsan Arabzadeh, Fernando Montani, Mathew E Diamond, and Stefano Panzeri. Information-thewee theoretic methods for studying population codes. Neural Networks, 23(6):713-727, 2010.
Rick L. Jenison and Richard
Reale. The shape of neural dependence. Neural Computation, 16 (4):665-672 672, 2004, doi: 00116220899766664222966666 08997660432286065 322860659
URL. https::/doi.orrg 1162/
Harry Joc. Asymptotic efficiency of two-stage estimation method for copula-based models. Jour-
nal of Multivariate Analysis, 94(2):401 - 419, 2005. ISSN 0047-259X. doi: http://doi.rg//00
j016//jivva.2004.06
URL http: //t:://wwwww wew.nceiieeeet..e..c....... com/ http://ccccce/rrice/e......
S004 7259XX40091289
Harry Joc, Dependence modeling with copulas. CRC press, 2014.
James J Jun, Nicholas Steinmetz, Joshua Siegle, Daniel Denman, Marius Bauza, Brian Barbarits, bert K Lec, Costas A Anastassiou, Alexandru Andrei, Çağatay Aydm, et al. Fully integrated silicon probes for high-density recording of neural activity. Nature, 551(7679):232225) 236, 2017.
Adam Kohn, Ruben Coen-Cagli, Ingmar Kanitscheider, and Alexandre Pouget. Correlations and neuronal population information. Annual review of neuroscience, 39:237-256, 2016.
Alexander Kraskov, Harald Stögbauer, and Peter Grassberger. Estimating mutual information. Phys. Rev. E. 19:066118, Jun 2004, doi: 10.1 10.100/Pyyeeeee.6..0666 69.661.8.... hRI.ttss: https:///inkk aps .rr./org//oo. brg/doi/10
1103/PhysRevE. 69. 066138
David Lopez-Paz. Jose Miguel Hernández-Lobato, and Ghahramani Zoubin. Gaussian process vine copulas for multivariate dependence. In International Conference on Machine Learning, pp. 10-18, 2013.
David Luengo, Luca Martino, Mónica Bugallo, Víctor Elvira, and Simo Särkkä. A survey of monte carlo methods for parameter estimation. EURASIP Journal on Advances in Signal Processing, 2020(1):1-62, 2020.
10

---

Under review as a conference paper at ICLR 2021
Alexander Mathis, Pranay Mamidanna, Kevin M Cury, Taiga Abe, Venkatesh N Murthy, Mackenzie Weygandt Mathis, and Matthias Bethge. DacpLabCut: markerless pose estimation of user-defined body parts with deep learning. Nature neuroscience, (1) 1e:::12112211 2018.
Arno Onken Stefano Panzeri. Mixed vine copulas as joint models spike counts local field potentials. In Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS' 16, pp. 1333-1341, Red Hook. NY, USA. 2016. Curran Associates Inc. ISBN 978 1711588888819
Arno Onken, Steffen Grünewälder, Matthias HJ Munk, and Klaus Obermayer. Analyzing short-term noise dependencies of spike-count counts in macaque prefrontal cortex using copulas and the flashlight transformation. PLoS computational biology, 5(1 2009.
Travis A. OrBie, Karthik Kashinath, Nicholas R. Cavanaugh, William D. Collins, and John O'Brien. A
fast and objective multidimensional kernel density estimation method: fastKDE. Computational Statistics & Data Analysis, 101:148- 160, 2016. ISSN 0167-9473, doi: https://dii.rr/10000010006 cxda 2016.02..1 URL http: sc encedi rect com science /ttp:///ccc.eee........//////1//
Janelle MP Pakan, Stephen P Currie, Lukas Fischer, and Nathalie I. Rochefort. The impact of visual cues, reward, and motor feedback on the representation of behaviorally relevant spatial locations in primary visual cortex. Cell reports, (4110)252212522 2018a.
Janelle MP Pakan, Valerio Francioni, and Nathalie L Rochefort. Action and learning shape the activity of neuronal circuits in the visual cortex. Current opinion in neurobiology, 52:88 -78889977777 2018b.
Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang. Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer. Automatic differentiation PyTorch. In Advances in Neural Information Processing Syxtems, 2017.
Rodrigo Quian Quiroga and Stefano Panzeri. Extracting information from neuronal populations: information theory and decoding approaches. Nature Reviews Neuroscience, 173-18173-185, 2009.
Wang A Panzeri S Harvey CD Safaai, H. Characterizing information processing of parietal cortex projections using vine copulas. In Bernstein Conference 2019. American Physical Society, 2019. nodee org. attpg//bsstatts////
doi: 10ee6333-8fff--2000
URL htttp::/abst a27 2bbb96211a
Shreya Saxena and John Cunningham. Towards the neural population doctrine, 2019. ISSN 18736882.
Babak Shahbaba, Bo Zhou, Shiwei Lan, Hernando Ombao, David Moorman, and Sam Behseta. A semipara- metrie bayesian model for detecting synchrony among multiple neurons. Neural Computation, 26(9):2025 2051, 2014. doi: 10. 1021/CO/) JaAANECOO_aa,00000 00631.. URL pattpa:///oon org/10.11/.ooo www.11111111...... PMID: 24922500.
Maoz Shamir and Haim Sompolinsky. Nonlinear population codes. Neural computation, 16(6):1105-1133: 2004.
Hideaki Shimazaki, Shun-ichi Amari, Emery N Brown, and Sonja Grün State-space analysis of time-varying higher-order spike correlation for multiple neural spike train data. PLoS computational biology. 8(3), 2012.
Abe Sklar. Fonctions de reprtition dimensions et leursmarges. Publ. Inst. Statis. Univ. Paris. 8:229-231, 1959.
Jiaming Song and Stefano Ermon. Understanding the limitations of variational mutual information estimators, 2019.
Benjamin Staude, Stefan Rotter and Sonja Grün. Cubic: cumulant based inference of higher-order correlations massively parallel spike trains. Journal of computational neuroscience. 29(1-2):3273350 2010.
Ian Stevenson and Konrad P. Kording. Neuroscience, 14(2):139 142, feb 2011. ISSN 1097-6256. doi: 10. 10.1033/m..731. URL http:// hett:://wwk
How advances in neural recording affect data analysis. Nature
nature. ccm art i c as/nn.22311
Carsen Stringer. Marius Pachitariu, Nicholas Steinmetz, Charu Bai Reddy, Matteo Carandini, and Kenneth D Harris. Spontaneous behaviors drive multidimensional, brainwide activity. Science, 36446443))eeay78999 2019.
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial Intelli- gence and Statistics, pp. 567-574, 2009.
11

---

Under review as a conference paper at ICLR 2021
Ke Wang. Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact gaussian processes on a million data points. Advances Neural Information Processing Systems, pp.
2019.
Sumio Watanahe. A widely applicable Bayesian information criterion. Journal of Machine Learning Research, 840Mar)88778877 2013.
A APPENDIX METHODS
A.1
GOODNESS-OF-FI
measure the accuracy of the density estimation with the proportion of variance explained R². We com-
pare the empirical conditional CDF ecdf(u2luu [n = y) vs. estimated conditional CDF ccdf(ugluu = y) and calculate:
[eed(uug/us hu
codf( v) ecdf(uz2uu = y) T2
= 1
(AI)
12
where 172 (w) quantifies the portion of the total variance of u2 that our copula model can explain given y. and u5- F(v2) = 0.5. The sum was calculated for 1/2 = 0.05m, n = 0... 20.
Next, we select all of the samples from certain interval of the task (z E [I1, æ2)) matching one of those shown in Figure in the paper. We split these samples E [0, 1] into 20 equally sized bins: {1:) For each bin 4, we calculate [AT]. We evaluate eod u2 |M y1 ccdf(uglu) E I() using copula model from the center of mass of the considered interval of x: Fy = mean(z) for samples x E [x1, We use the average measure:
E p(u] C4,)
R² (mean(u I:)),
(A2)
to characterize the goodness of fit for a bivariate copula model. Since un is uniformly distributed on [0, 1]. the probabilities each bin p(u E 4:) are equal to 1/20, and resulting measure H² just an average H²² from all bins. The results were largely insensitive to the number of bins (e.g. 20 vs. 100).
A.2 VARIATIONAL INFERENCE
Since our model has latent variables with GP priors and intractable posterior distribution, the direct maximum likelihood Type-Il estimation is not possible and an approximate inference is needed, We used stochastic variational inference (SVI) with single evidence lower bound (Hensman al. 2015):
N
LELBO
R ((A)
[log [i]] (T]g [p(u)],
(A3)
implemented as Variatioonaie in GPyTorch (Gardner et al. 2018). Here N is the number of data samples, u are the inducing points, q(u) is the variational distribution and q(f) = p(f)u)y(u)duu
Following the Wilson and Nickisch (2015) approach (KISS-GP), we then constrain the inducing points to a regular grid, which applies deterministic relationship between and u. As a result, we only need to infer the variational distribution q(u), but not the positions of u. The number of grid points is one of the model paper-parameters:: grid_ssi dizie.
Equation A3 enables joint optimization of the GP hyper-parameters (constant mean μı and two kernel parame- ters: scale and bandwidth) and parameters the variational distribution q (mean and covariance at the inducing points: u ~ N M(H, x 1, Σ.)) ) (Hensman et al. 22055). We have empirically discovered by studying the conver- gence synthetic data, that the best results are achieved when the learning rate for the GP hyper-parameters (base.) 1r) is much greater than the learning rate for the variational distribution parameters (var_ _rr Table
For both the neuronal and the synthetic data, we use a standard normal prior p(u) ~ N/ I) for a variational distribution. Note, that the parametrization for mixture models was chosen such that the afore- mentioned choice of the variational distribution prior with zero mean corresponds to priori equal mixing coefficients d, = 1/M for 1. M. In experiments with the simulated and real neuronal data, we observed that the GP hyper-paramett optimisation problem often had minima (which is common situation, see Figure 5.5 on page 116 in Williams and Rasmussen (2006). One of those corresponds to a short kernel lengthscale (A) and low noise (ming ")), which we interpret as overfitting. To prevent overfitting, we used À ~ N (0.5, 0.2) prior on RBF kernel lengthscale parameter that allows the optimizer to approach the minima from the region of higher X, ending up in the minimum with larger lengthscale.
Priors
12

---

Under review as a conference paper at ICLR 2021
Optimization We use the Adam optimizer with two learning rates for GP hyper-parameteess (base.1r) and variational distribution parameters (v)r_1r). We monitor the loss (averaged over 50 steps) and its changes in the last 50 steps: A loss mean (loss [ :
mean (loss [-50: }) If the change be- comes smaller than chock_wa.e then we evaluate the model WAIC and check if it is lower than - -WAICzot- If it is higher, we consider that either the variables are independent, or the model does not match the data. Either way, this indicates that further optimisation counterpraductive. If the WAIC the optimisation until the change loss 50 steps Aloss becomes smaller than 1 csss coss.tel 1 (see Table Al].
WATCCats we proceed with
Effective learning rates for different families The coefficients in the GPLink functions for different copula familics are also a part of model hyper-parameters. The choice of these coefficients affects the gradients of the log probability function. Since GPLink functions are nonlinear, they affect the gradients in various parameter ranges to a different extent. This results in variable convergence rates depending on the true copula parameters.
To address the problem of setting up these hyper-parameters, have created the tests on synthetic data with different copula parameters. Using these tests, we manually adjusted these syper-parameters such that the GP parameter inference converged in around 1000-2000 iterations for every copula family and parameter range. We have also multiplied the GP values corresponding to the mixture coefficients by 0.5, to effectively slow down the learning of the mixture coefficients ø compared to the copula coefficients 0, which also facilitates the convergence.
Hyper-parane parameter selection The hyper-parameten, of our model (Table Al were manually tuned, often considering the trade off between model accuracy and evaluation time. A more systematic hyper-parameter search might yield improved results and better determinc the limits of model accuracy.
Table Al: Hyper-parameters of the bivariate Copula-GP model
Haper-prr
Value
Description
base_lr
10-2
Learning rate for GP parameters
var.lr
10-3
Learning rate for variational distribution
grid.size
128
Number of inducing points for KISS-GP
waic_tor
0.005
Tolerance for WATC estimation
loss
10-4
Loss tolerance that indicates the convergence
check_wa
0.005
Loss tolerance when we check WAIC
and GPLink parameters listed in Table 1.
A.3 BAYESIAN MODEL SELECTION
In model selection, we are aiming to construct a model with the lowest possible WAIC. Since our copula prob- ability densities continuous, their values exceed I and the resulting WAIC is typically negative. Zero WAIC corresponds the Independence model (pdf 1 on the whole unit square). We also set up tolerance (WAICtw = 0.005), and models with WATC e [-WAICrod; WAICtor] are considered indistinguishable from the independence model.
Since the total number of combinations of 10 copula elements (Fig.1) is large. exhaustive search for the optimal model not feasible. In our framework, we propose model algorithms for constructing close-to-optimal copula mixtures: greedy and heuristic.
A.4 MODEL SELECTION ALGORIHTM
The greedy algorithm (Algorithm starts by comparing WAIC of all possible single-copula models (from Table 1, in all rotations) and selecting the model with the lowest WAIC. After that, we add one more copula (from another family or in another rotation) to the first selected copula, and prepend the element that yields the lowest WAIC of the mixture. We repeat the process until the WAIC stops decrcasing. After the best model is selected, we remove the inessential elements using the reduce (.) function. This function removes those elements which have average concentration of 10% everywhere on [0, 1]. This step is added improve the interpretability of the models and computation time for entropy estimation (at a small accuracy cost) and can, in principle, be omitted.
The greedy algorithm can be improved by adding model reduction after each attempt to add an element. In this case, the number of elements can increase and decrease multiple times during the model selection process, which also must be terminated if the algorithm returns the previously observed solution. Even though it
13

---

Under review as a conference paper at ICLR 2021
Algorithm 1: Greedy algorithm for copula mixture selection
M, Mold + []:[];
S. + [Independence, Gauss, Frank, 4 x Clayton, x Gumbel]; // 4x includes all rotations
// while every update of the model yields new best while WAIC(M) WAIC(Maa) and size(S.) do
Mold t M
select from S. such that WAIC(prepend(e, M)) is minimal; . + prepend(c, M);
remove from Se:
end
Mbett f- Ceduce(MMui))
10 return Moesti
Algorithm 2: Heuristic algorithm for copula mixture selection
G <- [Gauss];
if WAIC(G) > -waicttol then
return [Iddependence]:
end
Mct < Independence, Gauss, 4 x Clayton]; MG [Ideppedencce, Gauss, Gumbel]; Mhest, Maorst + (Mcr.Mo.) sorted by WAIC; WAIC(G)) WAIC(MMeul) then
return G; 10 end
for i x- 3... size(Mznn))) (Meez) do
M +- Mhea with i-th element replaced by Mworst [i];
12
Maest +- M if WAIC(M)) < WAIC(M
end
15 Mbest +- reduce(MMea):
if Gauss E Mbest then
17
M +- Mhes with Gauss replaced by Frank;
Moest < M if WAIC(M) WAIC(Moas):
end
Gauss often gets confused with pairs e.g. Clayton 0 Gumbel ()
20 sizee
then
for +- 1 1... (size Mbeer)
do for j t- (i + 1) ... (Mbur) do
22
M Mhea with i-th and j-th elements removed;
24
M < prepend Gauss, M);
if WAIC(M) < WAIC(M then
Mbeat + M;
break:
end
end end
30
31 end
32 Mbest +- reduce(MMaa):
33 return Moest
14

---

Under review as a conference paper at ICLR 2021
complicates the algorithm, it reduces the maximal execution time (observed on the real neuronal data) from ~90 minutes down <040 minutes.
The hcuristic algorithm focuses on the tail dependencies (Algorithm 2. First, we try a single Gaussian copula. If variables are not independent, we next compare 2 combinations of elements, which are organized as follows: an Independence copula together with Gaussian copula and either 4 Clayton or 4 Gumbel copulas in all 4 rotations (0", 90", 180", 270%). We select the combination with the lowest WAIC. After that, we take the remaining Clayton/Gumbel copulas one by one and attempt switch the copula type (Clayton to Gumbel or vise versa). If this switching decreases the WAIC, we keep a better copula type for that rotation and proceed to the next element.
Here we make the assumption, that because Clayton and Gumbel copulas have most of the probability density concentrated in one corner of the unit square (the heavy tail), we can choose the best model for each of the 4 corners independently. When the best combination of Clayton/Gumbel copulas selected, we can (optionally) reduce the model.
We have not yet used a Frank copula in heuristic algorithm. We attempt to substitute the Gaussian copula with a Frank copula (if it is still part of the reduced mixture, see lines 16-19 in Alg. 2. Sometimes, Gaussian copula can be mistakenly modeled as Clayton & Gumbel or two Gumbel copulas. So, as final step (lines 20-31, Alg. 2). we select all pairwise combinations the remaining elements, and attempt to substitute each of the pairs with a Gaussian copula, selecting the model with the lowest WAIC. Despite a large number of steps in this algorithm, the selection process takes only up to 25 minutes (in case all elements in all rotations are required).
The procedure was designed after observing the model selection process on variety of synthetic and real neuronal datasets.
A.5
VINE COPULAS
Vine models provide a way to factorize the high- nigh-dimensional copula probability density into hierarchical set of bivariate copulas (Aas et al.] 2009). There are many possible decompositions based on diffcrent assumptions about conditional independence of specific elements in a model, which can be classified using graphical models called regular vines (Bedford and Cooke 2001
2002) A regular vine can be represented using hierarchical conditional distribution function (e.g. F F(ua (ux) )) and each edge corresponds to bivariate copula (e.g. e(u2, Ua |ua)) The copula models from the lower trees are used obtain new conditional distributions (new nodes) with additional conditional dependencies for the higher trees, e.g. a codf of a copula c(u,, uz| (s/wu) and a marginal conditional distribution F(uz (us) from the 1st tree provide a new conditional distribution F(us|uu Mi, m22) for a 2nd tree. Therefore, bivariate copula parameters are estimated sequentially, starting from the lowest and moving up the hierarchy. The total number of edges in all trees (= the number of bivariate copula models) for an m-dimmensional regular vine equals m mm - 1)/2.
set of trees, where each node corresponds
The regular vines often assume that the conditional copulas c(u,, uyl(uz) themselves are independent of their conditioning variables [ ]. but depend on the them indirectly through the conditional distribution functions (nodes) (Acar et al. 2012)) This is known the simplifying assumption for vine copulas (Haff et al. 2010). which, if applicable, allows to escape the curse of dimensionality in nigh-dimensional copula construction.
this study, we focus on the canonical vine or C-vine, which has a unique node in each tree, connected to all of the edges in that tree. For illustration, see, for example, Figure 2 in Aas et al. (209). The C-vine was shown to be good choice for neuronal datasets (Onken & Panzeri 2016). as they often include some proxy of neuronal population activity as outstanding variable, strongly correlated with the rest. This variable provides a natural choice for the first conditioning variable in the lowest tree. In the neuronal datasets from Henschke et al. (2020). this outstanding variable is the global fluorescence signal in the imaged field of view (global neuropil).
To construct a C-vine for describing the neuronal and behavioral data from Henschke el al. (2020). we used a heuristic element ordering based on the sum of absolute values of Kendall's 's of given element with all of the other elements. It was shown by Czado al. (2012) that this ordering facilitates C-vine modeling. For all of the animals and most of the recordings (14 out 16), including the one used in Figure 3, the first variable after such ordering was the global neuropil activity. This again confirms, that a C-vine with the global neuropil activity as first variable is an appropriate model for the dependencies in neuronal datasets.
A.6 ALGORITHMIC COMPLEXITY
In this section, we discuss the algorithmic complexity of the parameter inference for a C-vine copula model.
The parameter inference for each of the bivariate Copula-GP models scales as Onn), where n is the number of samples, since use a scalable kemnel interpolation KISS-GP (Wilson and Nickisch 2015). As we mentioned in Sec. A.S full m-diimensional C-vine model requires m(m - 1)/2 bivariate copulas, trained sequentially.
15

---

Under review as a conference paper at ICLR 2021
As result, the (n) GP parameter inference has to be repeated m(m 1)/2 times, which yields O(n m²² complexity.
In practice, the computational cost (in terms of time) of the parameter inference for cach bivariate model varies from tens of seconds to tens of minutes. The heuristic model selection is designed in such way, that discards independent variables in just around 20 seconds (line 3 in Alg. 2). As a result, most of the models are quickly skipped and further considered can be neglected. When the model is evaluated, the Independence components are also efficiently "skipped" during sampling, as ppcf function not called for them. The Independence models also add zero to C-vine log probability, so they are also 'skipped" during log probability calculation. They also reduce the total memory storage, no GP parameters, which predominate the memory requirements, are stored for these models.
Independence models, and their contribution the total computational cost
In conditional C-vine trained on real neuronal dataset with 109 variables, 5253 out of 5886 (89%) bivariate models were Independence, which leaves only 633 non- non-Independence models.
practice, this means that the algorithmic complexity the model is much better than naive theoretical prediction (n - m²), hased on the structure of the graphical model. Suppose that the actual number of the doneendeeeeeece models N,I in vine model is much smaller than m(m 1)/2 and can be characterized
by an effective number of dimensions m ~ /NnI. In this case, instead of the ((m)) scaling with the number of variables, the complexity highly depends on the sparsity of the dependencies in the graphical model and scales with as O(n Na1) ~ O(n mers).
Therefore, the method is especially efficient the datasets with low effective dimensionality m.ff, such as neuronal data. The number of variables m itself has little effect on the computational cost and memory storage.
B MORE VALIDATION ON SYNTHETIC DATA
Computing infrastructure We developed our framework and ran the majority of our experiments (de- scribed both in the paper and Supplemental Material) on an Ubuntu 18.04 LTS machine with 2 Intel(R) Xeon(R) Gold 6142 CPU @ 2.60GHzz and Ix GeForce RTX 2080 GeForce RTX 2080 Ti GPUs. For training C-vine models, we used another Scientific Linux 7.6 machine with Intel(R) Xeon(R) Silver 4114 CPU @ 2.20GHz and 8 x GeForce RTX 2080 Ti GPUs
Code availability
Code will be made available on GitHlub upon paper acceptance.
B.1
MODEL SELECTION FOR BIVARIATE COPULAS
Synthetic data
We generate artificial data hy sampling from a copula mixture, parametrized in two different
ways:
mixing concentrations of all copulas were constant and equal 1/N (N = number copulas), but copula parameters were parametrized by the phase-shifted sinus functions:
0, A; sin
N + 2TT + B,,
x E [0,1]
(A4)
where is the index of the copula a mixture, m For Clayton and Gumbel copulas, the absolute value of the sinus was used. The amplitudes A, were chosen to cover most of the range of parameters, except for extremely low or high Øs for which all copula families become indistinguishable (from independence or deterministic dependence, respectively).
2. copula parameters were constant, but mixing concentrations á were parametrized by the phase- shifted sinus functions (same as Eq. A4 with A; = B, = 1/N and m = 2). Such parametrization ensures that the sum of mixing concentrations remains equal to one (Σ'' á = 1). Yet, each é turns to zero somewhere along this trajectory, allowing us to discriminate the models and infer the correct mixture.
We tested the ability of the model selection algorithms to select the correct mixture of copula models, the same the from which the data was generated. We generated 5000 samples with equally spaced unique inputs on [0,1].
Identifiability tests
Both model selection algorithms were able to correctly select all of the 1-component and most of the 2- component models on simulated data. For simulated data with larger numbers of components (or 2 very similar components), the WATC of the selected model was either lower (which is possible given limited number of samples) or close to the WAIC of the correct parametric model. In other words, the difference hetween the
16

---

Under review as a conference paper at ICLR 2021
WAIC of the correct model and of the best selected model never exceeded the WAIC no! 0.05, which we set up a criteria for passing the test: AWAIC < WAICtest tot. Since all the tests were passed success- fully, we conclude that both algorithms are capable of finding optimal or close-to-optimal solutions for copula mixtures.
more detailed report on the model identifiability tests Tables A3 A7 below illustrate the search for the best model. The copula model names in these tables are shortened to the first two letters, e.g. Gumbel becomes 'Gu', Frank becomes 'Fr'. The information in these Tables provides some intuition on the model selection process and the range of WATCs for the correct incorrect models. The final selected models are shown in bold.
Table A3 demonstrates that both greedy and heuristic algorithms can identify the correct single copula model. Some key intermediate models (M in Alg. 22) with their WATCs are listed in the table, along with the total duration of simulations (T, in minutes) RTX 2080Ti for both algorithms.
A4ble|A4hhove the identification of the mixtures with 2 components, where the copula parameters were con- stant (independent of æ) and mixing concentrations & were parameterized by the phase-shifted simus functions (Eq. A4). All of these models were correctly identified with both algorithms. The mixtures with 2 components, where the copula parameters varied harmonically (as in Eq. A4) but the mixing concentrations & were con- stant, were harder identify. Table A5 shows that few times, each of the algorithms selected model that was better than the true model (WAlCbest WAICerue 0). The greedy algorithm made one mistake, yet the model it selected was very close to optimal. Such misidentification happens due to the limited number of samples in a given synthetic dataset.
show the model selection for 3 component models. Again, in Tables A4 A5 either @ or & was constant. Here, the model selection algorithms could rarely identity the correct model (due to overcompleteness the mixture models), but always selected the one that was very close optimal: WAIC WAIC < WATCtest.to-
Tables A6 A7
Note, that WAICrest.tof is different from waic.tol. We have set waic.tol I for comparison against Inde- pendent model to such small value (10x smaller than WAICresrod) because we want to avoid making false assumptions about conditional independences in the model. Also note, that the WAIC of the true model de- pends on the particular synthetic dataset generated each test. Therefore, the final WAIC the left and in the right columns Tables A3
A3 A7 can slightly different (yet, right within WATCtest.tol) tot).
B.2 ACCURACY OF ENTROPY ESTIMATION
Gaussian entropy
Estimation error
Gaussian entropy
Estimatiio error
100
A
altty) - Hose Her, bit
B
C
Mogs bit
True
Copula-GP
10
Mly), hit
40
6
May
BI-KSG
4
20
2-
0
0.0
0.5
1.0
0.0
0.5
1.0
10
20
10
20
corr. coef.
P
con.. coef.
dim y
dim y
Figure A4: Accuracy of the entropy estimation for multivariate Gaussian distributions. A Entropy of the 20-dimensional multivariate Gaussian distribution for different correlation coefficients p. B Es- timation error for the entropy shown in A. C Entropy of the multivariate Gaussian distributions with p 0.99 and varying dimensionality. D Estimation error for the entropy shown in C.
In this section, we consider fixed copula mixture model with known parameters 0 and test the reliability of the entropy estimation with Monte Carlo (MC) integration. We test the accuracy of the entropy estimation on a multivariate Gaussian distribution, with cov(gk, M;) = + (1 p) d.j. where d.j is Kronecker's delta and E [0, 0.999). Given known Gaussian copula, we estimate the entropy with MC integration and compare it the analytically calculated true value. We up tolerance AH 0.01 (dim y). As result, for every correlation P (Fig. A44 and every number dimensions dimy (Fig. A4C-D)) the Copula-Gi GP provides an accurate result, within the error margin. In Figure AA44 BI- KSG estimates (Gao et al., (2016) obtained on the dataset with 10k samples are shown for comparison. This experiment 1) validates the MC integration; 2) validates the numerical stability of the prohability density function of the Gaussian copula up to specified maximal p 0.999 (for p > 0.999 the model indistinguishable from the deterministic dependence in ug).
17

---

Under review as a conference paper at ICLR 2021
B.3 VALIDATION ON NON-GAUSSIAN COPULA
Copula-GP guaranteed to produce unbiased entropy estimates when the true dependency matches the para- metric model (i.e. is Gaussian, Clayton, Frank or Gumbel, or any linear mixture of those), assuming that the implementation of log likelihood and sampling is correct, which is ensured by the tests and correctly estimated parameters.
Figure Z tested the entropy and mutual information estimates the Gaussian copula. Here, we follow the same procedure and generate samples from a C-Vine, in which each bivariate element a Clayton copula. varied the parameter / of the Clayton copula linearly from to 2.
The observed result is qualitatively identical to the validation on Gaussian copula model (see Fig. A Copula- GP measures both entropy and mutual information exactly (within integration tolerance, see Fig. AS while the performance of the non-parametric methods on this dataset is lower.
Clayton with corr-G.
10-
Thue
1.50
CHAM iniegnctal
BHAS
KOU
MINE102
Ba (3lm)
#votuako 98m
Figure A5: Conditional entropy Hy]x)) (left panel) and mutual information [(x, y) (right panel) measured by different methods on synthetic data, generated from Clayton copula. The performance of various estimators is qualitatively similar to Fig. A with the Gaussian copula.
B.4 PERFORMANCE ON UCI BENCIIMRR DATASET
Here we compare our method with the similar method from Lopez- et 2013 Hernández-Lobald ett al. (2013). Apart from differences in implementation (SVI vs. EP, (un)/conditional marginals), our Copula- GP method has two conceptual improvements 1) we use non parametric conditional marginals, instead of unconditional: 2) we linear mixtures of copula elements with qualitatively different tail dependencies. Using the UCT shuttle dataset (Frank et al. 2) Copula-GP model which was either allowed to use only one copula element, or a mixture of any number of elements (selected with the model selection method, see Appx.[A.4] We used vine copulas truncated at 3 or 5 trees, ran them on the training set and evaluated the log-likelihood on the test set (higher is better).
compare the performance of the Lopez-Paz. model vs.
As a result, Copula-GP restricted to the use of only one copula element produced similar results to Lopez- Paz 2013), while our mixture model considerably outperformed both of these single element models.
Table A2: Average test log-likelihood for Copula- GP vs. Lopez-Paz et al. 2013
Trees
Lopez-Paz et
Single copula
Mixture
3
3.645:10.227 4.755+:0.89
3.759 4.470=10.215
4.326::0.170 5.115:10.10
9
C INTERPRETATION OF COPULA-GP MODELS
C.I CONDITIONAL VS. UNCONDITIONAL MARGINALS
In this section, we explain the difference between copula models with conditional vs. unconditional marginals.
Consider pair of neurons, which both fire strongly in a certain context (e.g. in the reward zone of the VR environment, Fig. A6A) We simulate the activity of these neurons using Poisson model (Fig. (Fig.A66)) and then convolve the simulated spike count data with the exponential kernel, order to generate continuous data similar to the recordings shown in Fig. The resulting generated activity is shown in Fig. AAC.
Let us then use unconditional marginals to project simulated neuronal recordings onto unit cube (Fig. A6 D.
left). The color of the data point shows the position X, and one can see that the orange and red points (X 0.5) concentrated in the upper right corner of the plot. If we plot histogram for one of the neurons and only
18

---

Under review as a conference paper at ICLR 2021
B
C
2.0
10
ta
20
1.0
0.5
0.0
0.2
Unconditional marginals
E
Copula parameters for 1 vs
Mixing param for 1 vs O
1.0
pormalired
2.0
Clayton
Clayton
1.00
DOC
Dlyyan 0"
1.0
Cayton 180"
0.0
0.5
0.25
0.0
0.5
1.0
Y: for X=0.55
0.0
0.2
0.4
0.0
1.0
G
Conditional marginals
F
Positio in VR. [cm]
Position VR, [cm]
0.0.0.6
0.6-0.7
0.7-0.9
0.9-1.0
1.0
1.0
X
0.5
0.5
0.0
0.0
0.0 0.0
0.0
0.5
1.0
for X>0.5
1.0
Figure A6: Toy model of two independent neurons with the same conditional firing rate: A average firing rate depending on X; raster plot with spike counts for one of the neurons; simulated calcium traces from 10 trials for one of the neurons; D data-points projected onto unit cube by the probability integral transform based on unconditional marginals; colored points show joint neural data-point (Y, recorded at X, where X is color-coded (as in A); the histogram shows conditional marginal distribution p > 0.5); E mixture model fit for the data in D (with unconditional marginals); density plots comparing predictions of the mixture model (shades) with the actual data (black outlines); the fit is poor due to non-uniformity of empirical marginal distributions; same as D, but based on conditional marginals, estimated using fastKDE.
the right
orange-red points, which corresponds p(Y X 0.5). we see that far from being uniform (Fig. A6D, right).
In previous works on conditional copulas (Heránndez -oriaddezz-obbttoo et al. 2013). the marginals were assumed to be unconditional. Such statistical copula-based models would ignore the non-uniformity of the empirical copula marginals, and still attempt describe the data with some heavy-tailed copula (Fig. A6E). However, uniformity of copula marginals one of the main model assumptions, and it does hold in this case. As result, the model poor (Fig. F)
Modelling with conditional marginals, as our Copula-GP model, is radically different. The non- parametric conditional marginal model (estimated with fastKDE) reflects the fact that both neurons are typically firing strongly around X 2s 0.7, and would transform the marginals accordingly (Fig. A6G. left). As result, all the model assumptions, including uniformity copula marginals, will strictly hold (Fig. A6G, right). Addi- tionally, our copula model would focus solely on the joint trial-to-trial to-taall variability (i. on noise correlations) of these neurons, regardless of the absolute values of their firing rates. In this example, the model finds no noise correlations (best fit for Fig. A6G is Independence). which is the correct result for the simulation of two independent neurons.
Therefore, conditional marginals make our model more flexible and capable holding the modelling assump- tions even on complex neuronal data with highly variable marginal statistics. Also, for the neuronal data in particular, the model becomes more interpretable, the conditional copula c(u*(z) focuses only the joint trial-to-trial variability of neural activity. Therefore, the statisties corresponding to the network becomes tirely separated from the single unit statistics. This makes our Copula-GP model particularly well suited for neuroscience applications.
C.2 MODEL PARAMETERS FOR THE BIVARIATE NEURONAL AND BEHAVIORAL EXAMPLES
In this section, provide visualisations for the parameters of the bivariate copula models from Figure 3C-F and discuss the interpretability of these models.
19

---

Under review as a conference paper at ICLR 2021
A
Inter-neuronal noise correlations
Copula parameters e
Mixing parameters
0-60
60-120
120-140
140-160
1.0
Navrnne
0.5
E 0.s
Clayton
Pearson's mho
-0.5
0.0
50
100 (cm)
150
100 Position VR, [em)
50
150
Neuron: 3
Poxition
Behavioral modulations
Copula parameters
Mixing parameters
0-50
60-120
120-140
140-160
0.25
E.poemstieed
Neuronn 60
Clayton 0
0.00
Gumbel 270
Pearsonss the
-0.75
0.0
300
150
Ldks
Fosttion [cm]
Posttion [cm]
Figure A7: Parameters of the copula mixture models. From left to right: copula probability densities (same as Fig.3C-D); a list of selected copula elements; copula parameters 0: mixing concentrations s. These plots are provided for: A the noise correlation example; B the behavioral modulation example.
Figure A7 shows the probability density of the joint distribution of two variables and the parameters of a corresponding Copula-GP GP mixture model. The plots on the left repeat Fig.3C-D and represent the true density (outlines) and the copula model density (blue shades) for each part of the task.
In the noise correlation example (Fig. (A7A), we observe the tail dependencies between the variables (i.c. concentration of the probability density in a corner of the unit square) around [0-60] cm and [140-160] cm of the virtual corridor. There is only one element with tail dependency in this mixture: Clayton 90 copula. On the right-most plot in Fig. AA74, see the mixing concentration for the elements of the mixture model. The concentration of Clayton 90" copula (orange line) is close to 100% around 20 cm and 150 cm, which agrees with our observations from the density plots.
The confidence intervals (=2r) for the parameters approximated with Gaussian processes are shown with shaded arcas in parameter plots. These intervals provide a measure of uncertainty in model parameters. For instance, when the concentration of the Gaussian copula in the mixture is close to 0% (r around 20 and 1.50 cm), the confidence intervals for the Gaussian copula parameter (0, blue shade) in Fig. A77 become very wide (from almost 0 to 1). Since this copula element is not affecting the mixture for those values of æ, its parameter has no effect on the mixture model log probability. Therefore, this parameter is not constrained to any certain value. In similar manner, we see that the variables are almost independent between 60 and 120 cm (see density plots the left in Fig. A7). Both copula elements can describe this independence. As result, the mixing concentrations for both elements have high uncertainty in that interval of a. Yet, Gaussian copula with a slightly positive correlation is still a bit more likely to describe the data in that interval.
The copula parameter plot in Fig. A7 A also shows Pearson's p. which does not change much in this example and remains close to zero. This illustrates, that the traditional linear noise correlation analysis would ignore (or downplay) this pair of neurons as the ones with dependence. This happens because the Pearson's p only captures the lincar correlation and igmores the tail dependencies, whereas our model provides a more detailed description of the joint bivariate distribution.
In the behavioral modulation example (Fig. 37 we observe more complicated tail dependencies in the den- sity plots. The best selected model supports this observation and provides a mixture model with 3 components, 2 of which have various tail dependencies. The Clayton 0 copula (orange) describes the lower tail dependence observed in the second part of the virtual corridor with gratings (around [60-120] em, see Fig. 3A for task structure). This dependence can be verbally interpreted as follows: when there is no licking, the Neuron 60 certainly silent; but when the animal is licking, the activity of Neuron 60 is slightly posirively correlated with the licking rate.
These examples illustrate, that by analysing the copula parameters and the mixing concentrations of the Copula- GP mixture model, one can interpret the changes in the bivariate dependence structure. Just like traditional tuning curves characterize the response of a single neuron, our mixture model characterizes the 'tuning' of the dependence structure between pairs of variables to given stimulus or context. Knowing the qualitative properties of the copula elements that constitute copula mixture, one can focus on the dominant element of the copula mixture for every given conditioning variable and describe the shape of the dependence.
20

---

Under review as a conference paper at ICLR 2021
C.3 ABLATION STUDY
After the model selected, one can ablate elements one by one and check the WAIC of these models. Following this procedure, one can test whether all of the elements are important in mixture.
We performed the ablation on the models from Fig. A7 originally selected model. Ablation of the Clayton copula in Fig A7B increased WAIC the most (from -0.043 -0.033), while both elements in Fig A7 were equally important (WAIC after ablation was the same (within tolerance)).
All models after ablation had higher WAIC than the
C.4 DETECTION OF HEAVY TAILS
Since Copula- GP framework based on a linear mixture model, is possible to calculate the probability that a certain data aata-ppoit was generated from heavy-tailed component of the distribution (c.g. Clayton). Given a mixture model with estimated parameters, described by:
K
(ula)) x)) Σ (,()))((;;(((()))
where M, CDF, (w,l)), is conditioning variable and y is vector of neuronal and/or behavioral record- ings. Then, for a data-point yz, this probability is equal to:
p(clayton (z, p(y]clayton, 2) p(clayton/x)
Ejeclayton (1; (x))
p(y|2))
ΣeET By (ux
By thresholding p(clayton|z, y) and p[y[z). select those data pitaappoints that constitute the heavy tail of the distribution, which is described by Clayton copulas in our copula mixture model.
APPENDIX REFERENCES
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational gaussian process classi- fication. 2015.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS- GP). In International Conference on Machine Learning, pages 1775-1784, 2015.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2.
MIT press Cambridge, MA, 2006.
Tim Bedford and Roger M Cooke. Probability density decomposition for conditionally dependent random variables modeled by vines. Annals of Mathematics and Artificial intelligence, 32(1-4):245-266 2001.
Tim Bedford and Roger M Cooke. Vines: A new graphical model for dependent random variables. Annals of Statistics, pages 1031-1064 2002.
Flif Acar, Christian Genest, and Johanna NešLehová. Beyond simplified pair-copula constructions. Journal of Multivariate Analysis, 110:74-90 00,799,,,20 2012.
Ingrid Hobæk Haff, Kjersti Aas, and Amoldo Frigessi. On the simplified copulc construction- simply useful or too simplistic? Journal of Multivariate Analysis, 1296))299-11331 2010.
David Lopez-Paz, Jose Miguel Hernández-Looaal Lobato, and Ghahramani Zoubin. Gaussian process vine copulas for multivariate dependence. In International Conference on Machine Learning, pages 10-18, 2013.
José Miguel Hernández-Lobate James R Lloyd. and Daniel HerrándezzLoball ditional copulas with applications to financial time series. In C. J. C. Burges, L. Bottou. M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors. Advances in Neural Information Processing Sys- tems pages 7736-1744, Curran Associates, Inc., 2013. URL http:/ /ptp::///appaaaaaae nips wttppper/ccppp/////p// 5084- 1844-gussiian pdf
Gaussian process con-
Andrew Frank, Arthur Asuncion, et al. UCI machine learning repository, 2010. URL http://archive. ics. uci. edu/ml, 15:22, 201
21

---

Under review as a conference paper at ICLR 2021
Table A3: The model selection histories for 1-element mixtures
True Model
Greedy
Heuristic
Search attempts
WAIC
Search attempts
WAIC
T
Ga
Ga
-0.1619
25m
Ga
-0.1513
3 m
GaFr
-0.1610
InGaGu Gu Gu
-0.1499
-0.1619
C C1
-0.1498
Ga
-0,1513
Fr
Fr
-0.1389
57 m
Ga
-0.1400
3 m
FrC1
-0.195
InGaGu lu Gu Gu
-0,1391
Gu
-0.1396
c C1
-0.1391
-0.1396
-0.1509
Fr
-0.1389
C1
CI
-0,5225
37
Ga
-0.3825
5
Gu
-0.5226
InGaGu Gu
-0.4943
Gu
-0.5225
-0.5303
CT
-0.5224
C1
-0.5311
Gu
Gu
-0.6267
43
Ga
-0,5555
7m
Gu Gu
-0.6268
InGaGu Gu
-0.5988
Gu G
-0.6267
InGaCI (1
-0.5946
Gu
-0.6230
GaGu Gu
-0.6040
-0,6050
Ci
CI
-0.5389
m
-0.3922
5
CI70
-0.5389
Gu C1
-0.5047
CT
-0.5389
-0.5409
C1
+0,5410
Go
Gu
-0.6137
55
Ga
-0.5501
7
Gu Gu 270
-0.6144
InGaGu Gu 270 Gu Gu
-0,5893
Gu
-0.6145
169
-0.5831
-0.6144
GaGu Gu
-0.5887
Gu
-0.6137
+0,5950
C1
1
-0.5566
3
Ga
-0,3932
7 m
11
00.5822
InGaGu Gu
-0.4956
c1
-0.5882
-0.5493
180
-0.5565
C1
-0.5489
Gu :30
Gu 180
-0.6131
43 m
Ga
-0.5533
6 6m
Gu 1
-0.6164
InGaGu
-0.6091
-0.6163
18
-0.6045
Gu 1su
-0.6131
Gu
-0.6154
2
C7
-0.5434
23 m
+0.3909
5
Gu
-0.5433
InGaGu Gu
-0.5094
C7
-0.5434
-0.5535
C/270
-0.5548
270
G7
-0,5928
51 m
Ga
-0,5763
-0.5934
InGaGu 2700
-0.6277
-0.5935
InGaCI" 1180
-0.6179
"
-0.5931
Gu
-0.6300
Gu
-0.5928
22

---

Under review as a conference paper at ICLR 2021
Table A4: The model selection histories for 2-element mixtures with constant 0 and variable ø
True
Greedy
Heuristic
Model
Search attempts
WAIC
T
Search attempts
WAIC
T
Gu
Ga
-0.1877
101
Ga
-0.1922
11
Ga
GaGu
-0.2855
Gu Gu
-0,3070
GaGu
-0.2855
InGaCI C127
-0.2996
GaGu 00 R
-0.2856
Gu
-0.3082
GaGu 90 CI 270 2r
-0.2856
G7 Gu
-0.3076
GaGu 90 c1 FrGu
-0.2856
GaGu
-0.3391
C1
Gu 20 Ga
-0.2854
Ga
Fr
-0.1635
87
Ga
-0.1600
5
C1
FrC1270
-0.2707
GuG
-0.2687
FrCI
-0.2747
C GaC1270
-0,2835
FrCI#" 180
-0,2782
-0.2845
FrCl c7
-0.2781
GaC1270
-0.2821
Gu 180
Gu
-0.1681
99
Ga
-0.1534
8 8
Fr
Gu
-0,2099
InGaGu 2u7 oo
-0.1993
-0.2101
C
-0.1977
Gu C
-0.2105
InGaGu 180
-0.2074
Gu " In
-0.2106
FrGu 350
-0.2104
-0.2099
Gu
FrGu 180
-0.2099
CTU
-0.1587
92
Ga
-0.1652
5
C
FrC1
-0.2600
InGaGu Gu Gu
-0.3142
20 FrCI' 270
-0.3173
-0.3430
-0.3176
c
-0.3448
C1 2
-0.3176
FrCI 2
-0.3175
c1"'r"
-0.3190
-0.2204
103
-0.1965
m
Gu 270
18C
-0.3488
InGaGu Gu Gu
-0.3591
-0.3874
-0.3688
Gu
-0,3877
GaGu
-0.3771
FrCI 180 Gu 2
-0,3878
Gu 270 21 180
-0.3772
FrCI [CC]180Ga 2770 9
-0.3878
Gu 90
Gu 270 250
-0.3888
Table A5; The model selection histories for 2-element mixtures with constant & and variable
True Model
Greedy
Heuristic
Search attempts
WAIC
Search attempts
WAIC
Gu
Gu
0..1419
60
Ga
-0.1538
10 m
Ga
-0,2022
InGaGu 180 270 Gu Gu
-0.2320
Gu
-0.2024
InGaCI 150
-0.2218
-0.2024
Gu"
-0.2321
FrGu 90
-0.2021
GaGu 90
-0.2326
WAICbest WAIC
-0.0013
Ga t
Go
-0.1495
56m
Ga
-0.1062
7 m
"r FTC
-0.1894
InGaGu "T
-0.1747
-0.1915
-0.1783
Gu FrC1 or
-0.1902
GaGu "C)
-0.1812
-0.1915
2aC
-0.1801
WAIC WAIC
0,0032
23

---

Under review as a conference paper at ICLR 2021
True Model
Greedy
Heuristic
Search attempts
WAIC
Search attempts
WAIC
Gu ts
Gu
0.1600
58
Ga
-0.1331
8 8m
Fr
Fr
-0.2191
InGaGu Gu 270 Gu Gu
-0.1944
Gu
-0.2195
CI C1
-0.1936
FrCI"""C"""""
-0.2190
GaGu GaGu C1
-0.1945
FrGu 180
-0.2190
GaGu 180
-0,1992
WAICbest WAIC
-0.0094
Cl C
2u
-0.0253
62 m
Ga
-0.0079
5 5
Gu C1 Gu 180 C1
-0.2383
InGaGu Gu Gu Gu
-0.1904
-0.2506
InGaCI 1800
-0.2330
Gu 180 9 C In
-0.2509
Cr
-0.2361
Gu C1 C
-0.2508
C1
-0.2586
27
-0.0242
69
Ga
-0.0083
27 270
Gu C
-0.2499
InGaGu 180 Gu Gu
-0,2277
f Gu
-0.2517
InGaCI
-0.2535
Gu CI Gu In
-0.2518
GaCI C1
-0.2549
Gu 1
180 Gu 2 InC1"
-0.2518
G CI 180
Gu InCI
-0.2518
CI C1180 Gu 270
-0.2500
WATCbest WAIC -0.0098
Table A6: The model selection histories for 3-element mixtures with constant & and variable ø
True Model
Greedy
Heuristic
Search attempts
WAIC
Search attempts
WAIC
Ga
Gu"
-0.1399
44 m
Ga
-0.1252
6m
C1
Gu Cu
-0.2494
InGaGu 180 Gu 27 G
-0.2481
Gu
"
-0,2519
I "1900
-0,2565
Gu Gu'
-0.2518
GaCI 11
-0.2564
Gu
-0.2494
WAICbest WAIC
-0.0036
WAICsest WAIC
0.0014
Fr
Fr
-0.0591
Ga
-0.0489
Ci
FrC190
-0.1460
InGaGu 180 2T
-0.1573
Gu
GrC
-0.1730
InGaCI
-0.1578
FrC 21100
-0.1736
Ga 11
-0.1621
FrC G CI
-0.1734
Gu C
-0.1731
WAIC
WAIC
0.0059
Fr
Fr
-0.0741
87 m
Ga
-0,0618
9
180
180
-0.1513
InGaGu 180 270 Gu Gu
-0.1567
27o
160
-0.1707
9 2
-0.1670
-0.1708
InGaGu 270 180
-0.1680
FrC1180
G 180
-0.1711
InG 180
-0.1695
FrCI 2u 9 Gm
-0.1710
InGu 110
-0.1735
Gu 270 F8
-0.1703
WAIC
WAIC Ctrue:
-0.0011
Gu
Gu
-0.1695
47 m
Ga
-0.1477
11 m
Gu
20
-0.3040
InGaGu 180 270 Gu
-0.2986
Gu Gu
-0.3234
-0.3033
Gu 180
-0.3233
GaGu C1 1180
-0.3054
Gu 180 91 Gu
-0.3234
GaGu 180 9n"
-0.3111
Gu "
-0.3113
24

---

Under review as a conference paper at ICLR 2021
Table A7: The model selection histories for 3-element mixtures with constant & and variable
True
Greedy
Heuristic
Model
Search attempts
WAIC
Search attempts
WAIC
Ga
Fr
-0.0177
66 m
Ga
-0.0142
13 m
Cl
FrGu
-0.1284
MiGaGut*"Ge 2 Gu
-0.1291
Gu
Gu
-0.1407
InGaCI
C
-0.1289
Gu
-0,1423
InCI Gu
-0.1317
Gu C
-0.1435
InCI"CI"
-0.1346
FrGu Gu G
-0.1432
-0.1301
Gu
-0.1451
GaCI C1
-0.1313
WAICbest WAICrruet
0.0132
WAICbest WAIC
0.0068
Fr
Fr
-0.0265
71m
Ga
-0,0192
9
C1
FrGu270
-0.1290
Gu
-0.1411
Gu
Gu
-0.1445
In | C1
-0,1429
FrGu2""Gu' C1
-0. 1451
InGaGu 180 C1
-0.1474
FrGu270 1800 C1
-0.1466
InGaGu" 180 C1
-0.1472
FuGu C In
-0.1468
InCI Gu
-0.1477
Gu
-0.1451
WATCbest WATCtrue'
0.0109
WATCbest WATC
0.0010
Fr
Fr
-0.0129
61
Ga
-0,0185
6
270
-0.1105
Gu
-0.1309
Gu
FrGu Gu
-0.1237
CI
-0.1326
Fr 280
-0.1254
InGu
-0.1393
FrGu 270 Gu CI Gu 180
1248
nnG
-0,1334
Gu Gu Gu 270 Fr
-0.1234
InGu Gu
-0.1326
WATCbest WAICtrue
0.0094
WATCbest WAIC
0.0088
Gu
Gu
-0.0756
55
Ga
-0,0454
7
Gu 180
Gu 90
-0.2380
InGaGu G8 27 Gu
-0,2476
C
Cr
-0.2556
InGaC1 91 1 c1
-0.2459
Gu"CI""CIGGa
-0.2591
Gu
-0.2493
Gu CI
-0,2590
-0,2559
Gu
-0.2555
CI
-0.2538
WAlCbest WAIC truei
0.0026
WAICbest WAICtruei
0.0006
25


Reviews:
[{"review_id": "wwjJdozAGUQ", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors exploit the expressive power of Copula mixtures to model time-varying multi-modal data, and employ Gaussian Processes to model the time-varying copula parameters. They demonstrate the efficacy of their method using information theoretic metrics on a synthetic dataset and a real-world joint neural-behavioral dataset from a neuroscience experiment.  Results demonstrate that the proposed techniques are comparable to the state of the art nonparametric methods, while being more scalable due to the use of stochastic optimization based methods that are commonly used with parametric methods. \n\nThe paper was quite thorough in the motivation, development and empirical analysis of the proposed technique. The use of Copulas, that are commonly employed in the Finance community, but are rare in statistical neuroscience, should interest the more theoretically inclined reader. Given the accelerating trend of collecting long-term joint neural and behavior in experimental neuroscience, the authors make an interesting and timely contribution to the statistical neuroscience literature.\n\nI found that the motivations of the paper were difficult to extract from the Introduction, as there was substantial use of jargon and the intuition behind the technical results were not accessible. To encourage the adoption of their methods in the neuroscience community, the authors should consider improving the readability of their manuscript by making it more friendly to the non-statistician reader. \n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "j0dDF_7j-Tk", "reviewer": "AnonReviewer4", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The authors develop a Gaussian process vine copula model, very much in the flavor of the modeling approach of Lopez-Paz et al. (2013). The improvement to the earlier work seems to be a framework for flexible copula modeling including a copula mixture model, approximate inference, model selection, and calculation of mutual information. The paper on its own is a modest improvement on existing work and is both an engineering accomplishment and has the potential for a useful model. The paper is exceptionally well-written and clear. I found it a breeze to read and I credit the authors for that. However, both the validation and the motivation for the model (namely characterizing the probabilistic relationships for neural and behavioral variables) seems particularly thin and could be substantially improved. The paper falls short on these points to such a degree that I am hesitant to recommend acceptance since it notably impacts my evaluation of the significance of the works\n\n## Major points:\n\nFirst, the authors claim (page 7, paragraph 2) that their model out-performs all non-parametric methods. This is by no means obvious form Fig2b,c. Moreover, If the model performs \"similarly\" to MINE, then it is maybe not worth using when a convenient technique already exists. What is the advantage of the present model?\n\nSecond, the authors show that the estimation of mutual information is only unbiased for a narrow range  of distributions (Gaussian or transformed Gaussian for small dimensions according to Fig 2) and fails for the heavy-tailed Student's-T. However, many neural and behavioral variables are themselves heavy-tailed and the authors did not demonstrate that the real data are sub-gaussian.\n\nThird, it seems like there was a wasted opportunity with this paper. The authors spent $\\approx$ 2/3 of a page discussing the estimation of mutual information without motivating why that was a good example metric that could be derived from their model. This modeling framework is an opportunity to determine virtually any expectation over the entire distribution and it is entirely possible that MI is neither all that interesting, nor does it play to the strengths of the model. They then describe changes  to the pairwise distributions of variables from the copula model but we didn't need the copula model to estimate.\n\nFinally, what are we to make of the real data results *vis-a-vis* the validation experiments with simulated data in Fig 2? Besides what I mentioned above regarding the tails of the real data, its not clear that the variance explained wouldn't behave differently. Could the authors report the variance explained for their simulation experiments?\n\n## Minor points: \nThe authors state (page 7, last paragraph) That the \"stimulus-related changes in the joint variability of the two neuronal signals are commonly described as _noise correlations_.\" But, isn't that the definition of _signal correlations_?\n\nthe authors state (page 8, paragraph 5) \n>The Copula-GP “estimated” (dashed line) almost perfectly matches the “integrated” result, which suggests that the model was able to tightly approximate both $p(u^x|x)$ and $p(u^x)$, and, as a result, $I(x, \\{u|^x_{i<N} \\})$.\n\nHowever, it is not clear at all that the later statement regarding $I(x, \\{u|^x_{i<N} \\})$ follows from the former. In fact Figure 2 demonstrates that the integrated result may not estimate $I(x, \\{u|^x_{i<N} \\})$ well at all.", "extra_scores": {}, "extra_text": {}}, {"review_id": "WCwtFmssM4z", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This manuscript models the conditional joint distribution over variables by using Copula models and copula vines.  The experimental data shows that when the observed variables are highly correlated that the proposed approach improves estimation of entropy over competing benchmark approaches (MINE and KSG) when the variables are highly correlated. Synthetic results demonstrate good improvements, and application to real scientific data seems promising.\n\nStrengths:\nApplication to real, novel scientific data shows potential utility of the model.\nSynthetic results show a good improvement over competing methods, albeit in a limited setup (highly correlated variables)\nMixtures of copulas seems an effective way to produce model complexity, and by linking it to a GP can make sure that it's smooth over the conditioned variable x.\n\nWeaknesses:\nNot all modeling steps are clear.  In particular, the interaction between the GP prior and the model selection step is undescribed.  Since there are multiple ways to do this, needs a full description.\nNo comparison to more neuroscience focuses techniques.\nExperimental setup and utility is not fully described.\nLacks ablation studies to elicit key model components.\n\nQuestions:\nThe choice to make the link function on the different copula families all dependent on the same $f$ seems like a strange and limiting choice.  Why was this choice made?\n\nThe choice to model calcium trace level and not neural spiking seems mathematically convenient for this method, but it's not clear to me that this is the correct scientific choice.  Typically, calcium imaging traces are preprocesses to extract spikes.  Is there a scientific rationale for using the raw data?  Or was this primarily motivated by avoiding discrete measurements that would be harder to model in the copula?\n\nPlease describe more of the scientific setup.  For example, why is there any relationship to licks outside of the reward?\n\nWhat is the scientific question on this neuroscience application?  Why is entropy relevant, rather than using one of the many predictive problems?\n\nUpdate after author response:\nMost of my methodological concerns have been address (except for the ablation studies).  The scientific application here is not super well-motivated.  It would improve the article greatly to show a greater utility (or at least, clearly describing future utility for answering scientific questions).  It is mentioned that \"A full application of the method to study the dependence of contextual signals in mouse visual cortex will be the focus of a follow-up publication.\"  That's vague; it would be nice to at least clearly discuss how this could be used to facilitate or enhance these scientific experiments.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "ItIEBJmeoia", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "- This is an interesting neuroscience application where Copula estimation has shown to be effective. \n\n- Having said that, I believe that the technical novelty is minimal. To the authors' credit, they did not claim a huge theoretical edge. They honestly reported the findings of the paper. \n\n- The paper is well written. Ideas flow very neatly. \n\n- p2: \"It was previously shown that such a combination of parametric copula models with GP priors ...\": Precisely.. Therefore, given the work by Hernandez-Lobato et al. (2013), as well as the works which have eventually built on top of it, novelty of the proposed method is minimal. I am not saying the application is not useful though.  \n\n- Along the same lines from above, equations (2), (3) and (4) which depict the core of the method, are all taken from seminal previous works of copula mixture estimation. \n\n- p3: \"Since none of the aforementioned families alone could describe such conditional dependency, we combined multiple copulas into a linear mixture model\": Is it possible to elaborate a little bit more on whether this is actually the best technical choice here? For instance are there any side effects (e.g. computational) resulting from using this linear mixture? \n\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Summary: The authors built on existing work of GP vine copula
models. Some modifications are made, to conditional marginals and
mixing. Applications to mutual information estimation are discussed
and evaluated, and the approach is applied to joint
neural/behavioral data.


Discussion:
Strengths mentioned in the reviews are that the
application is (from a neuroscience perspective) interesting, that
estimating mutual information is an important problem, and that the
paper is very well-written. Weaknesses are the limited novelty (from a
machine learning perspective), and weak empirical validation. 

The authors have responded in detail, and were able to clarify a
number of unclear points. Clearly, however, the main criticisms noted
above are hard to address in discussion.

Despite the paper being overall clearly written, I agree with
reviewers that it is hard to tell from
abstract and introduction where the paper is going (even after
modifications made by the authors in the course of the discussion); of
the fairly long abstract, just about half a sentence relates to where
the proposed model differs from previous work. 


Recommendation:
I recommend rejection. Despite some clearly positive aspects, the two main criticisms voiced
by reviewers are serious: Weak validation and minimal 
novelty from a machine learning perspective. I agree that the
neuroscience application may be interesting, but requires more validation.

If the authors want to pursue this work further, I would suggest to
perhaps consider first where to position the paper's focus.
Estimation of mutual information is
a problem that is both hard and important. Any progress here would be
welcome, and simple usefulness could offset any lack of model novelty,
but it would have to be carefully and comprehensively
evaluated. On the other hand, a focus on neuroscience applications would
require more emphasis on, and presumably more space in the paper for,
relevant experiments. 

Decision: Reject