arXiv:1711.05869v1 [stat.ML] 16 Nov 2017

Predictive Independence Testing,
Predictive Conditional Independence Testing,
and Predictive Graphical Modelling
Samuel Burkart∗ and Franz J. Király†
Department of Statistical Science, University College London,
Gower Street, London WC1E 6BT, United Kingdom
November 17, 2017

Abstract
Testing (conditional) independence of multivariate random variables is a task central to statistical inference and modelling in general - though unfortunately one for which to date there
does not exist a practicable workflow. State-of-art workflows suffer from the need for heuristic or
subjective manual choices, high computational complexity, or strong parametric assumptions.
We address these problems by establishing a theoretical link between multivariate/conditional
independence testing, and model comparison in the multivariate predictive modelling aka supervised learning task. This link allows advances in the extensively studied supervised learning workflow to be directly transferred to independence testing workflows - including automated tuning
of machine learning type which addresses the need for a heuristic choice, the ability to quantitatively trade-off computational demand with accuracy, and the modern black-box philosophy for
checking and interfacing.
As a practical implementation of this link between the two workflows, we present a python
package ’pcit’, which implements our novel multivariate and conditional independence tests, interfacing the supervised learning API of the scikit-learn package. Theory and package also allow
for straightforward independence test based learning of graphical model structure.
We empirically show that our proposed predictive independence test outperform or are on par
to current practice, and the derived graphical model structure learning algorithms asymptotically
recover the ’true’ graph. This paper, and the ’pcit’ package accompanying it, thus provide powerful, scalable, generalizable, and easy-to-use methods for multivariate and conditional independence
testing, as well as for graphical model structure learning.

∗ samuel.burkart.16@ucl.ac.uk
† f.kiraly@ucl.ac.uk

1

Contents
1 Introduction
1.1 Setting: testing independence . . .
1.2 Predictive independence testing . .
1.3 Graphical model structure learning
1.4 Principal contributions . . . . . . .
1.5 Paper overview . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

2 Statistical independence testing
2.1 About (in)dependence . . . . . . . . . . .
2.2 Statistical independence testing . . . . . .
2.3 The state-of-art in advanced independence
2.4 Issues in current methodology . . . . . . .

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

4
4
4
5
5
5

. . . . .
. . . . .
testing
. . . . .

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

7
7
7
8
10

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

3 Predictive independence testing
3.1 Mathematical setting . . . . . . . . . . . . .
3.2 Elicitation by convex losses . . . . . . . . .
3.3 Predictive uninformedness . . . . . . . . . .
3.4 Statistical dependence equals predictability
3.5 Conditional independence . . . . . . . . . .
3.6 Testing prediction error against baseline . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

11
11
12
13
14
18
19

4 Graphical Models
4.1 Informal definition of graphical models . . .
4.2 Types of graphical models . . . . . . . . . .
4.2.1 Bayesian Networks . . . . . . . . . .
4.2.2 Markov Networks . . . . . . . . . . .
4.3 Graphical model structure learning . . . . .
4.3.1 Score-based methods . . . . . . . . .
4.3.2 Independence testing based methods

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

22
22
22
22
24
25
25
26

5 Predictive inference algorithms
5.1 Predictive conditional independence testing . . . . . . . . .
5.1.1 False-discovery rate control . . . . . . . . . . . . . .
5.1.2 Improving the prediction functionals . . . . . . . . .
5.1.3 Supervised learning for independence testing . . . .
5.2 Predictive structure learning of undirected graphical models

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

.
.
.
.
.

27
27
28
28
29
29

6 pcit package
6.1 Overview . . . . . . . . . . . . . . . .
6.1.1 Use cases . . . . . . . . . . . .
6.1.2 Dependencies . . . . . . . . . .
6.2 API description . . . . . . . . . . . . .
6.3 Function signatures . . . . . . . . . . .
6.4 API design . . . . . . . . . . . . . . .
6.4.1 Sklearn interface . . . . . . . .
6.4.2 Wrapper for Sklearn estimators
6.5 Examples . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

31
31
31
31
31
32
34
34
34
35

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

2

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

7 Experiments
7.1 Performance tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.1.1 Performance of conditional independence test . . . . . . . . . . . .
7.1.2 Performance of structure learning algorithm: Error rates . . . . . .
7.1.3 Performance of structure learning algorithm: Variance of estimated
7.2 Experiments on real data sets . . . . . . . . . . . . . . . . . . . . . . . . .
7.2.1 Sklearn data sets: Boston Housing and Iris . . . . . . . . . . . . .
7.2.2 Key short-term economic indicators (UK) . . . . . . . . . . . . . .

. . . .
. . . .
. . . .
model
. . . .
. . . .
. . . .

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

36
36
36
37
40
41
41
41

8 Discussion

43

Appendices

45

A Best uninformed predictors: classification
A.1 Misclassification loss is a probabilistic loss . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 Logarithmic loss is a proper loss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 Brier loss is a proper loss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45
45
46
46

B Elicited statistics for regression losses
B.1 Squared loss elicits the mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
B.2 Quantile loss elicits the quantile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47
47
48

3

1. Introduction
1.1. Setting: testing independence
The study of dependence is at the heart of any type of statistical analysis, and independence testing
is an important step in many scientific investigations, be it to determine if two things are related, or
to assess if an intervention had the desired effect, such as:
• When conducting market research, one might be interested in questions such as “Are our advertisement expenditures independent of our profits?” (hopefully not), or the more sophisticated
version “conditional on the state of the market, are advertisement expenditures independent of
our profits?”, which, if found to be true, would mean we are unlikely to increase profits through
an increase in our advertising budget (subject to the usual issue that inferring causality from
data requires an intervention or instrument).
• When collecting data for a medical study on the occurrence of an outcome Y , one might ask
“In the presence of data about attributes A for the subjects, should we still collect data for
attributes B ”. If Y is independent of attributes B given attributes A, additionally collecting
information about attributes B will not improve the knowledge of the state of Y (subject to the
usual issue that this conclusion is valid only for patients sampled in the same way from the same
population).
The difficulty of the independence testing task crucially relies on whether the following two complications present:
• Multivariate independence testing. This concerns the type of values which the involved
variables, in the second example the values which attributes A and attributes B may take: if the
domain of possible values consists either of a single number (continuous variable), or one class
out of many (categorical variable), the hypothesis test is “univariate”, and powerful methodology
exists that deals well with most scenarios, subject to some constraints. Otherwise, we are in the
“multivariate hypothesis testing” setting.
• Conditional independence testing. Whether there are conditioning random variables which
are to be controlled for, in the sense of testing independence conditional on a possible third
attributes C. If so, we are in the “conditional hypothesis testing” setting.
For the task which is neither multivariate nor conditional, well-recognized and universally applicable
hypothesis tests (such as the t-test or chi-squared test) are classically known. The multivariate setting
and the conditional setting are less well studied, and are lacking approaches which are general and
universally accepted, due to difficulties in finding a suitable approach which comes with theoretical
guarantees and is free of strong model assumptions. The setting which is both multivariate and
conditional is barely studied. The three major state-of-the-art approaches are density estimation,
copula and kernel based methods. Most instances are constrained to specific cases or rely on subjective
choices that are difficult to validate on real-world data. A more detailed overview of the state-of-art
and background literature is given in Section 2.

1.2. Predictive independence testing
The methodology outlined in this paper will consider multivariate and conditional independence testing from a new angle. The underlying idea for the test is that if two random variables X, Y are
independent, it is impossible to predict Y from X - in fact it will be shown that these two properties
are equivalent (in a certain quantitative sense). The same applies to conditional independence tests:
two random variables X, Y are conditionally independent given Z, if adding X as predictor variable
4

above Z will not improve the prediction of Y . In both cases, the predictive hypothesis test takes the
form of comparing a good prediction strategy against an optimal baseline strategy via a predictive
loss function. By determining if losses stemming from two predictions are significantly different, one
can then test statistically if a variable adds to the prediction of another variable (potentially in the
presence of a conditioning set), which by the equivalence is a test for (conditional) independence.

1.3. Graphical model structure learning
Probabilistic graphical models are a concept that heavily relies on independence statements for learning and inference. Most structure learning algorithms to date are, as a result of the lack of a scalable
conditional independence tests and additional combinatorial issues, constraint-based, or make heavy
assumptions on underlying distributions of a sample. This paper will leverage the predictive independence test into a new routine to estimate the undirected graph for the distribution underlying a
sample, based on conditional independence testing, allowing it to make only weak assumptions on the
underlying distribution.

1.4. Principal contributions
The new approach to multivariate and conditional independence testing outlined in this paper improves
concurrent methodology by deriving an algorithm that
• features a principled model selection algorithm for independence testing by linking the field of
independence testing to the field of predictive modelling, thus filling a gap in state-of-the-art
methodology,
• additionally allowing independence testing to directly benefit from the well-understood and
efficiently implemented theory for model selection and parameter tuning in predictive modelling,
• is comparatively fast and scalable on a wide variety of problems and
• deals with the multivariate and conditional independence testing task in a straightforward manner.
Additionally, an algorithm leveraging the newly derived test into a scalable independence testingbased graphical model structure learning algorithm is outlined, which overcomes issues in the field by
offering a test for undirected graph structure learning that
• has polynomial time complexity, O(p2 ) in the number of variables, p,
• is exact and
• offers stringent methodology to control the number of type 1 errors in the estimated graph.

1.5. Paper overview
Section 2 will provide an overview of the most important tasks in and approaches to statistical independence testing and outline the issues in current methodology. Section 3 will then propose and derive
a novel approach to independence testing, the predictive conditional independence test (PCIT). After,
section 4 will introduce the concept of a graphical model and survey the current structure learning
algorithms. Section 5 will then state the relevant algorithms of the Python implementation, which
is outlined in section 6. Section 7 provides performance statistics and examples for the conditional
independence test as well as the graphical model structure learning algorithm. Lastly, section 8
will describe the advantages of using the PCIT for independence testing and outline drawbacks and
directions for further research.

5

Authors’ contributions
This manuscript is based on SB’s MSc thesis, submitted September 2017 at University College London
and written under supervision of FK, as well as on an unpublished manuscript of FK which relates
predictive model selection to statistical independence. The present manuscript is a substantial reworking of the thesis manuscript, jointly done by SB and FK.
FK provided the ideas for the independence tests (paper) in section 3 and the usage of them for
graphical models, SB and FK jointly conceived the ideas for the graphical model structure learning algorithm. Literature overview is due to SB with helpful pointers by FK. Python package and
experiments are written and conducted by SB, partially based on discussions with FK.

Acknowledgements
We thank Fredrik Hallgren and Harald Oberhauser for helpful discussions.
FH has briefly worked on a thesis about the same topic under the supervision of FK, before switching
to a different topic. While FH did, to our knowledge, not himself make any of the contributions found
in this paper, discussions with him insipred a few of them.
HO pointed FK in generic discussions about loss functions towards some prior results (elicitation of
median, and elicitation as a concept as defined by Gneiting and Raftery [17]) which helped making
some of the proven statements more precise (such as the Q-losses being faithful for the univariate real
case).

6

2. Statistical independence testing
2.1. About (in)dependence
Statistical (in)dependence is a property of a set of random variables central to statistical inference.
Intuitively, if a set of random variables are statistically independent, knowing the value of some will
not help in inferring the values of any of the others.
Mathematically: let X, Y and Z be random variables taking values in X, Y and Z respectively. As
usual, we denote for A ∈ X by P (X ∈ A) the probability of X taking a value in A, and for A ∈ X, C ∈ Z
by P (X ∈ A|Z ∈ C) := P (X ∈ A, Z ∈ C)/P (Z ∈ C) be the conditional probability of X taking a
value in A when Z is known/observed to have taken a value in C.
Definition 2.1. X and Y are called marginally independent (of each other) if for all A ⊆ X and
B ⊆ Y (where the below expression is defined) it holds that
P (X ∈ A, Y ∈ B) = P (X ∈ A)P (Y ∈ B).
This formulation allows for X and Y to be defined over sets of random variables that are a mixture of
continuous and discrete, as well as being univariate or multivariate (thus implicitly covering the case
of multiple univariate random variables as well).
Definition 2.2. X and Y are called conditionally independent (of each other) given Z if for all
A ⊆ X, B ⊆ Y, and C ⊆ Z (where the below expression is defined) it holds that
P (X ∈ A, Y ∈ B|Z ∈ C) = P (X ∈ A|Z ∈ C)P (Y ∈ B|Z ∈ C).
For absolutely continuous or discrete X and Y , Definition 2.1 straightforwardly implies that marginal
independence is equivalent to the joint distribution or mass function factorizing, i.e., it equals the
product of the marginal distributions’ probability or mass function. The analogue result also shows
for the conditional case.
We would like to note that the mathematical definition independence is symmetric, i.e., having the
property from Definition 2.1 is unaffected by interchanging X and Y . In contrast, the intuitive motivation we gave at the beginning however, namely that knowing the value of X does not yield any
additional restriction regarding the value of Y , is non-symmetric. This non-symmetric statement of
statistical independence is commonly phrased mathematically as P (Y ∈ B|X ∈ A) = P (Y ∈ B)
which is directly implied by Definition 2.1 (and that of the conditional).
This non-symmetric characterization of statistical independence is also morally at the heart of the
predictive characterization which we later give in the supervised learning setting: Y cannot be predicted from X, phrased in a quantitative way that connects checking of statistical independence to
estimation of a supervised generalization error difference which eventually allows for a quantitatve
testing of the hypothesis of independence.

2.2. Statistical independence testing
Testing whether two random variables are (conditionally) independent is one of the most common
tasks in practical statistics (perhaps the most common one after summarization) and one of the most
important topics in the theoretical study of statistical inference (perhaps the most important one).
The most frequently found types of quantitative certificates for statistical independence are phrased in
the (frequentist) Neyman-Pearson hypothesis testing paradigm. In terms of the methodological/mathematical idea, the main distinction of independence tests is into parametric (where the data is assumed
to stem from a certain model type and/or distribution) and non-parametric (”distribution-free”) tests.
The hypothesis tested usually takes one of the three following forms:
7

|=

Marginal independence, X

Y

Conditional independence, X

|=

Marginal independence of two random variables is the topos of much of classical statistics. Well
known classical test statistics for the continuous univariate case are the Pearson correlation coefficient,
Kendall’s τ and Spearman’s ρ. For discrete variables, Pearson’s χ-squared statistic can be used. From
theoretical results on these statistics’ asymptotic or exact distribution, univariate independence tests
may be derived. For X being discrete and Y being continuous, t-family tests or Wilcoxon family tests
may be used.
Testing with at least one of X or Y being multivariate is much more difficult and less standard, current
methodology will be outlined in the next section.
Y |Z

Testing for conditional independence is an inherently more difficult problem than marginal independence [7]. Hence, few viable options exist to date to test if two sets of variables are conditionally
independent given a conditioning set. The most common instances are either based on parametric
model assumptions (e.g., linear) or on binning Z and X, then comparing if the distribution of Y
conditioned on Z changes, when additionally conditioning on X. Strategies for of conditional independence testing may demand much larger sample sizes than marginal independence tests due to the
explicitly modelling sub-domains of the values which Z takes.
d

Equality of distribution of two unpaired samples, X1 = X2 ,
d

where = indicates equality in underlying distribution. An unpaired two-sample test for equality
of distribution tests whether two unpaired i.i.d. samples from X1 and X2 are sampled from the
same underlying distributions. The connection to independence testing is not obvious but may be
established by pairing each draw from X1 or X2 with a draw from a variable Y taking values in {1, 2},
indicating whether the draw was from X1 or X2 , i.e., taking value i if the draw was from Xi .

2.3. The state-of-art in advanced independence testing
The task of independence testing (in various settings) has been tackled from many different angles.
As soon as one considers observations that are multivariate, or the conditional task, there is not one
universally agreed-upon method, but many different approaches. The most prominent ideas used for
multivariate and/or conditional independent testing will be presented in this section together with
their advantages and shortcomings.
Density estimation
The classical approach. In case of existence, the joint probability density function contains all the
necessary information about independence structures for a set of random variables in the multivariate
case. While many univariate tests are based on or may be interpreted as being based on some sort of
density estimation, for the multivariate and conditional case, density estimation is a difficult problem.
One example is shown in [12], where density estimation-based information-theoretical measures are
used to conduct multivariate marginal independence tests.
Copulas
An approach that is widely used in finance and risk management. Copulas are multivariate probability
distributions with uniform marginals. To use copulas for independence tests, one transforms the
marginals of a multivariate distribution into uniform distributions, and the resulting copula contains
all information about the independence structure. A simple example using the empirical cdf can be
found in [16]. Copulas have mainly been used for marginal multivariate independence tests, such as in
8

[34]. The task of testing for conditional independence has been attempted through the use of partial
copulas [6], however strong assumptions (linearity) are made on the type of relationship between the
variables that are to be tested for independence and the conditioning set. Two-sample testing has
largely remain unaddressed, since its application in the financial sector are less relevant. [33] describes
a two-sample test for the estimated copulas, and hence, the independence structures, but the test does
not extend to a comparison of the joint distributions as a whole.
Kernel methods

|=

A relatively recent approach using optimization techniques in reproducing kernel Hilbert spaces to
answer independence queries based on various test-statistics. Kernel methods were first used for
marginal multivariate independence testing via a quantitative measure of dependence, the HilbertSchmidt Independence Criterion (HSIC), which is defined as the squared Hilbert Schmidt norm of
the cross-covariance operator (the covariance between two difference spaces), see [19]. [20] further
expands on the theoretical foundations of the HSIC’s distribution under the null, resulting in a hypothesis test for the independence of two variables X and Y using the HSIC criterion. [18] outlines
a non-parametric test for independence, and demonstrates its performance on samples of random
variables with up to three dimensions.
The conditional case is tackled in [37], which derives a test statistic based on conditional crosscovariance operators, for which the asymptotic distribution under the null-hypothesis X Y |Z is
derived. They state that the manual choice of kernel can affect type II and, more importantly, type I
errors, especially when the dimensionality of the conditioning set is large. Additionally, they found
that the performance of their method decreases in the dimensions of the conditioning set, constraining
the set of problems for which the test is viable for.
As for two-sample testing, [3] derive a test based on the Euclidean inter-point distances between the
two samples. A different approach to the same test statistics with the additional use of characteristic
functions was made by [13]. While the difficulty of density estimation is reduced when using the
empirical distributions, the data requirements tend to be much larger, while additionally imposing a
sometimes taxing computational complexity. Both of these methods can be related to kernel functions, which was picked up in [22], which proposed kernel-based two-sample test by deriving three
multivariate tests for assessing the null-hypothesis that the distributions p and q generating two samples X ∼ p, Y ∼ q are equal, p = q. The criterion introduced is the “Maximum Mean Discrepancy”
(MMD) between p and q over a function space F.
M M D = sup (EX∼p [f (X)] − EY ∼q [f (Y )])
f ∈F

That is, the maximum difference between the expected function values with respect to the function
space (in their case, a Reproducing Kernel Hilbert Space).
Predictive independence testing: precursors
There are a few instances in which independence testing has already been approached from the perspective of predictability: most prominently, Lopez-Paz and Oquab [28] present a two-sample test
(on equality of distribution) based on abstract binary classifiers, e.g., random forests. The presented
rationale is in parts heuristic and specific to the two-sample setting. Their main idea is as follows: if
the null-hypothesis of the two samples being drawn from the same distribution is true, no classifier
assigning to a data point the distribution from it was drawn should fare significantly better than
chance, on an unseen test set.
The work of Lopez-Paz and Oquab [28] may be seen as a special case of our predictive independence
test presented later, for Y taking values in a two-element set, and the loss being the misclassification
loss. It is probably the first instance in which the abstract inability to predict is explicitly related to
independence, albeit somewhat heuristically, and via the detour of two-sample testing for equality in
9

distribution.
In earlier work, Sriperumbudur et al. [35] already relate the kernel two-sample test to a specific classifier, the Parzen window classifier, evaluated by the misclassification loss. Thus this ideas of Sriperumbudur et al. [35] may in turn be seen as a precursor of the later work of Lopez-Paz and Oquab [28]
which abstracts the idea that any classifier could be used to certify for equality of the two samples.

2.4. Issues in current methodology
Density estimation
While the probability density function is in some sense optimal for measuring dependence, since
it contains all the information about the random variable, its estimation is a difficult task, which
requires either strong assumptions or large amounts of data (or both). Furthermore, due to the
curse of dimensionality, density estimation for 3 or more dimensions/variables is usually practically
intractable, hence its practical usefulness is often limited to testing pairwise independence rather than
testing for full (mutual) independence.
Copulas
Leaving aside issues arising by practitioners misunderstanding the method (which had a strong contribution to the 2007/2008 financial crisis, but no bearance whatsoever on the validity of the method
when applied correctly), copula-based independence testing is largely a heuristics driven field, requiring many subjective manual choices for estimation and testing. Above all, copula methods require a
user to subjectively choose an appropriate copula from a variety of options, such as for example the
survival copula, the multivariate Gaussian or the multivariate Student’s t-copula [9]. Additionally,
two-sample testing is to date largely unaddressed in the field.
Kernels
As for the copula-based methods, kernels require many subjective and manual choices, such as the
choice of kernel function, and its hyper-parameters. While in a theoretical setting, these choices are
made by using the ground truth (and generating artificial data from it), in practice it is difficult to
tune the parameters and make statements about the confidence in results. Additionally, the cost of
obtaining the test statistic and its asymptotic distribution may be high for many of the proposed approaches. There are attempts at resolving these issues, [36] and [23] outline heuristics and strategies
to minimize the heuristic choices, and [21] and [10] propose more computationally efficient strategies.

While all methods have their merits, and there is vast research on specific problems, providing specific
solutions to the individual challenges, what is missing is an approach that is scalable and powerful
not just for specific cases, but for the general case, that automatically finds good solutions to real
problems, where the ground truth is not known and cannot be used for tuning.

10

3. Predictive independence testing
This section will explore the relationship between predictability and dependence in a set of variables. It
will be shown that one can use supervised learning methodology to conduct marginal and conditional
independence tests. This distinction between marginal and conditional dependence is not theoretically
necessary (since the marginal case can be achieved by setting the conditioning set to the empty set),
but is made to highlight specific properties of the two approaches. First, equivalence of independence
statements and a specific supervised learning scenario will be shown. After, a routine leveraging this
equivalence to test for conditional independence will be proposed.

3.1. Mathematical setting
The independence tests will be based on model comparison of supervised learning routines.
Supervised learning is the task where given i.i.d. samples (X1 , Y1 ), . . . , (XN , YN ) ∼ (X, Y ) taking
values in X × Y, to find a prediction functional f : X → Y such that f (Xi ) well approximates a
target/label variables Yi , where ”well” is defined with respect to a loss function L which is to be
minimized in expectation.
Definition 3.1. A (point prediction) loss functional is an element of [Y×Y → R], i.e., a function with
range Y × Y and image R. By convention, the first argument will be considered the proposed/predicted
value, the second the true value to compare to. By convention, a lower loss will be considered more
beneficial.
A loss functional L : Y × Y → R is called:
(i) convex if L(., y) is lower bounded and E [L(Y, y)] ≤ L (E[Y ], y) for any Y-valued random variable
Y and any y ∈ Y, and
(ii) strictly convex if L(., y) is lower bounded and E [L(Y, y)]
L (E[Y ], y) for any non-constant
Y-valued random variable Y and any y ∈ Y.
More formally, a good prediction functional possesses a small expected generalization error
L (f ) := E[L(f (X), Y )].
In usual practice and in our setting f is estimated from training data
D = {(Xi , Yi )}N
i=1 ,
hence it is in fact a random object. The generalization error may be estimated from test data
T = {(Xi∗ , Yi∗ )}M
i=1

∗
∗
where (X1∗ , Y1∗ ), . . . , (XM
, YM
) ∼ (X, Y ),

which we assume i.i.d. and independent of the training data D and f , as
b
L (f ) =

M
X

[L(f (Xi∗ ), Yi∗ )].

i=1

Independence of training data D and test data T is required to avoid bias of the generalization error
estimate (“over-fitting”). In the following sections, when estimating the generalization error, it will
be assumed that f is independent of the test data set, and we condition on D, hence treat f as fixed
(i.e., conditional on D).
In the following, we assume that Y is a vector space, and Y ⊆ Rq , which naturally includes the setting
of supervised regression. However, it also includes the setting of supervised classification through a
special choice of loss:

11

Remark 3.2. Deterministic classification comprises the binary case of Y = {−1, 1}, or more generally
Y being finite. In neither case are additions or expectations defined on Y, as Y is just a discrete
set. Supervised classification algorithms in the deterministic case are asked to produce a prediction
functional f : X → Y, and are usually evaluated by the misclassification loss or 0/1-loss
L : Y × Y → R : (y, y∗ ) 7→ 1[y 6= y∗ ],
where 1[y 6= y∗ ] is the indicator function which evaluates to 0 if y = y∗ and to 1 otherwise. Hence,
Definition 3.1 of convexity does not directly apply to the misclassification loss as expectations are not
defined on Y.
However, by allowing all algorithms to make predictions which are Y-valued random variables (instead
of constants), one finds that
E[L(Y, y)] = L0 (p, y)

where L0 : (pY , y) 7→ 1 − pY (y)

and pY is the probability mass function of Y . Identifying Y-valued random variables with their probability mass functions, one may replace
(i) Y by the corresponding subset Y0 of R#Y−1 which is the set of probability vectors (the so-called
probability simplex). For example, Y = {−1, 1} would be replaced by [0, 1].
(ii) the misclassification loss by the probabilistic, convex (but not strictly convex) loss L0 : Y0 ×Y0 → R,
where the observations in the second Y0 argument are always pmf describing constant random
variables.
A further elementary calculation (see appendix A.1 for an explicit derivation) shows that L0 is always
minimized by making deterministic predictions:
argmin E[L0 (p, Y )] ∩ [Y → {0, 1}] is non-empty,
p∈Y0

i.e., the L0 -best classifier may always be chosen to be a deterministic one, i.e., one that always predicts
a probability mass functions with probabilities 0 or 1.
This exhibits deterministic classification as a special case of probabilistic classification with a special
choice of loss function.

3.2. Elicitation by convex losses
Loss functionals are canonically paired with summary statistics of distributions:
Definition 3.3. Let L : Y × Y → R be a convex loss functional. For a Y-valued random variable Y ,
we define
µL ([Y ]) := argminE [L(y, Y )]
y∈Y

where [Y ] denotes Y as a full object (i.e., a measurable function), rather than its value.
Following Gneiting and Raftery [17], we term the functional which maps Y-valued distributions to
sub-sets of Y the eliciting functional associated to L. We call µL ([Y ]) the summary of Y elicited by
L.
Note that well-definedness, i.e., existence of the minimizer, is ensured by convexity of L (and the
implied continuity). If L is strictly convex, there is furthermore a unique minimizer, in which case we
will exchangeably consider µL to be a functional with target Y.
Well-known examples of elicited summaries are given in the following:
Lemma 3.4. The following relations between losses and elicited statistics of real-valued random variables hold:
12

(i) the (strictly convex) squared loss L : (y, y∗ ) 7→ (y−y∗ )2 elicits the mean. That is, µL ([Y ]) = E[Y ]
for any Rn -valued random variable Y .
(ii) the (convex but not strictly convex) absolute loss L : (y, y∗ ) 7→ |y − y∗ | elicits the median(s).
That is, µL ([Y ]) = median[Y ] for any R-valued random variable Y .
(iii) the (convex but not strictly convex) quantile-loss (or short: Q-loss) L(y, y∗ ) = α · m(y∗ , y) + (1 −
α) · m(y, y∗ ), with m(x, z) = min(x − z, 0), elicits the α-quantile(s). That is, µL ([Y ]) = FY−1 (α)
for any R-valued random variable Y , where FY−1 : [0, 1] → P (R) is the set-valued inverse c.d.f. of
Y (with the convention that the full set of inverse values is returned at jump discontinuities rather
than just an extremum).
Proof. (i) after substitution of definition, the claim is equivalent to the statement to the mean being
the minimizer of squared distances. A more explicit proof is given in Appendix B.1.
(ii) follows, by setting α = 12 , from (iii).
(iii) This is carried out in Appendix B.2.
In the supervised setting, the best possible prediction functional can be exactly characterized in terms
of elicitation:
Proposition 3.5. Let L be a (strictly) convex loss. Then, it holds that
argmin εL (f ) = [x 7→ µL [Y |X = x]] .
f ∈[X→Y]

That is, the best possible prediction as measured by L is predicting the statistic which L elicits from
the conditional random variable Y |X = x.
Proof. The prediction functional ω : x 7→ µL [Y |X = x] is well-defined, hence it suffices to show that
εL (f ) ≥ εL (ω) for any prediction functional f : X → Y.
Now by definition of µL , it holds that
E [L(ω(X), Y )|X] ≤ E [L(f (X), Y )|X] .
Taking total expectations yields the claim.
Intuitively, the best prediction functional, as measured by a convex loss L, always predicts the statistic
elicited by L from the conditional law [Y |X = x].

3.3. Predictive uninformedness
We will now introduce the notion of an uninformed baseline which will act as a point of comparison.
Definition 3.6. A prediction functional f : X → Y is called uninformed if it is a constant functional,
i.e., if f (x) = f (y) for all x, y ∈ X. We write uα for the uninformed prediction functional uα : x 7→ α.
We will show that, for a given loss function, there is one single uninformed baseline that is optimal.
Lemma 3.7. Let L be a (strictly) convex loss, let µ := µL ([Y ]) be the/a statistic elicited by L (see
Definition 3.3). Then, the following quantities are equal:
(i) inf{εL (f ) : f is an uninformed prediction functional}
(ii) εL (uµ )
That is, uµ is achieves the lowest possible (L-)loss amongst uninformed prediction functionals and
prediction strategies.
Proof. Note that ε(uα ) = E[L(α, Y )]. It follows hence by definition of µ that ε(uµ ) ≤ ε(uα ) for any
(constant) α ∈ Y. I.e., uµ is the best uninformed prediction functional.
13

Lemma 3.7 motivates the definition of the best uninformed predictor:
Definition 3.8. We call uµ , as defined in Lemma 3.7, the (L-)best uninformed predictor (even though
it may not be unique, the choice, when possible, will not matter in what follows).
We call a prediction functional (L−)better-than-uninformed if its expected generalization loss is strictly
smaller than of the (L−)best uninformed predictor. More formally, a prediction functional f is Lbetter-than-uninformed if εL (f ) ≤ εL (uµ ).
For convenience, we further introduce some mathematical notation for best predictors:
Notation 3.9. Let L be a (strictly) convex loss. We will write:
(L)

(i) ωY
(ii)
(L)

(L)
ωY |X

:= [x 7→ µL [Y ]] for the/a (L-)best uninformed predictor as defined in Definition 3.8.
:= [x 7→ µL [Y |X = x]] for the/a (L-)best predictor as considered in Proposition 3.5.
(L)

ωY and ωY |X are unique when L is strictly convex, as per the discussion after Definition 3.3. When
multiple choices are possible for the minimizer, i.e., if L is convex but not strictly convex, an arbitrary
choice may be made (not affecting subsequent discussion). The superscript L may be omitted in
situations where the loss is clear from the context.
An important fact which we will use in testing is that if a better-than-uninformed prediction functional
exists, then ωY |X is an example:
Proposition 3.10. Fix a (strictly) convex loss. The following are equivalent:
(i) ωY |X is L-better-than-uninformed.
(ii) ωY |X is not L-uninformed.
(iii) There is an L-better-than-uninformed prediction functional.
Note that equivalence of (i) and (ii) is not trivial, since there are prediction functionals which are not
better-than-uninformed but not uninformed.
Proof. The equivalence of (i) and (iii) follows directly from Lemma 3.7, and noting that (i) implies
(iii).

(i)⇔(ii): By Proposition 3.5, ε ωY |X ≤ ε(f ) for any f , in particular also any uninformed f . By
Lemma 3.7 and the above, the inequality is strict if and only if ωY |X is better-than-uninformed.

3.4. Statistical dependence equals predictability
We continue with a result that relates - or more precisely, equates - better-than-uninformed predictability with statistical dependence. As such, it shows that the choice of constant functions as a
proxy for uninformedness was canonical.
We start with the more intuitive direction of the claimed equivalence:
Proposition 3.11. Let X, Y be random variables taking values in X, Y. Let L be any convex loss
functional. If X, Y are statistically independent, then: There exists no L-better-than-uninformed
prediction functional for predicting Y from X. More precisely, there is no prediction functional f :
X → Y and no convex loss L : Y × Y → R such that f is L-better-than-uninformed.
Proof. Assume X, Y are statistically independent. Let L be a convex loss function, let f : X → Y be
a prediction functional. Then, by convexity of L,
E [L(f (X), Y )|Y ] ≥ L(E[f (X)], Y )|Y.

14

Since X is independent of Y (by assumption) and f is not random, it holds that E[f (X)]|Y = E[f (X)],
i.e., E[f (X)]|Y = y, as a function of y, is constant. Writing ν := E[f (X)], we hence obtain that
L(E[f (X)], Y )|Y = L(ν, Y ) = L(uν , Y ).
After taking total expectations, it hence holds by the law of total expectation that
E [L(f (X), Y )] ≥ E [L(uν , Y )] ,
meaning that f is not better-than-uninformed w.r.t. L. Since f and L were arbitrary, the statement
holds.
Proposition 3.11 states that independence implies unpredictability (as measured per expected loss). A
natural question to ask is whether the converse holds, or more generally which converses exactly, since
the loss functional L in (ii) and (iii) of Proposition 3.10, was arbitrary, and independence in (i) is a
statement which remains unchanged when exchanging X and Y , while predictability is not symmetric
in this respect. Hence the weakest possible converse would require unpredictability w.r.t. all convex
losses and w.r.t. either direction of prediction, however much stronger converses may be shown.
Before stating the mathematically more abstract general result for the converse direction, we first
separately present special cases for the three important sub-cases, namely deterministic classification,
probabilistic classification, and regression.
Theorem 1. As in our setting, consider two random variables X, Y , taking values in X and Y.
Assume Y = {−1, 1}. The following are equivalent:
(i) X and Y are statistically independent.
(ii) There exists no better-than-uninformed prediction functional predicting Y from X. More precisely, there is no prediction functional f : X → Y such that f is L-better-than-uninformed for
the misclassification loss L : (p, y) 7→ 1 − p(y).
Regarding the specific loss L, see Remark 3.2 for the identification of a class prediction with a
probability-1 probabilistic prediction.
Proof. By Proposition 3.11, the only direction which remains to be proven is (ii)⇒(i): (i)⇒(ii) follows directly from substituting the specific L into the implication between statements with the same
numbers in Proposition 3.11.
(ii)⇒ (i): We prove this by contraposition: we assume X and Y are statistically dependent and
construct a better-than-uninformed prediction functional.
By the equivalence established in Remark 3.2, the best uninformed predictor ωY is predicting the
most probable class, w.l.o.g. 1 ∈ Y, and its expected generalization loss is one minus its generative
frequency ε(ωY ) = P (Y = −1).
Since X and Y are statistically dependent, there is a positive probability X0 ⊆ X (measurable with
positive probability measure) such that P (Y = 1|X ∈ X0 ) ≥ P (Y = −1|X ∈ X0 ) (the definition yields
6=, but w.l.o.g. replacing X0 by its complement yields ≥). An elementary computation shows that the
prediction functional f : X → Y : 2 · 1[x ∈ X0 ] − 1 is better-than-uninformed.
The proof of Theorem 1 shows that for a converse, one does not necessarily need to consider predictability of X from Y as well. However, the deterministic binary classification case is somewhat
special, since the misclassification loss is insufficient in the multi-class case, and in general a single
loss function will be unable to certify for independence. In order to formulate these negative results,
we define shorthand terminology for stating correctness of the converse.

15

Definition 3.12. Fix a (label) domain Y, let L be a set with elements being convex loss functionals
in [Y × Y → R]. We call L faithful for Y if for every statistically dependent pair of random variables
X, Y taking values in X, Y, there is L ∈ L and a prediction functional f : X → Y such that f is
L-better-than-uninformed. If L is a one-element-set, we call its single element faithful for Y.
If L is faithful for Y, and L is endowed with a measure µ such that no µ-strict sub-set of L is faithful
for Y, then we call L (µ-)strictly faithful.
If Y is canonically associated with a prediction task (such as classification for finite Y or classprobabilities Y, regression for continuous Y), the reference to Y may be replaced by a reference to
that task.
In this terminology, Theorem 1 states that the misclassification loss is faithful for the label domain Y =
{−1, 1}, or equivalently that the misclassification loss is faithful for deterministic binary classification
(and strictly faithful, since any set smaller than an one-element set, as by the counting measure,
is empty and by definition not faithful). We can now state some negative and positive results on
deterministic multi-class classification and probabilistic classification:
Proposition 3.13. (i) The misclassification loss L : (p, y) 7→ 1 − p(y) is not faithful for deterministic multi-class classification, i.e., for Y being finite and containing 3 or more elements.
(ii) The log-loss L : (p, y) 7→ − log(p(y)) is (strictly) faithful for probabilistic classification.
(iii) The squared probabilistic classification loss L : (p, y) 7→ (1 − p(y))2 is (strictly) faithful for
probabilistic classification.
P
(iv) The Brier loss L : (p, y) 7→ (1 − p(y))2 + y0 6=y p(y 0 )2 is (strictly) faithful for probabilistic
classification.
Proof. (i): It suffices to construct a counterexample for each possible Y, i.e., every finite Y with 3 or
more elements. Let X = {−1, 1}, and let y ∈ Y be arbitrary. Define X, Y such that P (Y = y|X =
1) = P (Y = y|X = −1) := 0.9 and P (Y = y 0 |X = 1) 6= P (Y = y 0 |X = −1) for some class y 0 6= y.
This choice is possible as Y has 3 or more elements. The best uninformed predictor always predicts
y, with expected generalization loss 0.1, while it cannot be outperformed, see e.g. the discussion in
Remark 3.2.
(ii)-(iv): for faithfulness, it suffices to show that ε(ωY |X ) = ε(ωY ) implies statistical independence of
X, Y .
Explicit computations in each case (see Appendices A.2 and A.3) show that argminp E[L(p, Y )] = pY ,
where pY is the probability mass function of Y , implying that ωY = [x 7→ [y 7→ P (Y = y)]]. By
conditioning the same statement on X, this also implies that argminp E[L(p, Y )|X] = [y 7→ P (Y =
y|X)], thus argminf E[L(f (X), Y )] = [y 7→ P (Y = y|X = x)]. In particular, ωY |X = [x 7→ [y 7→
P (Y = y|X = x)]] is the unique minimizer of the expected generalization loss. Thus, ε(ωY |X ) = ε(ωY )
only if both functions are identical, i.e., P (Y = y|X = x) = P (Y = y) for all x, which is one possible
definition of X and Y being statistically independent.
For regression, the usual loss functions are unable to certify for independence anymore:
Proposition 3.14. The following convex loss functions (taken as single-element sets) are not faithful
for univariate regression, i.e., for the label domain Y = R
(i) the squared loss L(y, y∗ ) = (y − y∗ )2
(ii) the absolute loss L(y, y∗ ) = |y − y∗ |
(iii) the distance loss L(y, y∗ ) = d(y, y∗ )2 , for any metric d : Y × Y → R
Proof. It suffices to construct, for each convex loss functional L, a counterexample where X, Y are
statistically dependent, but no prediction functional predicting Y from X is better-than-L-uninformed.
In each case, one may construct two Y-valued random variables Y1 , Y2 with distinct laws such that
µL [Y1 ] = µL [Y2 ] - for example, an arbitrary non-constant Y1 and the Y2 being constant µL [Y1 ]. Further

16

setting X = {−1, 1} and defining a non-constant X-valued random variable X, together with an Yvalued random variable such that Yi = Y |X = i for i ∈ X yields an example of a statistically dependent
pair X, Y of random variables where the constant prediction of µL [Y1 ] = µL [Y2 ] is not only the L-best
uninformed prediction functional, but also the L-best prediction functional.
Using equivalence of (i) and (iii) in Proposition 3.10 proves the claim.
The previously introduced quantile losses form a strictly faithful set of losses for univariate regression:
Theorem 2. The set of Q-losses is strictly faithful for (univariate) regression.
More precisely, the set L = {Lα : α ∈ [0, 1]}, where Lα (y, y∗ ) = α · m(y∗ , y) + (1 − α) · m(y, y∗ ), with
m(x, z) = min(x − z, 0) (see Lemma 3.4 (iii)), endowed with the Lebesgue measure through canonical
identification of Lα with [0, 1], is strictly (Lebesgue-measure-)faithful for Y = R.
Proof. For faithfulness, we show that impossibility of (Lα )-better-than-uninformed prediction of a
R-valued random variable Y from an X-valued random variable X implies statistical independence of
X and Y .
Thus, assume that there is no α ∈ R and no prediction functional f such that f is Lα -betterthan-uninformed. By equivalence of (ii) and (iii) in Proposition 3.10 and negation/contraposition,
(Lα )
ωY |X
is uninformed hence constant. By Lemma 3.4 (iii), µLα [Y |X = x] is the α-quantile of the
conditional random variable Y |X = x, which by the previous statement does not depend on x. Since
α was arbitrary, none of the quantiles of Y |X = x depends on x, i.e., the cdf and hence the laws
of all conditionals Y |X = x agree, which implies (by one common definition/characterization of
independence) that X and Y are statistically independent.
Strict faithfulness follows by following through the above argument after removing a positive-measure
open set U ⊆ [0, 1] from the indices, i.e., Lα for α∈ U from L. As U has positive measure, we may pick
u ∈ U and X = {−1, 1} as well as conditional cdf such that P (Y ≤ u|X = 1) 6= P (Y ≤ u|X = −1)
while P (Y ≤ x|X = 1) = P (Y ≤ X|X = −1) for all x 6∈ U . By the above argument, predicting the
α-quantile of Y is the Lα -best prediction functional from X for any α 6∈ U , and it is furthermore a
uninformed prediction strategy, thus L is not faithful after removing Lα , α ∈ U .
In the light of Theorem 2, it may be interesting to ask for a theoretical characterization of a strictly
faithful set of losses for univariate regression (e.g., does it need to be infinite?), or what may be a set
of strictly faithful losses for multivariate regression.
A-priori, it is even unclear whether there is a set of (not necessarily strictly) faithful losses for general
prediction tasks, which the following result answers positively:
Theorem 3. Assume that Y may be identified (including the taking of expectations) with a sub-set
of Rn , for some n, e.g., in multivariate regression, or simultaneous prediction of multiple categorical
and continuous outputs. Then, the set of convex losses is faithful for Y.
Proof. Consider random variables X and Y taking values in X and Y, are statistically dependent.
By definition of (in)dependence, this is equivalent to there existing X0 ⊆ X, Y0 ⊆ Y (measurable
with positive joint probability measure) such that P (Y ∈ Y0 |X ∈ X0 ) 6= P (Y ∈ Y0 ). By taking
(or not taking) the complement of X0 within X, we may assume without loss of generality that
P (Y ∈ Y0 |X ∈ X0 ) P (Y ∈ Y0 ).
Define g : Y → R ; x 7→ x2 (where we have used the identification of Y with Rn . Define
L : Y × Y → R ; (y, y∗ ) 7→ g(y) · 1(y∗ ∈ Y0 ) + g(y − α) · 1(y∗ 6∈ Y0 ),
where 1(y∗ ∈ Y0 ) is the indicator function for y∗ ∈ Y0 , and α ∈ Y \ {0} is arbitrary. Define f : X →
Y ; x 7→ 0 if (x ∈ X0 ), otherwise α. An elementary computation shows that εL (f ) is better-thanuninformed, hence we have proved non-(ii).

17

For the general case, it seems interesting to ask what would be a strictly faithful set of losses, how
such sets may be characterized, or whether they even exist (which seems neither obviously nor directly
implied by the existence of a faithful set of losses).
Due to the constructive nature of respective proofs (or, more precisely, the semi-constructivity of
the proof for multi-variate regression in Theorem 3), model comparison procedures suggested by
Theorems 1 and 2 on univariate classification and regression will be used in the testing procedures.
For convenience we briefly repeat the main results used in testing as a corollary:
Corollary 3.15. Consider two random variables X, Y , taking values in X and Y, where Y is finitely
supported pmf (classification) or where Y ⊆ Rq (regression). The following are equivalent:
(i) X and Y are statistically dependent.
(L)

(L)

(ii) ε(ωY |X ) εL (ωY ) for L the log-loss/Brier-loss (classification), resp. for some convex loss L
(regression). I.e., the L-best predictor is L-better-than-uninformed, for some L.
(L)

(iii) there exists a prediction functional f : X → Y such that εL (f )
εL (ωY ) for L the logloss/Brier-loss (classification), resp. for some convex loss L (regression). I.e., there exists an
L-better-than-uninformed prediction functional, for some L.

|=

Since Corollary 3.15 (i) is the alternative hypothesis in an independence test, we may use a test for
the equivalent hypothesis Corollary 3.15 (iii), i.e., comparing the performance a learnt f with the
best uninformed baseline, as an independence test. Any choice of f is sufficient for the methodology
outlined in the latter Section 5.1. Since the null-hypothesis is that X Y , choosing a bad f or failing
to detect multivariate dependence (a specific case) only decreases the power of the test, while the type
1-error remains unaffected.

3.5. Conditional independence
The last section established a theoretical foundation for marginal independence testing, now this foundation will be expanded by adding a framework for testing conditional independence.
The statement one would like to make thus connects two random variables, X and Y , that are conditionally independent given a third variable, Z, taking values in Z, with the expected generalization
loss from predicting Y from Z and from the set {X, Z}.
In slight extension of our setting, we consider prediction functionals in [X × Z → Y], i.e., we separate
the features in two parts corresponding to X and the conditioning Z.
We generalize the main definitions from Section 3.3 to the conditional case:
Definition 3.16. A prediction functional f : X × Z → Y is called conditionally uninformed if it does
not depend on the first argument, i.e., if f (x, z) = f (y, z) for all x, y ∈ X and z ∈ Z. By notational
convention, a functional g : Z → Y is identified with the conditionally uninformed prediction functional
(x, z) 7→ g(z).
In the text, it will be implicitly understood that conditioning will happen on the second argument.
We define straightforward conditional generalizations of baselines and best predictors:
Definition 3.17. We define the following prediction functionals in [X × Z → Y :
(L)
The best conditionally uninformed prediction ωY |Z : (x, z) 7→ µL ([Y |Z = z])
(L)

The best conditional prediction ωY |X,Z : (x, z) 7→ µL ([Y |X = x, Z = z])
(L)

(L)

It is briefly checked that ωY |Z and ωY |X,Z have the properties their names imply:
Lemma 3.18. Let L be a (strictly) convex loss. The following holds:
(L)

(i) εL (ωY |Z ) = min{εL (f ) : f is a conditionally uninformed prediction functional}
18

(L)

(ii) εL (ωY |X,Z ) = min{εL (f ) : f ∈ [X × Z → Y]}
Proof. (i) Conditioning on the event Z = z, Lemma 3.7 for the unconditional case implies the equality
conditional for the event. Since z was arbitrary, it holds without the conditioning.
(ii) This is directly implied by Proposition 3.5, by substituting the joint random variable (X, Z) for
the X in Proposition 3.5.
With these definitions, the conditional variant of Theorem 3 can be stated:
Theorem 4. As in our setting, consider three random variables X, Y, Z, taking values in X, Y, Z,
where Y is finitely supported pmf (classification) or where Y ⊆ Rq (regression). The following are
equivalent:
(i) X and Y are statistically dependent conditional on Z.
(L)

(L)

(ii) εL (ωY |(X,Z) )
εL (ωY |Z ) for L the log-loss/Brier-loss (classification), resp.
loss L (regression)

for some convex
(L)

(iii) there exists a prediction functional f : X × Z → Y such that εL (f )
εL (ωY |Z ) for L the
log-loss/Brier-loss (classification), resp. for some convex loss L (regression)
(iv) there exists a prediction functional f : X × Z → Y such that for all conditionally uninformed
prediction functionals g : (X×)Z → Y, one has εL (f )
εL (g) for L the log-loss/Brier-loss
(classification), resp. for some convex loss L (regression)
The set of losses in which existence is required in (ii) may be restricted to a set of losses which is
faithful for the unconditional setting (such as: quantile losses for Y = R as per Theorem 2), as in
Section 3.4, without changing the fact that the stated equivalence is correct.
Proof. The proof is an analogue of that of Corollary 3.15. It may be obtained from following the
whole proof through while in addition conditioning on Z = z, and noting that (i) holds if and only if
the conditionals X|Z = z and Y |Z = z are statistically dependent for some z.
Our conditional independence test is based on statement (iv) in Theorem 4. Unlike in the uncondi(L)
tional case, there is no easy way to estimate ωY |Z directly. Thus, in the algorithm in Section 5.1,
the statement will be replaced by a slightly weaker statement, where automatic model selection will
(L)
(L)
determine an f as well as a g, without guarantees that f = ωY |X,Z or that g = ωY |Z are optimal.
This is in line with the usual approach to and caveats of supervised learning - if there were a universal
(L)
way to estimate ωY |Z , that would amount to a universally perfect supervised learning strategy.

3.6. Testing prediction error against baseline
Theorems 3 or 4 in the previous sections establish a new basis for (marginal or conditional) independence testing. Namely, the theorems relate testing of (marginal or conditional) dependence between
X and Y to testing predictability of Y from X. Namely, Theorems 3 or 4 state the following: If there
exists a significantly better-than-(conditionally-)uninformed prediction strategy f , i.e,
εL (f )

εL (g),
(L)

(1)
(L)

where g is a suitable baseline (an approximation of ωY or ωY |Z ), then we may conclude by Theorem 3
or 4 that X and Y are not (marginally or conditionally) independent.
Thus the problem of independence testing is reduced to the problem of testing whether there exists
a prediction functional f which outperforms the baseline g as measured by some convex loss function
L.
We make a few remarks about the strategy and logic of our proposed significance test.

19

• Neither the proposed functional f nor the baseline g are a-priori known in the usual practical
setting, and neither is an L which may make a difference apparent. However, both f and g
may be seen to approximate a best prediction functional which is unknown, thus as instances
of supervised learning. Hence, some choice has to be made, in absence of the “perfect learning
strategy”. A bad choice will contribute to a loss of power only if type I error control is achieved.
(L)

• The baseline g approximates ωY in the marginal case. For frequently used L, the baseline g
predicts constants which are the mean or the median or other well-studied statistics of a sample
from Y , hence g may have a beneficial asymptotic.
• Significance tests to compare prediction strategies are studied in [29], which proposes amongst
others a paired sample test of prediction residuals.
The values εL (f ), εL (g) may be estimated by respective empirical estimates, by the usual estimator
M
1 X
Li (f )
εbL (f ) :=
M i=1

where Li (f ) = L(f (Xi∗ ), Yi∗ ),

and similar εbL (g) for g. Since the test data (Xi∗ , Yi∗ ) are independent of each other and of f, g, by the
central limit theorem, one has
√
d
M (b
εL (f ) − ε(f )) −
→ N(0, Var[L(f (X), Y )]), M → ∞,
conditional on f and g being fixed. That is, the empirical mean of the loss residuals for prediction
strategy f and loss L is asymptotically normal with mean ε(f ) and variance Var[L(f (X), Y )])/N .
Instead of directly estimating εL (f ) and εL (g) with confidence intervals and then comparing, one
notes (as also in [29]) that the samples of loss residuals Li (f ) and Li (g) are inherently paired, which
eventually leads to a more powerful testing procedure.
We hence consider the difference in the i-th loss residual, Ri := Li (g) − Li (f ). The Ri are still
i.i.d. (conditional on f, g) and also have a normal asymptotic, usually with smaller variance than
either of f, g in empiry. An effect size of the difference is obtained with normal confidence intervals,
and one may conduct a one-sided paired two-sample test for the null
H0 : E(ε(g) − ε(f )) ≤ 0 against HA : E(ε(g) − ε(f ))

0.

Note that the test is one-sided since we test for f to outperform the baseline g.
To assess this null-hypothesis, two simple tests, one parametric and one non-parametric, are implemented.
• T-test for paired samples [29]: A parametric test assuming that the sample mean of loss
residual differences Ri is normally distributed with mean 0 and standard deviation σ. Under
the assumptions of the test and the null hypothesis, the normalized empirical mean t := m
d
u/b
σ,
where m
d
u and σ
b are the sample mean and standard deviation of the Ri , respectively, follows a
t-distribution with M degrees of freedom, where we have used the fact that under H0 it holds
that E[R] = 0.
• Wilcoxon signed-rank test: If the sample mean of loss residual differences Ri is not normally
distributed, an increase in power can be achieved by testing if the rank differences are symmetric
around the median. Since the Ri in question are differences in loss residuals, there is no ex-ante
reason to assume normality, in particular for a “good” non-baseline method one may expect
that most are positive, hence the normality assumption may be too strong. Hence the nonparametric Wilcoxon class test should be a better general approach to compare loss residuals
(see pages 38-52 of [11] for the Wilcoxon tests).
20

One subtlety to notice is that instead of testing for the alternative hypothesis “there exists f such
(L)
that εL (f ) ≤ εL (ωY )” which certifies for independence, by the above strategy we are testing for
the alternative “there exists f such that for all/the best g it holds that εL (f ) ≤ εL (g)” where g is a
(L)
fixed estimate of ωY . However, g is dependent on the training data and itself uncertain. Thus, for a
(L)
conservative estimate of the significance level, uncertainty in g estimating ωY needs to be taken into
account.
(L)
In the unconditional case, ωY will be a constant predictor of an elicited statistic of the training
labels, such as the mean, median or one minus the majority class frequency, with known central limit
theorems that may be used to adapt the Student or Wilcoxon type tests (e.g., by explicitly increasing
the variance in the t-test).
In either case, one may also bootstrap the distribution of g and make use of it as follows: Since the
training set is independent of the test set, the function g, is as a random variable which is a function
of the training data, also independent of the test data. Thus instead of using a fixed g in prediction,
one can ensure a conservative test by using a bootstrap sample of pseudo-inputs for the predictions,
a different sample of g per test data point that is predicted.
Since the exact form of the correction appears an unresolved sub-problem of the predictive model
validation and model selection process in general, and since bootstrapping of g is computationally
costly, we have implemented the “fixed g” variant (potentially liberal) in the package which can be
easily adapted to potential corrections based on explicit asymptotics.
We will present an algorithmic implementation in Section 5.1 as an algorithmic conditional testing
routine, which is later applied to graphical model structure learning, in Section 5.2

21

4. Graphical Models
This section will briefly review graphical models and graphical model structure learning, which we
will later address using the conditional independence test outlined in Section 3.

4.1. Informal definition of graphical models
A probabilistic graphical model is a graph-based description of properties of a joint probability distribution (possibly but not necessarily a full specification thereof). It usually consists of a directed
or undirected graph, together with a fixed convention how to translate the graph structure into a
full probability distribution, or more generally, into conditional independence statements about a
distribution of not further specified type.
Graphs
A graph G is an ordered pair (V, E) where V ⊆ E × E. The set V are interpreted as vertices (nodes)
of the graph, and the elements in E are interpreted as edges (links) of the graph (ordered for directed
graphs, unordered for undirected graphs). G is usually visualized by drawing all vertices V and edges
between all pairs of vertices in E. In graphical models, the vertex set is identified with a collection of
(possibly associated) random variables X = [X1 , ..., Xn ], and the edges encode some independence or
conditional independence information about the components of X. Two popular choices are discussed
in Section 4.2 below.
Probabilistic graphical models
Graphical model theory is usually concenred two main tasks, structure learning (of the graph) and
inference (on a given graph). Inference is concerned with estimating parameters or statistics of the
parametric distribution assuming the independence structure prescribed by a fixed graph, from a finite
sample drawn from X. Structure learning is concerned with inferring the graph from such a sample,
thus inferring the encoded independence relations. For an extended introduction to graphical models,
the reader might refer to [26] or [2].
Graphical model structure learning is usually considered the more difficult of the two tasks, due to the
combinatorial explosion of possible graphs. Manual approaches involve an expert encoding domain
knowledge in presence or absence of certain edges; automated approaches usually conduct inference
based on parametric distributional assumptions combined with selection heuristics [26].

4.2. Types of graphical models
We review two of the most frequently used types of graphical models.
4.2.1. Bayesian Networks
A Bayesian network is a graphical model which states conditional independence assumptions without making parametric distributional assumptions. The conditional independence assumptions are
encoded in a directed acyclic graph over a set of random variables. Acyclicity of the graph implies
that each node Xi has a (potentially empty) set of descendants,
Definition 4.1. A node Xj is a descendant of Xi in the graph G if there is a directed path from Xi
to Xj , where a path is any connection of links that lead from Xi and Xj .
That is, if there is a path following the arrows in G, going from Xi to Xj . Since the graph is acyclic,
no cycles exist, and following the arrows in the graph, the same variable can never be reached twice.
Examples of such a networks are shown in Figure 1.

22

Difficulty

Intelligence

Difficulty

Intelligence

GPA

Grade

Grade

Grade

Letter

Letter

Letter

(a)

(b)

(c)

Figure 1: Example of expert based graphical model structure learning (backtracking), adapted
from [26].
(a) The quality of an academic reference letter is determined by the grade
(b) The students intelligence and course difficulty determine the grade.
(c) Knowing a students GPA gives additional information about the state of intelligence
Bayesian Network graphs arise in a natural way when considering the factorization properties of a
probability distribution. Assume we are interested in a probability distribution P over the Dif f iculty
of a course, a students Intelligence, the Grade a student achieved in a course, and the quality of
a reference Letter received from the tutor of the course, P (Dif f iculty, Intelligence, Grade, Letter).
Further assume the following is true
|=

• The quality of the Letter is solely determined by the Grade a student received in the course.
That is, Letter {Dif f iculty, Intelligence}|Grade
• The Grade of a student depends on the Dif f iculty of the course and his Intelligence
• Dif f iculty and Intelligence are not causally influenced by any other variable in the graph
This gives a natural way to order the variables for factorization. Dif f iculty and Intelligence are
so-called root nodes, hence their order is irrelevant, however Grade depends on both of them and
Letter depends on Grade, giving the ordering: {Letter, Grade, Dif f iculty, Intelligence}, which will
now be denoted as {L, G, D, I}. An ordering {X1 , ..., Xn } implies that for ∀i, j ∈ {1, ..., n} and i < j,
Xj can not be a descendant of Xi . The distribution is then factorized, according to the ordering,
P (L, G, D, I) = P (L|G, D, I)P (G, |D, I)P (D|I)P (I),
a standard way to factorize distributions. Working in the independence statements gives
P (L, G, D, I) = P (L|G)P (G, |D, I)P (D)P (I),
|=

|=

since L {D, I}|G and D I. Returning to Figure 1 (b) shows that this factorized distribution exactly
matches the arrows in the graph, which start at the parents (the conditioning variables) and lead to
the children (the variable that the conditional distribution is over).
To understand the independence statements encoded in a Bayesian Network, one needs to first distinguish between active and passive trails [26, p. 71]. Let a trail be any path between Xi and Xj on G
and a v-structure a structure such that Xi−1 → Xi ← Xi+1 , where Xi is a descendant of both Xi−1
and Xi+1 .
23

A

C

E

B

D

F

G

Figure 2: Markov Network
Definition 4.2. A trail between Xi and Xj is active given a conditioning set Z, if for every v-structure
Xi−1 → Xi ← Xi+1 along the trail, Xi or one if it’s descendants is in Z and no other variables on
the trail are in Z.
In Figure 1 (c) that means, the trail from Letter to GPA is active only if conditioned on either the
empty set or Difficulty, and the trail between Intelligence and Difficulty is only active if Grade or
Letter is in the conditioning set.
|=

Definition 4.3. [2] If there is no active trail in the graph G between the nodes X and Y given Z,
then X Y |Z in any distribution consistent with the graph G
4.2.2. Markov Networks
A Markov Network, or Markov Random Field, is a graphical model which encodes conditional independence statements in an undirected graph over a set of random variables. While the Bayesian
Network is defined in terms of probability distributions, the Markov Network is usually specified in
terms of factor products which are unnormalized probability distributions. The scopes of these factors
determine the complexity of the system, if the scope of a factor covers all the variables in the graph,
the graph is unconstrained, whereas any constraint to smaller scopes decreases its complexity.
Figure 2 shows the Markov Network for the (normalized) factor product
P (a, b, c, d, e, f, g) =

1
φ1 (a, b, c)φ2 (b, c, d)φ3 (c, e)φ4 (e, f )φ4 (f, g)φ5 (e, g),
Z

with the normalizing constant (partition function) Z. Small letters denote realizations of the capitallettered nodes. If P is a probability mass function,
X
Z=
φ1 (a, b, c)φ2 (b, c, d)φ3 (c, e)φ4 (e, f )φ4 (f, g)φ5 (e, g),
a,b,c,d,e,f,g

|=

|=

where φ(s) is such that φ : S → R, ∀s ∈ S.
If P is a probability density function over continuous variables, Z would be attained by integrating
over the support of the variables in the scopes of P .
Since the links are undirected, independence statements arising from Markov Networks are symmetric.
A path between Xi and Xj , i 6= j is called active given Z, if no node on the path is in Z. To determine
whether Xi and Xj are conditionally independent given Z in all distributions consistent with G, one
again considers all paths going from Xi to Xj . It holds that Xi Xj |Z if there is no path between Xi
and Xj that is active given Z.
So to attain a Markov Network that is consistent with a probability distribution, one has to consider
the pairwise conditional independence properties of the variables that the nodes in the graph represent.
In general, if in a Markov Network there is no link between a variable Xi and Xj , then Xi Xj |X \
{Xi , Xj } in any distribution P over X consistent with the graph. These are called the pairwise
Markov-independencies of the graph [15, ch. 17.2].
The independencies that can be encoded using Bayesian and Markov Networks differ. Bayesian
24

X

Y

X

Y

Figure 3: Two DAGs equivalent with respect to independence statements
Networks are natural to express (hypothesized) causal relationships, distributions where an ordering
can be attained and thus have nice factorization properties, while Markov Networks are natural
for expressing a set of pairwise Markov independencies. Additionally, Bayesian Networks can be
transformed into an undirected graph by a (potentially forgetful) process called ”moralization” [26,
p. 135].

4.3. Graphical model structure learning
There are two dominant approaches to structure learning, independence testing based and score-based
methods, however they both suffer to a varying extent from the same underlying problem: combinatorial explosion. Given a probability distribution P (with associated graph G) over a multivariate
random variable X = [X1 , ..., Xn ], for undirected networks, between any pair of variables Xi and Xj
potential
there can be a link or no link in G. Since the links are undirected, there are thus p(p−1)
2
edges in the graph, and an exponential number of possible graphs. If one wants to test if Xi and Xj
are independent given Z, Z ⊂ X, one needs to test if any path between the two is active, leading to a
potentially very large number of tests. A graph’s complexity can be decreased by upper bounding the
number of edges for each node, however, it was shown that, for directed graphs where each node can
have d parents, for d > 1, the problem a finding an optimal graph is NP-hard [26, p. 811]. Generally,
the space of graphs is reduced by making structural (what type of relationships can be captured in the
graph) and parametric (constraints on the distribution) assumptions. An additional problem is posed
when there exists more than one optimum to the structure learning method, as a result of the fact
that different graphs can be equivalent (and thus not identifiable) in the respective search space. An
example of this are the two graphs shown in Figure 3, which arises from the fact that if {X, Y } ∼ P ,
P can be decomposed either into P (X|Y )P (Y ) or P (Y |X)P (X), which are equivalent. There are
different ways to approach these problems, some of which are outlined below.
4.3.1. Score-based methods
Score-based methods attempt to find a (usually local) optimum of a score function in the space of
possible graphs. Examples of this can be the Kullback-Leibler divergence scoring function for directed
graphical models, which can be used to find the optimal Chow-Liu tree graph [2, p.219]. While this
is a convex algorithm, finding the global optimum, it does so under the heavy constraint that each
node only has one parent at maximum (resulting in a tree-structured graph). For undirected graphs
over discrete variables, gradient descent can be used directly on the log likelihood to find the local
optimum [2, p. 221], however each gradient evaluation requires calculation of the partition function
Z (by summing over all discrete states), which makes this algorithm very expensive. Many of the
algorithms for score-based methods are designed for discrete variables, and not applicable when the
variables are continuous. Performance is usually evaluated by calculating the likelihood score on a
hold-out sample. One other area not mentioned above is Bayesian model averaging (e.g. [14]), which
seeks to improve the above methods by averaging over different model structures. It can be viewed
as an ensembling method for the score-based methods mentioned above. State-of-the-art score based
methods are oftentimes very expensive or make very strong assumptions on the underlying distribution
P , such as tree-structure of the graph [26, ch. 18].

25

4.3.2. Independence testing based methods
Unlike the score-based models, independence testing based models conduct independence tests between
variables or sets of variables locally and then aggregate the information to produce the graphical model.
These sets of independence tests are then used to infer the structure of the graph in a straightforward
manner, based on the properties of the respective type of graphical model. An example of this is the
PC-Algorithm [2, e.g. p. 214] which is used in attempts at causal discovery. Another algorithm for
recovering the undirected skeleton is given by Algorithm 3.3 of [26].
There are probably two main reasons that this latter independence testing based method is not used
for graphical modelling:
(a) is that it relies on a conditional independence test where the conditioning is on multiple variables,
i.e., all except two. As outlined in Section 2, this is a largely unsolved problem expect when
there is at most one conditioning variable, i.e., if there are three or less variables in total.
(b) It is hard to regulate the bias-variance trade-off, unlike for score-based methods where this may
be achieved through the use of constraints such as an upper bound on the number of parents.
However, these two issues may be addressed through our novel predictive conditional independence
testing framework:
(a) The predictive conditional independence test described in the subsequent Section 5.1, based on
the theoretical insights in Section 3, allows for efficient conditional independence testing with
variables of any dimension, and thus provides a predictive approach to learning graphical models.
(b) The link to the supervised learning workflow allows for direct adoption of well-known strategies to
trade-off bias and variance and error estimation from supervised learning, such as regularization
and cross-validatied tuning, to the graphical modelling workflow..
In addition, our suggested algorithms have further desirable features:
• The intimate link to supervised learning also allows for a controlled trade-off between power of
the algorithm and time complexity. This trade-off is readily available when using the conditional
independence testing routine described in Section 5.1, since the user can choose estimation
methods of varying power/complexity in order to influence the total run time needed by the
algorithm.
• Section 5.1.1 will also introduce a false-discovery-rate control routine to provide a tool that lets a
user control the proportion of false-positives (erroneously found edges in the estimated graphical
model structure), regardless of the size of the graph.
Note on causal inference
Graphical models, and specifically Bayesian Networks, are a natural way to express hypothesized
causality, since the arrows seem to express causation. However, when actually learning graphical
models from data, causality may only be inferred by making strong (and often incorrect) assumptions
on the underlying distribution, or by collecting data in a manner that allows for causal interpretation,
namely in proper experimental set-ups (e.g., following Pearl’s do-calculus, or standard experimental
study design). As generally for graphical model structure learning, all algorithms outlined in this paper
are not to be interpreted as expressing causality, but merely as producing a collection of statements
about association which certify for causality only in combination with the right experimental set-up.

26

5. Predictive inference algorithms
This section will first introduce the proposed predictive conditional independence testing routine
(PCIT), which is based on the methodology outlined in Section 3, and important subroutines related
to performance and error control of the test. After, an algorithm to leverage the independence test into
a graphical model structure learning routine that addresses the issues outlined in 4.3.2 is presented.

5.1. Predictive conditional independence testing
Algorithm 1 implements the results from Section 3 to test if a set of variables Y is independent of
another set of variables X, given a conditioning set Z (optional, if not provided the algorithm tests
for marginal independence). It will later be used as a subroutine in the graphical model structure
learning, but can also be used in it’s own right for tasks such as determining if a subset X would add
additional information when trying to predict Y from Z.
Algorithm 1 Predictive conditional independent test (PCIT)
Split data into training and test set
for all variables y ∈ Y do
on training data:
4:
find optimal functional f for predicting y from Z
5:
find optimal functional g for predicting y from {X,Z}
1:
2:
3:

6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

on test data:
calculate and store p-value for test that generalization loss of g is lower than f
end for
if symmetric test needed then
exchange X and Y, repeat above process
end if
p values adjusted ← Apply FDR control to array of all calculated p-values
return p values adjusted

When the test is used as a marginal independence test, the optimal prediction functional for f is
the functional elicited by the loss function, that is, the mean of y for continuous outputs, and the
class probabilities of y for the discrete case (Appendix ??). The link to supervised learning further
allows/forces the user to distinguish between two cases. Independence statements are symmetric, if
X is independent of Y , then Y is independent of X, in both the marginal and conditional setting.
The same cannot be said in supervised learning, where adding X to predicting Y from Z might result
in a significant improvement, but adding Y does not significantly improve the prediction of X, given
Z (as can be seen in the asymmetry of OLS). So if a user is interested in a one-sided statement, the
algorithm can be run for a single direction, for example to evaluate if a new set of variables improves
a prediction method, and is thus worth collecting for the whole population. If one is interested in
making a general statement about independence, the test is “symmetrized” be exchanging X and Y ,
and thus testing in both directions, and FDR control be applied to the union of the p-values.
It is important to distinguish between the two types of tests and null-hypotheses in this test. On one
hand (in the symmetric case) for each variable in X and Y , it will be assessed, if adding X to the
prediction of y ∈ Y from Z results in an improvement (and vice versa). The null-hypothesis of this
“prediction-null” is that no improvement is achieved. After all these p-values are collected, a multiple
testing adjustment is applied (Section 5.1.1), after which the original null-hypothesis, that X and
27

Algorithm 2 The Benjamini-Hochberg-Yekuteli Procedure
j
Require: Set {p(i) }m
i=1 s.t. pj : p-value for observing Xj under H0
1: Sort the p-values in ascending order, p(1) ≤ ... ≤ p(m)
Pm
2: Let q = α/( i=1 1/i) for chosen confidence level α
i
3: Find the k s.t. k = max(i : p(i) ≤ m
q)
j
4: Reject H0 for j ∈ 1, ..., k

Y are conditionally independent, is assessed. We reject this “independence-null”, if any one of the
“prediction-nulls” can be rejected after adjustment. The p-value of the “independence-null” hence
reduces to the lowest p-value in all the “prediction-nulls”. As such, the null of the independence test
is that all the “prediction-nulls” are true. False discovery rate control, the chosen multiple testing
adjustment, is appropriate since it controls the family-wise error rate (FWER), the probability of
making at least one type 1 error, if all the null-hypotheses are true [4].
5.1.1. False-discovery rate control
To account for the multiple testing problem in Algorithm 1, the Benjamini - Hochberg - Yekutieli
procedure for false-discovery rate (FDR) control is implemented [5]. In their paper, they state that
traditional multiple testing adjustments, such as the Bonferroni method, focus on preserving the
FWER. That is, they aim to preserve the probability of making any one type 1 error at the chosen
confidence level. As a result, tests are usually very conservative, since in many multiple-testing scenarios the p-values are not independently distributed (under the null), and thus the power of these
tests can be significantly reduced.
In their 2001 paper, they propose to control the false discovery rate instead, “[..] the expected proportion of erroneous rejections among all rejections”, as an alternative to the FWER. The FDR allows
for more errors (in absolute terms) when many null-hypothesis are rejected, and less errors when few
null-hypotheses are rejected.
Algorithm 2 shows the procedure outlined in [5]. They showed, that this procedure always controls
the FDR at a level that is proportional to the fraction of true hypotheses. As hinted at before, while
this algorithm controls the false-discovery rate, in the special case where all null-hypotheses in the
multiple testing task are assumed to be true, it controls the FWER, which then coincides with the
FDR. This is especially useful, since both scenarios occur in the graphical model structure learning
routine described in Section 5.2.
For the choice of optimal false discovery-rate for an investigation, even more so than in the classical
choice of appropriate type 1 error, there is no simple answer for which rate might serve as a good
default, and it is highly application dependent. If the goal of the procedure is to gain some insight
into the data (without dire consequences for a false-discovery), a user might choose a FDR as high
as 20%, meaning that, in the case of graphical model structure learning, one in five of the discovered
links is wrong on average, which might still be justifiable when trying to gain insight into clustering
properties of the data. This paper will still set the default rate to 5%, but deviate willingly from this
standard whenever deemed appropriate, as should any user of the test.
5.1.2. Improving the prediction functionals
In practice, when assessing the individual “prediction-nulls” in Algorithm 1, the power of the test
(when holding the loss function constant) depends on the capability of the the prediction method to
find a suitable functional g that outperforms the baseline f . That means, a general implementation of
the independence test needs to include a routine to automatically determine good prediction functionals for g and f . The implementations in the pcit package presented in Section 6 support two methods
for the automatic selection of an optimal prediction functional. Both methods ensemble over a set
28

Stage 1

Regression

Classification

ElasticNetCV

BernoulliNB

GradientBoostingRegressor

MultinomialNB

RandomForestRegressor

GaussianNB

SVR

SGDClassifier
RandomForestClassifier
SVC

Stage 2

LinearRegression

LogisticRegression

Table 1: Prediction functionals used for Stacking/Multiplexing
of estimators, which are shown in Table 1. The prediction functionals refer to the estimator names
in sklearn1 . Some are constrained to specific cases (e.g. BernouilliNB, which only applies when the
classification problem is binary).
Stacking Stacking refers to a two-stage model, where in the first stage, a set of prediction function
is fit on a training set. In the second stage, another prediction function is fit on the outputs of the
first stage. If the prediction function in the second stage is a linear regression, this can be viewed as a
simple weighted average of the outputs in the first stage. In theory, stacking allows the user to fit one
single stacking predictor instead of having to compare many potential prediction functions based on
the model diagnostics, as in the second stage, better methods get more weight (in the expectation).
As an additional effect, improvement in the prediction accuracy through ensembling of predictors can
take place (see e.g. Section 19.5 of [1]). The used stacking regressor and classifier can be found in the
Python package Mlxtend [32].
Multiplexing When multiplexing over estimators, the training set is first split into a training and
validation set, a common procedure to find optimal hyperparameters. After, the predictors are fit
individually on the training set, and each predictors expected generalization loss is estimated on the
validation set. One then proceeds by choosing the predictor with the lowest estimate for the empirical
generalization loss, and refits it using the whole training data (including the former validation set),
and then uses the fitted estimator for prediction.
5.1.3. Supervised learning for independence testing
Algorithm 1 shows how to leverage supervised learning methodology into a conditional independence
test. This has a major advantage over the methods outlined in Section 2, as the supervised prediction
workflow is of utmost interest to many areas of science and business, and, as a result, a lot of resources
are going into development and improvement of the existing methods. By making a link between
predictive modelling and independence testing, the power of independence testing will grow in the
continuously increasing power of the predictive modelling algorithms.

5.2. Predictive structure learning of undirected graphical models
This section outlines a routine to learn the vertices in a directed graph (the skeleton) for a data
set by conducting a range of conditional independence tests with the null hypothesis of conditional
independence. Section 4.2 outlines the conditional independence statements of a Markov network. In
a directed graph, if variables Xi and Xj have no direct edge between them, they are conditionally
independent given all other variables in the network.
1 Details

can be found here http://scikit-learn.org/stable/modules/classes.html

29

Algorithm 3 Undirected graph structure estimation
for any combination Xi , Xj s.t. i 6= j do
X− ← X \ {Xi , Xj }
p vali,j ← p-value for test Xi Xj |X−
end for
p val adj ← Apply FDR control on p val matrix
return p val adj
|=

1:
2:
3:
4:
5:
6:

Algorithm 3 describes the skeleton learning algorithm for an input set X = [X1 , ..., Xn ], by considering
all possible pairs of variables in the data set, and testing if they are conditionally independent given
all other variables. The output p val adj is a symmetric matrix with entries i,j being the p-value
for the hypothesis that in the underlying distribution, Xi and Xj are independent, given all other
variables in X, and hence in the graph G describing it, there is no link between the vertices for Xi and
Xj . Ultimately, links should be drawn where the adjusted p-values are below a predefined threshold.
There are O(n2 ) function evaluations in the for-loop, where n is the number of variables. Section 7
will provide experimental performance statistics for the algorithm and showcase applications on real
data sets.

30

6. pcit package
6.1. Overview
The Python2 package implementing the findings and algorithms of this paper can be found on https:
//github.com/SamBurkart/pcit and is distributed under the name pcit in the Python Package Index.
This section will first provide an overview of the structure and important functions of the package. As
this implementation can be thought of as a wrapper for scikit-learn (sklearn) estimators, this section
will then describe the sklearn interface and how the package interacts with the sklearn estimators.
Lastly, simple application-examples are given.
6.1.1. Use cases
The package has two main purposes, independence testing and structure learning. While univariate
unconditional independence testing is possible, it is not expected to outperform current methodology
in the simple tasks, but the difficult ones, such as
Multivariate independence tests, such as for hedging purposes in a multivariate test to determine
which financial products are independent from the ones already in the portfolio.
Conditional independence tests, such as assessing if a data set adds information to existing data
in a prediction task.
Graphical model structure learning, such as for finding clusters in the data as part of an exploratory data analysis or thoroughly investigating associations in the data (both in section 7.2.2).
For these tasks, the PCIT serves as a readily available tool that works without the need for manual
choices or hyperparameter tuning, and scales well in the dimensionality of the data.
6.1.2. Dependencies
The package has the following dependencies:
Scipy [25], for the calculation of p-values
Sklearn [31], for its estimators (predictors)
Mlxtend [32], for the implementation of stacking

6.2. API description
The package introduces three main routines, one for automated prediction (MetaEstimator), conditional independence testing (PCIT), and undirected graph structure learning (find neighbours). The
following section gives on overview, the function signatures can be found in Appendix 6.3.
MetaEstimator
The MetaEstimator class provides a user with a type-independent predictor that automates model selection for given training data, by automatically determining appropriate loss functions and prediction
functionals. It is initialized for sensible defaults, which should generally lead to good results, but can
be changed for specific tasks (such as the use of more complex base estimators for more powerful, but
also more computationally demanding, predictions). The ensembling methods used for the estimator
are described in Section 5.1.2. For regression, the square loss is used for training and calculation of
the residuals, for classification, the logistic loss serves as the loss function.
2 https://www.python.org/

31

PCIT
PCIT implements the conditional independence test in Algorithm 1 to test if two samples stem from
conditionally (or marginally) independent random variables. The MetaEstimator class is used as a
prediction functional, and hence the user can trade off between computational complexity and power
by adjusting the chosen MetaEstimator. That is, if speed is important, the used MetaEstimator should
be a combination of base estimators and ensembling method that is quick in execution, whereas if
computational resources are vast and a more powerful test is needed, the used base-estimators should
be highly tuned.

find neighbours
find neighbours implements Algorithm 5.2 to learn the undirected skeleton for an input data set X,
using the PCIT.

6.3. Function signatures
PCIT
Type: function

Inputs:

Outputs:

Name

Description (type)

Default

x

Input data set 1 ([n x p] numpy array)

y

Input data set 2 ([n x q] numpy array)

z

Conditioning set ([n x r] numpy array)

None (empty set)

estimator

Estimator object to use for test (MetaEstimator)

MetaEstimator()

parametric

Parametric or nonparametric test (bool), section 3.6

False

confidence

Confidence level for test (float [0,1])

0.05

symmetric

Conducts symmetric test (bool), section 5.1

True

p values adj

p-values for ”prediction nulls” of each y ∈ Y (list)

independent

tuple, first value shows if ”independence-null” is rejected (bool)
second value is p-value of ”independence-null” (float [0,1])

loss statistics

RMSE difference for baseline f and altern. g,
loss residuals with standard deviation, for each y ∈ Y .
Only applicable if Y continuous

Note: The variance of the difference is estimated by assuming 0 covariance between the residuals
of baseline f and alternative g, which generally leads to more conservative confidence intervals for
error residuals (due to the irreducible error, prediction residuals for different methods are generally
positively correlated).

32

MetaEstimator
Type: class
Methods
init

Name

Description (type)

Default

Inputs:

method

Ensembling method (’stacking’, ’multiplexing’ or None)

’stacking’

estimators

Estimators to ensemble over (2-tuple of lists of sklearn

None (default estim,

estimators [regression estim], [classification estim])

section 5.1.2)

method type

Task definition (’regr’, ’classif’, None)

None (auto selection)

cutoff

Cutoff for automatic selection of method type (integer)

10

Inputs:

y

Dependent variable ([n x 1] numpy array)

Outputs:

estimators

Appropriate set of estimators (list)

x

Independent variables ([n x p] numpy array)

y

Dependent variable ([n x 1] numpy array)

fitted

Fitted estimator (MetaEstimator)

x

Independent variables ([n x p] numpy array)

y

Dependent variable ([n x 1] numpy array)

fitted

Fitted uninformed baseline estimator (MetaEstimator)

get estim

fit
Inputs:

Outputs:
fit baseline
Inputs:

Outputs:
predict
Requires:

MetaEstimator has been fitted

Inputs:

x

Test set independent variables ([n x p] numpy array)

Outputs:

predictions

Predictions for test set ([n x 1] numpy array)

x train

Training set independent var. ([n x p] numpy array)

y train

Training set dependent variables ([n x 1] numpy array)

x test

Test set independent variables ([n x p] numpy array)

y test

Test set dependent variables ([n x 1] numpy array)

baseline

Should baseline be fitted (boolean)

resid

Residuals for prediction strategy ([n x 1] numpy array)

get resid
Inputs:

Outputs:

33

False

find neighbours
Type: function

Inputs:

Outputs:

Name

Description (type)

Default

X

Input data set ([n x p] numpy array)

estimator

Estimator object to use for test (MetaEstimator)

MetaEstimator()

confidence

False-discovery rate (float [0,1])

0.05

skeleton

Matrix, p-values for each indep test ([p x p] numpy array)

skeleton adj

Learnt graph after applying FDR control ([p x p] numpy array)

6.4. API design
The API is designed as a wrapper for Scikit-learn (sklearn), a package in the Python programming
language, that aims to provide a user with a consistent, easy-to-use set of tools to analyze data [31].
It is one of the most-used tools in today’s supervised learning community, which is why it is chosen as
the supervised prediction workflow to build on for the predictive conditional independence test. This
section will outline the advantages of basing the test on the sklearn package.
6.4.1. Sklearn interface
Sklearn is built around estimator objects, which implement a consistent set of methods. Estimators
provide a fit and, if applicable, a predict method, in order to be able to fit the estimator to training
data and then predict on a test set. Additionally, sklearn provides standardized approaches to model
selection and hyperparameter-tuning as well as data transformation and ensembling methods (see
[8] for a more thorough discussion). The methods can easily be combined to create more powerful
estimators. Defaults are chosen sensibly so that in most cases, an initial fit of a method to data
requires little manual parameter specification from the user’s side. While the highly automated and
simplified approach of sklearn lowers the bar of entry when aiming to generate knowledge from data,
it also comes with a downside. For most statistical applications that exceed fitting and predicting
from a data set, such as inference on the parameters and hypotheses about prediction accuracies,
the relevant subroutines are missing from the API. Relevant statistics can however be attained by
interfacing it with other packages (such as SciPy).
6.4.2. Wrapper for Sklearn estimators
As we saw in section 6.2, the conditional independence test described in Algorithm 1 uses the newly
defined MetaEstimator class to automate determining the optimal prediction functional for a given
task. It does so, by ensembling over a set of base estimators from sklearn. These are either chosen
to be the sensible defaults described in Table 1, or can be passed by the user as a tuple of lists
of sklearn base estimators. This is required since regression and classification tasks rely on vastly
different prediction functionals, and thus need to be specified separately. As a general rule, the passed
regressors need to have a fit and a predict method, whereas the classifiers need to have a fit and
a predict proba method. Requirements might be more stringent for certain types of data or certain
estimators, however specifying an unsuitable estimator class will result in an error message as specified
by the respective class, allowing the user to either remove the unsuitable class or proceed with the
sensible defaults. As mentioned before, this gives a user a flexible tool to trade off between power
and computational complexity. If in need of a fast method, one can use an algorithm that runs in
linear time, such as stochastic gradient descent linear regression, whereas if a test with high power is
needed, one can pass hyper-tuned estimators to the function, that take longer to run but generalizes
better for prediction on unseen data.
34

6.5. Examples
This section will provide some simple examples of how the code is used. For the following it is assumed
that data sets X, Y and Z, all of the size [number of samples × number of dimensions], are loaded as
numpy arrays, and have matching numbers of samples (sample indices in X, Y and Z correspond to
the same sample). After installing the pcit package, import the relevant objects:

Testing if X

Y |Z, using the default values:

I n d e p e n d e n c e T e s t . PCIT(X, Y, z = Z )
Testing if X
mators:

|=

1

|=

1 from p c i t import MetaEstimator , S t r u c t u r e E s t i m a t i o n , I n d e p e n d e n c e T e s t

Y |Z, with a custom MetaEstimator, multiplexing over a manually chosen set of esti-

1 from s k l e a r n . l i n e a r m o d e l import RidgeCV , LassoCV ,
2
SGDClassifier , LogisticRegression
3
4 r e g r e s s o r s = [ RidgeCV ( ) , LassoCV ( ) ]
5 c l a s s i f i e r s = [ SGDClassifier () , LogisticRegression ( ) ]
6
7 c u s t o m e s t i m = MetaEstimator . MetaEstimator ( method = ’ m u l t i p l e x i n g ’ ,
8
estimators = ( regressors , c l a s s i f i e r s ))
9
10 I n d e p e n d e n c e T e s t . PCIT(X, Y, z = Z ,
11
estimator = custom estim )

Learning the undirected skeleton of X:
1

S t r u c t u r e E s t i m a t i o n . f i n d n e i g h b o u r s (X)
Concrete outputs will be shown in section 7 and can be found on GitHub.

35

7. Experiments
This section will first evaluate the performance of the proposed algorithms, and then provide some
examples of applications on real world data sets. All performance tests can be found on Github3 .

7.1. Performance tests
This section will report on performance tests for the algorithms derived in Section 5. First the power
of the the predictive conditional independence routine is bench-marked against current state-of-the-art
methodology, then various tests on the directed graph structure learning algorithm are conducted.
7.1.1. Performance of conditional independence test
In this section the conditional independence routine will be bench-marked against the previous research, namely the kernel based approach for conditional independence testing, which is is shown
to be more powerful than other conditional independence testing algorithms in [37] (on multivariate
Gaussian data). The used code for the kernel test is taken from GitHub4 . To conduct a test using
data that is drawn from a distribution that more closely resembles real world data, as opposed to the
synthetic (Gaussian) data commonly used for performance tests, the UCI Wine Repository data set
[27] is used as follows:
• The columns ’Alcohol’, ’Malic Acid’ and ’Magnesium’ are randomly permuted (to make them
independent) and will serve as X, Y and noise arrays respectively
• Vector Z is created by individually sampling vectors X 0 , Y 0 and noise0 of size n with replacement
from X, Y and the noise vector, and then calculating
p
Zi = log(Xi0 ) × exp(Yi0 ) + u ∗ noise0i , i ∈ {1, ..., n},
where u is the sign, uniformly drawn from {−1, 1}
|=

This results in a scenario, where X Y , but X 6⊥⊥ Y |Z, and the signal to noise ratio is high for
small sample sizes. The test will be conducted by increasing the sample size from 100 to 5000, and
calculating the run times and power for both approaches. For each sample size, the PCIT is run 500
times, and the KCIT is run 200 times, since the KCIT is more computationally demanding than the
PCIT. The only exception is n = 5000 for the KCIT, which is run 25 times, since the time complexity
would be too high to draw more samples. Both methods are run for their default values, without
additional manual tuning, and at a confidence level of 5%. Each time, {X 0 , Y 0 , Z} is sampled, and
then the conditional independence tests are applied and the run time is recorded. If they reject
independence at a 5% level, the round counts as a success, 1, otherwise 0. The power and run times
are then calculated by averaging over the results, and standard errors for the power are attained by
realizing that the standard error of the power for a rerun number of B is the standard error of X
B,
where X ∼ Bin(B, θ), where θ is the observed power (the sample mean).
The results are shown in Table 2. The power at the higher end of the sample sizes seems to be similar
(it is important to note that the power of 1 for the KCIT for n = 5000 was achieved on 25 resamples
only), where as in the range 500 to 1000 samples, the proposed predictive conditional independence test
shows a significantly higher power. For small n, the KCIT seems to fare significantly better, however
both approaches have a very low power, and the PCIT especially shows the power levels below the
confidence levels, which might indicate a discrepancy between true type 1 error and expected type 1
error. Important to note is the very high computational complexity of the kernel-based approach for
a data set of size 5000, with a run time of approximately 80 minutes per test, while the predictive
3 https://github.com/SamBurkart/pcit/tree/master/tests
4 https://github.com/devinjacobson/prediction/tree/master/correlations/samplecode/KCI-test

36

n
PCIT

100

Power
Time (s)

KCIT

Power

200

500

1000

2000

5000

0.020

0.046

0.332

0.672

0.832

0.970

(0.006)

(0.009)

(0.021)

(0.021)

(0.017)

(0.007)

0.32

0.38

0.49

0.624

1.31

4.79

0.050

0.085

0.185

0.325

0.8

1

(0.015)

(0.019)

(0.027)

(0.033)

(0.028)

(∗)

Time (s)

0.57

1.25

9.8

44

383

4758

Stand. difference

-1.8

-1.78

4.25

8.85

0.97

*

Table 2: Performance statistics for the newly proposed predictive conditional independence test
(PCIT) and the kernel based approach (KCIT). The values in brackets show the estimated standard errors. The last row shows the standardized difference between the power estimates, PCIT KCIT
.
conditional independence test still has a very low run time of 4.8 seconds. This is to be taken with a
grain of salt, since the tests were run in different languages (PCIT in Python, KCIT in MATLAB),
but it is apparent that PCIT scales much better than KCIT, while the both converge to a power of 1.
7.1.2. Performance of structure learning algorithm: Error rates
This test will show that the graphical model learning algorithm is capable of recovering the true graph
with high probability as the sample size increases. No comparison with alternative tests is made, as
it would lead to infeasible run times.
The data used in the performance tests is generated as follows:
1. Sample a random positive definite, sparse, symmetric precision matrix M . The size of the
entries of this matrix (relative to the diagonal) determine the signal to noise ratio, and are thus
thresholded to 0, if below a certain value.

|=

2. Invert M , M 0 = inv(M) and use M 0 as covariance matrix to sample from a multivariate Gaussian
distribution. This has the effect, that the zero-entries in M express zero-partial correlations [15,
ch. 17.3] between respective variables (and hence, lack of an edge between the two nodes
in the graph). That is, for a multivariate Gaussian random variable X, X = [X1 , ..., Xp ],
Mi,j = 0 =⇒ Xi Xj |X \ {Xi , Xj }.
3. Sample D, a data set of size n, from the multivariate normal P = N (0, M 0 ). The choice of n
will allow to evaluate the algorithms performance for an increasing signal to noise ratio.
4. The undirected graph G consistent with the probability distribution P is given by M, where
edges occur between variables Xi and Xj , if Mi,j is non-zero.
Then, Algorithm 3 will be tested by getting an estimate Ĝ of the structure of an undirected graph
G induced by the distribution P from which the data set D was drawn. The performance evaluation
will be guided by three metrics:
• False-discovery rate: The FDR is the fraction of type 1 errors (edges in Ĝ that are not in G) over
the total number of identified edges in the learned Ĝ
• Power: Number of found edges (links in true graph G that were found by the structure learning
algorithm) over the total number of links in G

37

FDR

Time (sec)

No ensembling

3.09%

30

Stacking

3.03%

450

Multiplexing

2.75%

1000

Table 3: False-discovery rates and run times for a data set of 22000 for all used methods
• Time: Run time needed to estimate Ĝ
For the test, the number of random variables is chosen to be 10. This means, that each node in the
graph has up to 9 neighbours, and the total number of possible undirected links (before inducing
sparsity) is 45. The sparsity parameter is chosen in a way that generally between 10 and 15 of those
links exist in the true graph. The size and sparsity of the graph are chosen to produce estimators of
the metrics with reasonably low variance, but are otherwise arbitrary. The sample sizes range from
approximately 400 to 20000, increasing in steps of 10%. For each sample size, 10 tests are run to
decrease the variance in the estimators. The test is conducted for conditional independence testing
using stacking and multiplexing, as well as without using any ensembling method, which, since the
data is continuous, results in the usage of Elastic Net regularized regression.
If the algorithms work as expected, the FDR is expected to be at or below 5%. High power-levels
indicate a better performance of the algorithm, with respect to this specific task. At the very least,
the power level is expected to increase in the number of samples, suggesting that asymptotically the
routine will find the correct graph.
Table 3 shows the average FDR and run times for each of the three methods. The average FDR
seems to be similar across all three methods, whereas the computational complexities differ by a
large amount. No ensembling PCIT runs very quick, about 15 times faster than stacking, which
itself only takes about half as long as multiplexing. This is the case, since multiplexing requires the
calculation of performance measures for each used estimator. Figure 4 shows the power and FDR of
the algorithm for increasing sample size. The FDR for all 3 methods seem to be slightly higher for
small sample sizes than they are for large sample sizes, but they are generally around or below the
desired 5% in the expectation (the variances are quite high, as the number of possible reruns is low
due to the computational complexity of the multiplexing method). While it might seem surprising
that stacking and multiplexing are outperformed by the no-ensembling method, one has to remember
that the ensembling is used to choose the optimal prediction functional automatically from a set of
estimators. However, the data is multivariate Gaussian, for which ridge regression is optimal in terms
of generalization performance. While stacking and multiplexing are tasked to find this property in a
large set of estimators, the estimator used in the no ensembling case is Elastic Net regularized linear
regression, a generalization of ridge regression, and hence fares better since there is less variance in
finding the optimal estimator. For all three methods, the power increases roughly logarithmically in
the sample size, implying that for a test that recovers the truth with high probability, a large data
set or a more powerful set of estimators (see Section 7.2.2) might be needed for that specific task.
However, asymptotically, all three tests seem to recover a graph that is close to the truth, in this
specific scenario, unless the power starts to plateau before reaching 1 (which there is no indication
for). Since the power for the no ensembling case is biased by the fact that it uses an optimal prediction
functional, the power curves for stacking and ensembling provide a better estimate for the performance
on an unseen data set of an unknown distribution family. As graphical model structure learning is an
inherently hard problem (due to issues such as multiple testing, combinatorial explosion, as outlined
in Sections 4), it is promising that the algorithm finds an increasing fraction of the links while keeping
the ratio of false-discoveries constant.

38

Figure 4: Power (blue lines) and FDR (red lines) for increasing sample size for all 3 methods, showing
that the performance increases as expected for all 3 methods, when the sample size is increased.
The transparent blue and red areas denote the 90% confidence intervals, the dashed line shows the
expected FDR (0.05)

39

Figure 5: Frequencies of edge occurrence in the learned structure. Green denotes edges that occur in
less than 7% (as advocated by FDR) of the models, blue for edges that occur in more than a third of
the model, and red everything in between
7.1.3. Performance of structure learning algorithm: Variance of estimated model
As a second performance metric, this section will assess the consistency of the learned structures on
resamples from a data set D. Assuming that all the observation in D are identically distributed, the
structure learning method should arrive at the same conclusions on the resamples, less some variance
in the process. The Auto MPG Data Set5 containing various continuous and discrete car performanceattributes, from the UCI machine learning repository [27], is used to conduct the test. The data set
contains 398 instances of 8 variables. For the purpose of the experiment, data sets of the same size
will be sampled with replacement from the full data set 100 times. On each resample, a graph is
learned using the stacking estimator on a 10% FDR level. Each subsample contains about two thirds
of the original data points (with some instances repeated). This is a commonly used procedure to
estimate the sampling distribution of a statistics, and which will here allow us to assess the variance
in the learned graph structure. Figure 5 shows the results. On average, there are 6 links in the learned
structure, hence the FDR advocates about 0.6 type 1 errors per learned model. The green lines are
connections that are found less times than expected by the FDR and the blue lines are connections
that are found in a large fraction of the models (over one third of the resamples). The concerning
links are the ones in between, for which the null of independence is rejected more than occasionally,
but not in a reasonably large fraction of the learned graphs. There are only 2 links in the model for
which this occurs, so overall, the variance across learned graphs seems to be reasonably small and the
learning process is hence consistent.
5 https://archive.ics.uci.edu/ml/datasets/auto+mpg

40

Figure 6: Learned graphical model structures for the Boston Housing and the Iris data set

7.2. Experiments on real data sets
In this section, the outputs of the graphical model structure learning routine are shown on a selection
of real world data sets. It will outline some possibilities of the user to trade off between power and
computational complexity.
7.2.1. Sklearn data sets: Boston Housing and Iris
Setup: One receives a data set and is interested in initial data exploration. This involves finding
related sets of variables, and variables that lack obvious relationships with other variables.
Boston Housing and Iris are two standard data sets from sklearn. The housing data set contains 506
samples of 14 variables, originally collected to build a model for predicting house prices (descriptions
of the variables can be found online6 ). The Iris data set is a data set for predicting flower species
based on petal and sepal measurements with 150 observations in 5 dimensions.
Estimation will take place using the default stacking estimator. Since we are interested in initial
exploration, and finding interesting groups of variables, a large FDR (20%) was chosen. Note that,
unlike for other (mostly score-based) structure learning algorithms, the outcome of one experiment
(the presence of an edge in the model) does not influence other experiments, and hence, false discoveries do not compromise the estimated structure additionally.
The results are shown in Figure 6 (graphs drawn with NetworkX [24]). For the Boston housing data,
seemingly sensible variable groupings occur. In the top left, the variables related to industrialization
of a neighbourhood are shown, while on the right, demographic attributes form a cloud of related
variables. For the Iris data set, while length and width of both sepal and petal are related, as
expected, it seems that petal size has a higher association with class, and, in fact, width and length
are independent given the class.
Both of these analyses require no parameter tuning (the only non-default chosen for this experiment
was the adjusted confidence level) and take very little time (less than 15 seconds). The implementation
of algorithm 3 thus provides a quick, ready to use tool for initial exploratory data analysis.
7.2.2. Key short-term economic indicators (UK)
Setup: One is interested in finding the relationships within a set of variables and making local
conditional independence statements for sets of variables. The focus is on finding associations that
we are confident about.
The economic indicator data set contains monthly data of key economic indicators between February
6 http://scikit-learn.org/stable/modules/classes.html#module-sklearn.datasets

41

(a) Default Approach

(b) Bagged Support Vector Machines

Figure 7: Learned graphs for the Economic Indicator data set for the default approach (left) and a
more powerful approach using bagged SVM’s (right)
1987 and June 2017 for the UK from the OECD statistics database [30]. The economics indicators are
export and import figures, stock price and industrial production levels, overnight and 3 month Interbank interest rates, money supply (M4) and GDP. This rather small data set, around 369 instances
of 9 variables, will outline the possibility of a user to trade off between computational complexity
and power of the learning algorithm. As for most economic data sets, the signal to noise ratio can
be quite low. Figure 7a shows the structure learned by the default approach for confidence (false
discovery-rate) level of 5%, with a run time of 15 seconds. While the ground truth is not known
(economic theory aside), it is apparent that many links are missed by the fast default approach. This
becomes evident when considering a variable like the GDP, which (by definition) includes exports and
imports, and hence some connection between the variables should exist. If there is need for a more
powerful approach, a more powerful estimator can be chosen. Figure 7b show the learned structure for
an estimator that resamples the data 50 times and learns 50 different Support Vector Machines with
automatically tuned hyperparameters. The implementation of this is straightforward, since sklearn
provides hyperparameter and bagging wrappers, so an estimator can be tuned easily and passed to the
MetaEstimator. The graph shows the many edges that were found. While it is impossible to judge
the correctness of the graph, it seems that some reasonable groups of variables are found, such as
industrial production, exports and imports (in economic theory, countries seek to balance their trade
deficit), or the grouping of stock prices and the main drivers of inflation, money supply and interest
rates. Additionally, all the variables are connected in some way, which would be expected from a data
set of economic indicators. The computational complexity of this approach was rather high, with a
run time of about 20 minutes, however this shows how easily the user can trade off between power
and complexity,without having a large effect on the false discovery-rate (as shown in Section 7.1).

42

8. Discussion
This paper has introduced a novel way for multivariate conditional independence testing based on
predictive inference and linked to a supervised learning workflow, which addresses many of the current
issues in independence testing:
• Few subjective choices: By connecting the classical task of independence testing to supervised learning, and its well-understood hyperparameter tuning capabilities, there is no need
for the heuristics and manual tuning choices prevalent in current state-of-the-art conditional
independence tests
• Low computational complexity: By linking independence testing to one of the most-researched
fields in data science, predictive modelling, the high level of power in the state-of-the-artprediction methods can be directly benefited from, and any efficiency gains in the latter directly
benefits the former.
• High-dimensional problems: By connecting the problem to the task of supervised learning,
the method easily deals with the multivariate case, largely removing any issues concerning
dimensionality of the data
It is important to note, that some choices remain necessary, as is the case in statistical methodology in
general. Since it is not possible to test for an infinite amount of loss functions and predictive models,
one has to decide on a subset to conduct this test. The larger the subset, the higher the data requirements to arrive at a respective power level. How the outlined methodology differs from the methods
reviewed in section 2 is by outlining a principled way to choose from a set of subjective choices, and
by using a subset of all-star predictive modelling algorithms to ensemble over as a default, a test that
is expected to be able to deal with most of the usual scenarios, given a reasonable sample size.
To validate these claims, the test was bench-marked against current state-of-the-art conditional independence tests, and showed a similar or better performance in regions where the power of the tests
exceeded 10%.
Subsequently, the PCIT was applied in a method for learning the structure of an undirected graph
best describing the independence structure of an input data set. The combination of the new conditional independence test and the structure learning algorithm address some of the current issues in
graphical model structure learning:
• Low computational complexity: While current exact algorithms, such as the PC-algorithm,
often require a number of tests that is exponential in the number of variables, the proposed
algorithm runs in quadratic time in the size of the graph. Additionally, a straightforward
power-complexity trade off is provided
• Exact algorithm: Unlike many scored-based methods, the algorithm does not make any strong
constraints on the underlying distribution and is not subject to local optima
• False-discovery rate control: FDR control is added to provide a powerful tool to control the
fraction of type 1 errors when the number of variables is increasing
Performance tests showed that the false-discovery rate is as advertised and the power of the test
increases constantly in the number of samples. Additionally, consistency under perturbations in the
data set was demonstrated.
The algorithms have been distributed in the pcit package, providing users with an easy-to-use implementation to test for multivariate and conditional independence, as well as to perform graphical
model structure learning with little need for manual parameter tuning. The implementations are
particularly interesting for users
• assessing the value of collecting additional variables for a prediction task
43

• in need of a conditional independence test for multivariate data
• performing exploratory data analysis
• looking for a visual way to communicate the structure in their data set
There are a few ways in which the current implementation can be generalized to make for a more
powerful test. The power of the conditional independence test is directly linked to the power of the
supervised learning algorithm to find the optimal prediction functional and the correct loss. Currently, only the necessary minimum of 2 loss functions is implemented, one for regression tasks and
one for classification tasks, but this can easily be generalized by including more losses, and checking
if the baseline can be beaten with statistical significance in any of them. This would also strengthen
the argument when reasoning about variables being independent, when the null hypothesis cannot be
rejected. What’s more, the current methodology connected single-output prediction with FDR control
to make statements about the behaviour between the joint distributions. While this procedure results
in a logically sound test, the feasibility of multi-output predictors, predicting several outputs at once,
should be assessed, for appropriate multivariate loss functions.
On the other side, some extensions of the tests need to be conducted before a definitive conclusion
can be made as to its power. In terms of performance evaluation, the power of the proposed routine
was assessed for two specific cases only, for multivariate Gaussian data, and a simple conditional independence test. To get a better idea of the general feasibility of the algorithm, more scenarios need
to be analyzed and performance tests conducted. Additionally, the power of the test in the context
of alternative graphical model structure learning algorithms should be evaluated.

44

A. Best uninformed predictors: classification
This appendix collects proofs of number of elementary computations to obtain some of the explicit
statements about classification losses found in the main corpus.

A.1. Misclassification loss is a probabilistic loss
In this sub-section, we consider the alternative misclassification loss
L : Y0 × Y → [0, 1];

(p, y) 7→ 1 − p(y),

with Y being a discrete set, and Y0 being the probability simplex of probability mass functions on Y,
as considered in Remark 3.2. In said remark, it is claimed that rounding a probabilistic prediction to
a deterministic one never increases the generalization loss. We first prove this for the unconditional
case which is equivalent to constant predictions:
Lemma A.1. Let Y be an Y-valued random variable with probability mass function pY . There is
always a deterministic minimizer of the expected generalization loss, i.e., there exists a pmf p0 : Y →
{0, 1} such that
p0 = argmin E[L(p, Y )].
p∈Y0

Proof. Let p : Y → [0, 1] be any pmf. For its expected generalization loss, one has, by substituting
definitions,
X
X
E[L(p, Y )] =
p(y)(1 − pY (y)) = 1 −
p(y)pY (y).
y∈Y

y∈Y

Let y0 := argmaxpY (y) (if there are multiple maximizers, choose arbitrary). By construction, pY (y0 ) ≥
y∈Y
P
y∈Y p(y)pY (y). Thus, for p0 : y 7→ 1 if y = y0 ; 0 otherwise, one has
E[L(p0 , Y )] = 1 − pY (y0 ) ≤ 1 −

X

p(y)pY (y) = E[L(p, Y )].

y∈Y

Since p was arbitrary, this proves the claim.
Note that Lemma A.1 does not exclude that there are other minimizers of the expected generalization
loss (in general there will be an infinity), it only states that a deterministic minimizer may be found.
Proposition A.2. Let X, Y, Z be a random variables taking values in X, Y, Z. Then one can make
choices which always predict a deterministic class, for the following prediction functionals as considered
in Section 3:
(L)

(i) the best uninformed predictor ωY

: X → Y0

(L)

(ii) the best predictor ωY |X : X → Y0
(L)

(iii) the best conditionally uninformed predictor ωY |Z : X × Z → Y0
(L)

(iv) the best conditional predictor ωY |X,Z : X × Z → Y0
That is, in all four cases, choices may be made where the image is always a deterministic pmf p : Y →
{0, 1}.

45

Proof. (i) follows directly from Lemma A.1 which directly implies that the function y 7→ p0 is the best
constant predictor and hence the best uninformed predictor, where p0 is defined as in the statement
(L)
(or more constructively in the proof) of Lemma A.1, thus ωY : x 7→ p0 is a possible choice.
(ii) follows from noting that the Lemma A.1 and its proof remain still valid when considering the
conditional under X, i.e., defining a conditional
p0 : X → Y0 ; x 7→ argminE[L(p, Y )|X = x]
p∈Y0

(L)

and thus ωY |X : x 7→ p0 (x).
(iii) and (iv) follow analogously by additional conditioning on Z in the same way.

A.2. Logarithmic loss is a proper loss
In this sub-section, we consider the logarithmic loss (or cross-entropy loss)
L : Y0 × Y → R+ ;

(p, y) 7→ − log p(y),

with Y being a discrete set, and Y0 being the probability simplex of probability mass functions on Y,
as considered in Remark 3.2.
Proposition A.3. The expected generalization log-loss is minimized by the true distribution. I.e., let
Y be random variable taking values in Y, with probability mass function pY . Then,
pY = argmin E[L(y, Y )].
y∈Y0

Proof. Let p ∈ Y be arbitrary. Substituting definitions, it holds that
X
E[L(p, Y )] = −
pY (y) log p(y)
y∈Y

≥−

X

p(y) log p(y)

y∈Y

= E[L(pY , Y )],
where the inequality in the middle is Gibbs’ inequality.

A.3. Brier loss is a proper loss
In this sub-section, we consider the Brier loss (or squared classification loss)
X
L : Y0 × Y → R+ ; (p, y) 7→ (1 − p(y))2 +
p(y 0 )2 ,
y 0 6=y

with Y being a discrete set, and Y0 being the probability simplex of probability mass functions on Y,
as considered in Remark 3.2.
Proposition A.4. The expectedBrier loss is minimized by the true distribution. I.e., let Y be random
variable taking values in Y, with probability mass function pY . Then,
pY = argmin E[L(y, Y )].
y∈Y0

46

Proof. Let p ∈ Y be arbitrary. By applying the binomial rule, observe that
X
p(y 0 )2
L(p, y) = 1 − 2p(y) +
y 0 ∈Y

Substituting definitions, it holds that


X
X
pY (y) − 2pY (y)p(y) + pY (y)
p(y 0 )2 
E[L(p, Y )] =
y 0 ∈Y

y∈Y

=1−2

X

pY (y)p(y) +

=1−

2

pY (y) +

y∈Y

=1−

X

p(y)2

y∈Y

y∈Y

X

X

X
y∈Y

2

pY (y) +

y∈Y

X

X

pY (y)2 − 2

pY (y)p(y) +

y∈Y

X

p(y)2

y∈Y

2

(p(y) − pY (y)) .

y∈Y

Note that the first two terms (the 1 and the sum over pY (y)) do not depend on p, while the last term
is non-negative and minimized with a value of zero if and only if p = pY .
Since p was arbitrary, this proves the claim.

B. Elicited statistics for regression losses
This appendix collects explicit proofs of the elicited statistics for squared loss, and Q-loss (and the
absolute loss which is a special case).

B.1. Squared loss elicits the mean
In this sub-section, we consider the univariate regression case Y = R, and the squared loss
L : Y × Y → R; (y, y ∗ ) 7→ (y − y ∗ )2 .
Proposition B.1. The squared loss elicits the mean. I.e., let Y be random variable taking values in
Y. Then,
E[Y ] = argmin E[L(y, Y )].
y∈Y

Proof. Substituting definitions, it holds that
E[L(y, Y )] = E[(y − Y )2 ]
= E[y 2 − 2yY + Y 2 ]
= y 2 − 2yE[Y ] + E[Y 2 ]
= y 2 − 2yE[Y ] + E[Y 2 ] − E[Y 2 ] + E[Y 2 ]
= (E[Y ] − y)2 + Var(Y ),
which is the well-known derivation of the bias-variance decomposition.
The first term, (E[Y ] − y)2 , is minimized whenever E[Y ] = y, and the second term, Var(Y ), does
not depend on y. Thus, the sum of both terms is minimized (in y while fixing Y ) for the choice
y = E[Y ].

47

B.2. Quantile loss elicits the quantile
In this sub-section, we consider the univariate regression case Y = R, and the Q-loss (or quantile loss)
L : Y × Y → R; (y, y∗ ) 7→ α · m(y∗ , y) + (1 − α) · m(y, y∗ ),
with m(x, z) = min(x − z, 0).
Proposition B.2. The Q-loss elicits the α-quantile. I.e., let Y be random variable taking values in
Y with cdf F : R → [0, 1]. Then,
FY−1 (α) = argmin E[L(y, Y )].
y∈Y

Proof. We first assume that Y is absolutely continuous, i.e., Y has a probability density function
p : R → R+ and F is bijective. One then computes
EY [Lα (y ∗ , Y )] =

Z

y∗

(1 − α)(y ∗ − y)p(y)dy −

Z

(α)(y ∗ − y)p(y)dy

y∗

−∞
∗

∞

Z

∗

y∗

= y (F (y ) − α) + αE[Y ] −

yp(y)dy
−∞

∂EY [Lα (y ∗ , Y )]
!
= y ∗ p(y ∗ ) + F (y ∗ ) − α − y ∗ p(y ∗ ) = F (y ∗ ) − α = 0
∗
∂y
=⇒ P (y ∗ ) = α
∂ 2 EY [Lα (y ∗ , Y )]
= p(y ∗ ) ≥ 0
∂(y ∗ )2
Hence, the first order condition is a minimum, minimized by the α-quantile of Y , and thus the quantile
loss elicits the quantile.
For general Y , note that F always exists, and thus when p appears inside integrals, the integrals welldefined. The partial derivatives may not always be defined but is the same as the sign of sufficiently
small finite differences, thus the proof logic follows through for the general case. In case of jump
discontinuities of F , any monotone inverse F −1 may be chosen for the statement.

48

References
[1] C.C. Aggarwal. Data Classification: Algorithms and Applications. Chapman & Hall/CRC, 1st
edition, 2014.
[2] D. Barber. Bayesian reasoning and machine learning. Cambridge University Press, 2012.
[3] L. Baringhaus and C. Franz. On a new multivariate two-sample test. Journal of Multivariate
Analysis, 88:190–206, 2004.
[4] Y. Benjamini and Y. Hochberg. Controlling the false discovery rate: A practical and powerful
approach to multiple testing. Journal of the Royal Statistical Society. Series B (Methodological),
57:289–300, 1995.
[5] Y. Benjamini and D. Yekutieli. The control of the false discovery rate in multiple testing under
dependency. The Annals of Statistics, 29:1165–1188, 2001.
[6] W. Bergsma. Nonparametric testing of conditional independence by means of the partial copula.
ArXiv e-prints: 1101.4607, 2011.
[7] W.P. Bergsma. Testing conditional independence for continuous random variables. Eurandom,
2004.
[8] L. Buitinck, G. Louppe, M. Blondel, F. Pedregosa, A. Mueller, O. Grisel, V. Niculae, P. Prettenhofer, A. Gramfort, J. Grobler, et al. API design for machine learning software: experiences
from the scikit-learn project. arXiv preprint arXiv:1309.0238, 2013.
[9] U. Cherubini, E. Luciano, and W. Vecchiato. Copula methods in finance. John Wiley & Sons,
2004.
[10] K.P. Chwialkowski, A. Ramdas, D. Sejdinovic, and A. Gretton. Fast two-sample testing with
analytic representations of probability measures. In Advances in Neural Information Processing
Systems 28, pages 1981–1989. Nips, 2015.
[11] G.W. Corder and D.I. Foreman. Nonparametric Statistics for Non-Statisticians: A Step-by-Step
Approach. John Wiley & Sons, Inc., 2009.
[12] A. Dionisio, R Menezes, and D.A. Mendes. Entropy-based independence test. Nonlinear Dynamics, 44:351–357, 2006.
[13] V.A. Fernández, M.D.J. Gamero, and J.M. Garcı́a. A test for the two-sample problem based
on empirical characteristic functions. Computational Statistics & Data Analysis, 52:3730–3748,
2008.
[14] T.M. Fragoso and F.L. Neto. Bayesian model averaging: A systematic review and conceptual
classification. arXiv preprint arXiv:1509.08864, 2015.
[15] J. Friedman, T. Hastie, and R. Tibshirani. The elements of statistical learning, volume 1. Springer
series in statistics New York, 2001.
[16] C. Genest and B. Rémillard. Test of independence and randomness based on the empirical copula
process. Test, 13:335–369, 2004.
[17] T. Gneiting and A.E. Raftery. Strictly proper scoring rules, prediction, and estimation. Journal
of the American Statistical Association, 102:359–378, 2007.
[18] A. Gretton and L. GyĂśrfi. Consistent nonparametric tests of independence. Journal of Machine
Learning Research, 11:1391–1423, 2010.
49

[19] A. Gretton, O. Bousquet, A. Smola, and B. Scholkopf. Measuring statistical dependence with
hilbert-schmidt norms. In ALT, volume 16, pages 63–78. Springer, 2005.
[20] A. Gretton, K. Fukumizu, C.H. Teo, L. Song, B. Schölkopf, and A. Smola. A kernel statistical
test of independence. In Advances in neural information processing systems, pages 585–592, 2008.
[21] A. Gretton, K. Fukumizu, Z. Harchaoui, and B.K. Sriperumbudur. A fast, consistent kernel
two-sample test. Advances in Neural Information Processing Systems 22, pages 673–681, 2009.
[22] A Gretton, K.M. Borgwardt, M.J. Rasch, B. Schölkopf, and A. Smola. A kernel two-sample test.
Journal of Machine Learning Research, 13:723–773, 2012.
[23] A. Gretton, D. Sejdinovic, H. Strathmann, S. Balakrishnan, M. Pontil, K. Fukumizu, and B.K.
Sriperumbudur. Optimal kernel choice for large-scale two-sample tests. In Advances in Neural
Information Processing Systems 25, pages 1205–1213. Nips, 2012.
[24] A.A. Hagberg, D.A. Schult, and P.J. Swart. Exploring network structure, dynamics, and function
using NetworkX. In Proceedings of the 7th Python in Science Conference (SciPy2008), pages 11–
15, 2008.
[25] E. Jones, T. Oliphant, and P. Peterson. SciPy: Open source scientific tools for Python, 2001.
URL http://www.scipy.org/.
[26] D. Koller and N. Friedman. Probabilistic graphical models: principles and techniques. MIT press,
2009.
[27] M. Lichman. UCI machine learning repository, 2013. URL http://archive.ics.uci.edu/ml.
[28] D. Lopez-Paz and M. Oquab.
arXiv:1610.06545, 2016.

Revisiting classifier two-sample tests.

arXiv preprint

[29] Y. Nadeau, C.and Bengio. Inference for the generalization error. Machine Learning, 52:239–281,
2003.
[30] OECD. OECD statistics: Key short-term economic indicators. http://stats.oecd.org/, 2017.
Accessed: 2017-08-06.
[31] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning
Research, 12:2825–2830, 2011.
[32] S. Raschka. Mlxtend, 2016. URL http://dx.doi.org/10.5281/zenodo.594432.
[33] B. Rémillard and O. Scaillet. Testing for equality between two copulas. Journal of Multivariate
Analysis, 100:377–386, 2009.
[34] B. Schweizer and E. F. Wolff. On nonparametric measures of dependence for random variables.
The Annals of Statistics, 9:879–885, 1981.
[35] Bharath K Sriperumbudur, Kenji Fukumizu, Arthur Gretton, Gert RG Lanckriet, and Bernhard
Schölkopf. Kernel choice and classifiability for RKHS embeddings of probability distributions. In
NIPS, pages 1750–1758, 2009.
[36] W. Zaremba, A. Gretton, and M. Blaschko. B-test: A non-parametric, low variance kernel twosample test. Advances in Neural Information Processing Systems 26, pages 755–763, 2013.
[37] K. Zhang, J. Peters, D. Janzing, and B. Schölkopf. Kernel-based conditional independence test
and application in causal discovery. arXiv preprint arXiv:1202.3775, 2012.
50

