Faculty of Science, Technology and Communication

The Effect of Noise Level on Causal Identification with Additive Noise Models

Thesis Submitted in Partial Fulfillment of the Requirements for the Degree of
Master in Information and Computer Sciences

_Author:_

Benjamin Kap

_Supervisor:_
Prof. Thomas Engel
_Reviewer:_
Prof. Nicolas Navet
_Advisor:_
Dr. Marharyta Aleksandrova

May 2021

## Abstract

In recent years a lot of research has been conducted within the area of causal
inference and causal learning. Many methods have been developed to identify
the cause-effect pairs in models and have been successfully applied to
observational real-world data in order to determine the direction of causal
relationships. Many of these methods require simplifying assumptions, such as
absence of confounding, cycles, and selection bias. Yet in bivariate
situations causal discovery problems remain challenging. One class of such
methods, that also allows tackling the bivariate case, is based on Additive
Noise Models (ANMs). Unfortunately, one aspect of these methods has not
received much attention until now: what is the impact of different noise
levels on the ability of these methods to identify the direction of the causal
relationship. This work aims to bridge this gap with the help of an empirical
study. For this work, we considered bivariate cases, which is the most
elementary form of a causal discovery problem where one needs to decide
whether $X$ causes $Y$ or $Y$ causes $X$, given joint distributions of two
variables $X$, $Y$. Furthermore, two specific methods have been selected,
Regression with Subsequent Independence Test and Identification using
Conditional Variances, which have been tested with an exhaustive range of ANMs
where the additive noises’ levels gradually change from 1% to 10000% of the
causes’ noise level (the latter remains fixed). Additionally, the experiments
in this work consider several different types of distributions as well as
linear and non-linear ANMs. The results of the experiments show that these
methods can fail to capture the true causal direction for some levels of
noise.

###### Contents

  1. Abstract
  2. 1 Introduction
    1. 1.1 State of the Art
  3. 2 Additive Noise Models and Notation
    1. 2.1 Experiment Layout
  4. 3 Experiments on Regression with Subsequent Independence Test
    1. 3.1 RESIT with Different Noise Levels
      1. 3.1.1 Setup
      2. 3.1.2 Execution
      3. 3.1.3 Experimental Results
    2. 3.2 RESIT with Different Noise Levels and coupled estimation
      1. 3.2.1 Results
    3. 3.3 RESIT with Different Noise Levels without prior assumption
      1. 3.3.1 Setup
      2. 3.3.2 Execution
      3. 3.3.3 Results
    4. 3.4 RESIT with different means
      1. 3.4.1 Setup
      2. 3.4.2 Execution
      3. 3.4.3 Results
  5. 4 Experiments on Identification using Conditional Variances
    1. 4.1 Introduction
      1. 4.1.1 Algorithm
    2. 4.2 Setup
      1. 4.2.1 Execution
      2. 4.2.2 Results
  6. 5 Conclusion and Future Work

###### List of Figures

  1. 2.1 A simple directed acyclic graph showing causal relationships between variables.
  2. 2.2 Another directed acyclic graph showing causal relationships between variables.
  3. 3.1 Data is generated following $Y=X+N_{y}$ with uniform distributions only.
  4. (a) Regress $Y$ on $X$.
  5. (b) Regress $X$ on $Y$.
  6. 3.2 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{N}+\mathcal{N}$
  7. 3.4 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{N}+\mathcal{L}$
  8. 3.6 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{N}^{3}+\mathcal{U}$
  9. 3.8 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{U}+\mathcal{U}$
  10. 3.10 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{U}+\mathcal{L}$
  11. 3.12 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{U}^{3}+\mathcal{N}$
  12. 3.14 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{L}+\mathcal{L}$
  13. 3.16 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{L}+\mathcal{U}$
  14. 3.18 RESIT & different noise levels & decoupled estimation & $Y=\mathcal{L}^{3}+\mathcal{N}$
  15. 3.20 RESIT & different noise levels & coupled estimation & $Y=\mathcal{N}+\mathcal{N}$
  16. 3.22 RESIT & different noise levels & coupled estimation & $Y=\mathcal{N}+\mathcal{L}$
  17. 3.24 RESIT & different noise levels & coupled estimation & $Y=\mathcal{N}^{3}+\mathcal{U}$
  18. 3.26 RESIT & different noise levels & coupled estimation & $Y=\mathcal{U}+\mathcal{U}$
  19. 3.28 RESIT & different noise levels & coupled estimation & $Y=\mathcal{U}+\mathcal{L}$
  20. 3.30 RESIT & different noise levels & coupled estimation & $Y=\mathcal{U}^{3}+\mathcal{N}$
  21. 3.32 RESIT & different noise levels & coupled estimation & $Y=\mathcal{L}+\mathcal{L}$
  22. 3.34 RESIT & different noise levels & coupled estimation & $Y=\mathcal{L}+\mathcal{U}$
  23. 3.36 RESIT & different noise levels & coupled estimation & $Y=\mathcal{L}^{3}+\mathcal{N}$
  24. 3.38 $Y=X+N_{y}$. Contains all cases where $X$ and $N_{y}$ are drawn from the same type of distribution. Dashed lines are non-linear cases. Decoupled estimation (80% split).
  25. 3.40 $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{U}$ and $N_{y}\cancel{\sim}\mathcal{U}$. Dashed lines are non-linear cases. Decoupled estimation (80% split).
  26. 3.42 $Y=X+N_{y}$. Contains all cases where $X$ and $N_{y}$ are drawn from the same type of distribution. Dashed lines are non-linear cases. Coupled estimation.
  27. 3.44 $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{U}$ and $N_{y}\cancel{\sim}\mathcal{U}$. Dashed lines are non-linear cases. Coupled estimation.
  28. 3.46 Left: $X\sim\mathcal{N},N_{y}\sim\mathcal{U}$. Right: $X\sim\mathcal{U},N_{y}\sim\mathcal{L}$.
  29. 4.1 $Y=X+N_{y}$. Contains all linear cases.

##  Chapter 1 Introduction

Due to the technological and computational advances during the last decades,
scientist were able to tackle non-trivial problems from different research
areas successfully. One of these research areas is causality. The fundamentals
of causality is to determine causal relationship between two or more variables
in a system. For example, given altitude and temperature we want to answer the
question if temperature has an effect on altitude, or if altitude has an
effect on temperature. This is of particular interest since if such a causal
relationship is known then one can predict the effects on a system in case of
intervention or perturbation. One method to determine causal relationships in
a system is controlled experimentation (A/B tests) in which there are two
identical groups with only one variation. The only variable that is varied
(intervened on) is the potential cause. This procedure allows estimating
causal effect of this variable in the given system. For example, testing the
efficacy of medications is done within A/B tests. The control group receives
no medication or a placebo, and the intervention group receives the real
medication. The results often show the true effect (if any) of the medication
on the human health. However, such tests are often too expensive, unethical or
even technically impossible to execute. Therefore, it is of great interest to
determine causal relationships from observational data only (e.g., structure
learning).

###  1.1 State of the Art

Structure learning is the procedure to determine causal relationship
directions from observational data only and representing these as a (causal)
graph. The basic idea emerged from [1] as path analysis. This is used to
describe directed dependencies among a set of variables and includes various
models such as ANOVA, ANCOVA, etc. In his work, Wright made a distinction
between three possible types of causal substructures which were allowed in a
directed acyclic (no cycles) graph: 1) $X\to Y\to Z$, 2) $X\leftarrow Y\to Z$,
and 3) $X\to Y\leftarrow Z$. [2] developed an algorithm to recover directed
acyclic graphs from statistical data, which relied on this distinction of the
previously 3 mentioned substructures. In general, one can easily identify the
skeleton of a graph (that is the graph without arrows on the edges) and then
partially identify the arrows (partially, because the three substructures’
skeletons are identical but only 3) is distinguishable form others). [3, 4]
used Bayes networks to axiomatize the connection between causal structure and
probabilistic independence and formalized under what assumptions one could
draw causal knowledge from observational data only. Furthermore, they also
formalized how incomplete causal knowledge could be used for causal
intervention. Judea Pearl presented in his work [5] a comprehensive theory of
causality and unified the probabilistic, manipulative, counterfactual, and
structural approaches to causation. Judea also introduced precise mathematical
definitions of causal analysis for the standard curricula of statistics. From
the work [5] we have the following key point: if there is a statistical
association, e.g. two variables $X,Y$ are dependent, then one of the following
is true: 1) there is a causal relationship, either $X$ has an effect on $Y$ or
$Y$ has an effect on $X$; 2) there is a common cause (confounder) that has
effect on both $X$ and $Y$; 3) there is a possibly unobserved common effect of
$X$ and $Y$ that is conditioned upon data acquisition (selection bias); or 4)
there can be a combination of these. From there on a lot of research has been
conducted to develop theoretical approaches and methods for identifying causal
relationships from observational data only. Before we present some major works
from the last two decades, we introduce the common concept behind all these
approaches in a short formal manner.

In general, all these methods exploit the complexity of the marginal and
conditional probability distributions in some way (e.g., [6, 7]) and under
certain assumptions these methods are then able to solve the task of causal
discovery. Let $C$ denote the cause and $E$ the effect. In a system with two
or more variables we might have cause-effect pairs and then their joint
density can be expressed with $p_{C,E}(c,e)$. This joint density can be
factorized into either 1) $p_{C}(c)\cdot P_{E|C}(e|c)$ or 2) $p_{E}(e)\cdot
P_{C|E}(c|e)$. The idea is then that 1 gives models of lower total complexity
than 2 and this allows us to draw conclusions about the causal relationship
direction. Intuitively this makes sense, because the effect contains
information from the cause but not vice-versa (of course under the assumption
that there are no cycles aka feedback loops). Therefore, 2 has at least as
much complexity as 1. However, the definition of complexity is ambiguous. For
example, one can say that "$p_{C}$ contains no information about
$P_{E|C}(e|c)$" and then is able to draw partial conclusions about the causal
direction of the given system. This complexity question is often colloquially
referred to as "breaking the symmetry" (that is $p_{C}(c)\cdot
P_{E|C}(e|c)\neq p_{E}(e)\cdot P_{C|E}(c|e)$).

[8] addressed the problem of learning the structure of a Bayesian network in
domains which contain continuous variables. In their work they showed that in
probabilistic networks with continuous variables one can use Gaussian Process
priors to compute marginal likelihoods for structure learning.
[9] developed a model for causal inference using non-normality of observed
data and improved path analysis [1] using non-normal data. [10] proposed a
method on how to determine the complete causal graph of continuous data under
three assumptions: the data generating process is linear, no unobserved
confounders, and noise variables have non-Gaussian distributions of non-zero
variances. This method was not scale-invariant, but in a later work [11] this
problem had been addressed and a new method was proposed which was guaranteed
to converge to the right solution within a small fixed number of steps if the
data strictly followed the model.
[12] introduced a method based on comparing the conditional distributions of
variables given their direct causes for all hypothetical causal directions and
choosing the most plausible one (Markov kernels). Those Markov kernels which
maximize the conditional entropies constrained by their observed expectation,
variance and covariance with its direct causes based on their given domain are
considered as plausible kernels. [13] continued the work on kernels by using
the concept of reproducing kernel Hilbert spaces.
[14] generalized the linear framework of additive noise models to nonlinear
models. Additive noise models are models in which the effect is a function of
the cause and some random and non-observed additive noise term. For continuous
variables, earlier methods often assumed linear models for the independence
tests. However, if data contained non-Gaussian variables, then this can help
in distinguishing the causal directions and identify the causal graph.
[15] proposed a method for inferring linear causal relationships among multi-
dimensional variables by factorizing the joint distribution into products with
marginal and conditional distributions (as seen above with 1 and 2). Then, in
one of these products the factors (e.g., $P(E)$ and $P(C|E)$) satisfy non-
generic relations indicating that $E\to C$ is wrong.
[16] introduced a method which minimizes the statistical dependence between
the regressors and residuals. If residuals (the difference between actual
output and predicted output) are no longer dependent on the input, then
regression can successfully model the dependence of the output on the input.
This method does not need to assume a particular distribution of the noise
because any form of regression can be used (e.g., Linear Regression) and is
well suited for the task of causal inference in additive noise models.
[17] created a method to model observed data by using probabilistic latent
(hidden) variable models, which incorporate the effects of unobserved noise.
To analyze the joint density of cause and effect, the effect is modeled as a
function of the cause and some independent noise (not necessarily additive).
With general non-parametric priors on this function and on the distribution of
the cause the causal direction is then determined by using standard Bayesian
model selection.
[18] introduced a method to determine causal relationship in cyclic additive
noise models and state that such models are generally identifiable in the
bivariate, Gaussian-noise case. Their method works for continuous data and can
be seen as a special case of nonlinear independent component analysis.
[19] tested the stability of the identifiability of post-nonlinear models with
two variables and listed all cases in which this model is not identifiable
anymore. Furthermore, they showed how to approach multivariate cases with
post-nonlinear models.
[20] showed that even in deterministic cases (noise-free cases) there are
asymmetries that can be exploited for causal inference. Their method is based
on the idea that if $X\to Y$ then the distribution of $X$ and the function
mapping $X$ to $Y$ must be independent since they correspond to independent
mechanisms of nature.
[21] proposed a method which is based on the likelihood ratio under the linear
non-Gaussian acyclic model (LiNGAM, [22]) and therefore not resorting to
independent component analysis algorithm as previous methods did.
[23] proved full identifiability of linear Gaussian structural equation models
if all the noise variables have the same variance (full identifiability means
that not only the skeleton of the causal graph is recoverable but also the
arrows).
[24] proposed a method that can identify the directed acyclic graph from the
distribution under mild conditions. In contrast, other methods assumed
faithfulness and could only identify the Markov equivalence class of the graph
(Markov equivalence class refers to the class of graphs where all graphs have
the same skeletons).
[25] proposed to use penalized likelihood scores instead of independence
scores to determine the true causal graph from the Markov equivalence class.
[26] and [27] proved that linear Gaussian models where noise variables have
different error variance can be identifiable by ordering variables according
to the law of total variances and then performing independence tests between
variables.

Despite all the research in the past years one small but nonetheless important
aspect of causal discovery methods has not received much attention: can
different noise levels have an impact on the correctness of these methods. In
real world, observational data often differs in terms of noise level. Usually,
these levels do not differ significantly but it can occur that noise levels
change drastically from cause to effect. For example, if the data collection
process has a lot of interference (e.g., in outer space) then such noise
levels can differ a lot.

In this work we will focus on the Additive Noise Models (ANMs) as they are
well established and yielded many good results ([28]). ANMs are heavily based
on the presence of noise and thus our research question is then formulated as
follows: how do different noise levels of the additive noise impact the
correctness of ANM methods? For this work, two specific methods have been
selected to test this question: Regression with Subsequent Independence Test
(RESIT) [24] and Identification using Conditional Variances (Uncertainty
Scoring) [29]. Both methods yielded good results and can be used even when
variables have different distribution types (e.g., Laplace). The next chapter
introduces Additive Noise Models in a formal manner, followed by standard
notation definition. Chapters 3 and 4 contain the theory and experiments of
the first and the second methods respectively. The last chapter draws
conclusions and discusses possible future work.

##  Chapter 2 Additive Noise Models and Notation

The two selected methods for this work can be used for additive noise models,
although the second method is not limited to it ([26]). In general, such
causal discovery methods exploit the additivity of the noise in order to
determine the causal relationship between two or more variables. We will give
now a short definition of additive noise models for the bivariate case only.
For more details and for multivariate cases please refer to [14] and [24].

Let $X,Y\in\mathbb{R}$ be the cause and effect, respectively. Let there be
also $m$ latent (hidden) causes $U=(U_{1},\dots,U_{m})\in\mathbb{R}^{m}$. Then
the causal relationship can be modeled as

$\begin{cases}Y=f(X,U_{1},\cdots,U_{m})\\\ X\perp\\!\\!\\!\perp U,X\sim
p_{X}(x),U\sim p_{U}(u_{1},\cdots,u_{m})\end{cases}$

where $f:\mathbb{R}\times\mathbb{R}^{m}\to\mathbb{R}$ is a linear or nonlinear
function, and $p_{X}(x)$ and $p_{U}(u_{1},\cdots,u_{m})$ are the joint
densities of the observed cause $X$ and the latent causes $U$. We are assuming
that there is no confounding, no selection bias, and no feedback loop between
$X$ and $Y$ and therefore $X$ and $U$ are independent, which is denoted as
$X\perp\\!\\!\\!\perp U$. Since the latent causes $U$ are unobserved, their
influence can be summarized by a single noise variable $N_{y}\in\mathbb{R}$,
and the model can be defined as follows:

$\begin{cases}Y=f(X,N_{y})\\\ X\perp\\!\\!\\!\perp N_{y},X\sim
p_{X}(x),N_{y}\sim p_{N_{y}}(n_{y}).\end{cases}$

##### Notations

In the experiments, we are considering both linear and nonlinear additive
noise models:

$Y=\beta X+N_{y}\text{ with }\beta\in\mathbb{R}\text{, for the linear case}$

and

$Y=\beta X^{\alpha}+N_{y}\text{ with }\beta,\alpha\in\mathbb{R}\text{, for the
nonlinear case}.$

The arrow in "$X\to Y$" signifies $X$ has an effect on $Y$, or in other words,
$X$ is a cause of $Y$. $X$ and $N_{y}$ can be drawn from one of the following
distributions: the normal distribution denoted by the calligraphic letter
$\mathcal{N}$, the uniform distribution denoted by the calligraphic letter
$\mathcal{U}$, or the laplace distribution denoted by the calligraphic letter
$\mathcal{L}$. For example, throughout this work "X is drawn from a normal
distribution" is denoted as $X\sim\mathcal{N}$ or
$X\sim\mathcal{N}(\mu_{x},\sigma^{2}_{x})$ with $\mu_{x}$ standing for the
mean and $\sigma^{2}_{x}$ for the variance. The Structural causal model of a
directed acyclic graph (DAG) refers to the equations of the cause and effects.
For example, the equation $Y=f(X,N_{y})$ refers to the DAG in Fig. 2.1 and
also composes the entire structural equation model of that DAG (since there
are no other variables in that DAG). In Fig. 2.2 the structural equation model
looks as follows:

$\begin{cases}A=f(Y,N_{a})\\\ B=f(A,N_{b})\\\ X=f(B,N_{x})\end{cases}$

Figure 2.1: A simple directed acyclic graph showing causal relationships
between variables. Figure 2.2: Another directed acyclic graph showing causal
relationships between variables.

As we are only considering bivariate cases, we only have one equation and
therefore structural equation model, structural equation and equation are used
interchangeably. Lastly, we introduce the four following terms (with
abbreviations): parent (Pa), antecedent (An), descendent (De) and non-
descendent (Nd). These notations are used to describe the relationship between
two variables in a causal graph. Parent is the direct antecedent of a node
(e.g., $Pa(A)=Y$ in Fig. 2.2). A descendent of a (starting) variable is any
variable which can be reached from the starting variable following a directed
path (e.g., $De(Y)=\\{A,B,X\\}$ in Fig. 2.2). Non-descendent is the opposite,
i.e. cannot be reached following a directed path (e.g.,
$Nd(Y)=\\{N_{a},N_{b},N_{x}\\}$). Finally, the antecedent is any starting
variable for a descendent (e.g., $An(X)=\\{Y,A,B,N_{x},N_{a},N_{b}\\}$).

###  2.1 Experiment Layout

One additional part of the experiments is to see whether splitting data into
training and test data is beneficial or not. However, as some experiments did
not allow us to split data and/or the results were too small to be contained
in an own subsection it is not quite obvious from the contents table where
data splitting occurred or not. The following Table 2.1 thus shall provide an
overview where we applied splitting data and where not. Furthermore,
throughout this work splitting data is referred to as decoupled estimation and
not splitting data is referred to as coupled estimation, both terms originate
from [28].

Section | Type | Decoupled | Coupled
---|---|---|---
Section 3.1 | RESIT & different noise levels | ✓ |
Section 3.2 | RESIT & different noise levels |  | ✓
Section 3.3 |  RESIT & different noise levels & no prior assumption | ✓ | ✓
Section 3.4 | RESIT & different means | ✓ |
Chapter 4 | Uncertainty Scoring |  | ✓
Table 2.1: Experiments Layout for Decoupled and Coupled estimation.

We have uploaded all source codes and results on a Gitlab repository:
https://gitlab.com/Shinkaiika/noise-level-causal-identification-additive-
noise-models/

##  Chapter 3 Experiments on Regression with Subsequent Independence Test

This chapter is an empirical study which involves testing the performance of
Regression with Subsequent Independence Test (RESIT) methods [24] for additive
noise models (ANM) [14, 28, 30] given i.i.d data from a joint distribution to
try to estimate the corresponding directed acyclic graph (DAG). This chapter
is divided into four sections. The first three sections (Section 3.1 to
Section 3.3) contain the main work of this chapter. In Section 3.1 we test
RESIT with different noise levels and decoupled estimation. In Section 3.2 we
test RESIT with different noise levels and coupled estimation. In Section 3.3
we test RESIT with different noise levels and both decoupled and coupled
estimation, but without the prior assumption (which says that exactly one
causal direction must be present in the bi-variate case). In Section 3.4 we
perform tests on RESIT with different means (for both the cause variable and
the noise term). For all experiments we generate artificial data using linear
and non-linear functions. While both linear and non-linear data can be
identifiable in causal models, non-linearity helps in identifying the causal
direction as was shown by [14]. Therefore, we consider both cases in our
experiments in order to investigate how different noise levels affect causal
discovery for linear and on non-linear data. In all experiments we use the
equation $Y=X+N_{Y}$ for the linear cases and $Y=X^{3}+N_{Y}$ for the non-
linear cases. These two structural causal models have been selected
arbitrarily for simplicity. For the consistency of the identifiability of
linear and non-linear data in additive noise models, the reader is referred to
[28, 10, 14, 19]. Finally, in both sections we will also measure performance
difference when deploying decoupled estimation (splitting data into training
and test data) and coupled estimation (no splitting), see [28], [30].

###  3.1 RESIT with Different Noise Levels

The first section of experiments involves testing different noise levels for
the additive noise term and to see whether it has an impact on the accuracy of
RESIT methods or not. The RESIT method is based on the fact that for each node
$X_{i}$ the corresponding noise variable $N_{i}$ is independent of all non-
descendants of $X_{i}$. For example, in a DAG, if we have $Y=X_{1}+N_{1}$ then
$X_{1}\perp\\!\\!\\!\perp N_{1}$. Following this idea we visit every single
node in a specific order. To determine this order we do an iterative process
where in each step we determine the next node in the order. To determine the
next node we look for a sink node. This is done by regressing each single
variable on all other variables and then measuring the independence between
the residuals and those other variables, and finally selecting the single
variable which led to the least dependent residuals. After this iterative
process, the order is established and every node is visited again. We then
eliminate incoming edges until the residuals are not independent anymore. We
will restrict our experiments to bivariate cases only to reduce runtimes. The
authors in [24] generalized RESIT to multivariate cases and proved
identifiability. In our experiments we then have two variables, $X$ and $Y$,
and the task is to determine whether $X$ causes $Y$ ($X\to Y$) or $Y$ causes
$X$ ($Y\to X$).

####  3.1.1 Setup

For all empirical tests we assume $X$ to be a cause of $Y$, that is $X\to Y$.
In the sense of additive noise models, the equation is then:

$Y=\beta X+N_{y}$ $(Y=\beta X^{3}+N_{y}\text{ for the non-linear case})$

where

$\beta=1,$

and

$X\sim\begin{cases}\mathcal{N}(0,1)&\text{or}\\\
\mathcal{U}(-1,1)&\text{or}\\\ \mathcal{L}(0,1)\end{cases}$

and

$N_{y}\sim\begin{cases}\mathcal{N}(0,1\cdot i)&\text{or}\\\
\mathcal{U}(-1\cdot i,1\cdot i)&\text{or}\\\ \mathcal{L}(0,1\cdot
i)\end{cases}$

with $i$ being a scaling factor for the noise level in $N_{y}$. The goal is to
analyze how different standard deviations (boundaries for the uniform case) in
the noise term $N_{y}$ relative to the standard deviation (or boundaries for
the uniform case) in the $X$ term impact the RESIT method.

We apply the same algorithm as Algorithm 1 in [30] which requires inputs $X$
and $Y$, a regression method and a score estimator
$\hat{C}:\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}$ and outputs $dir$
(casual relationship direction). First, the data is split into training data
(80%) and test data (20%). [28] refers to this as "decoupled estimation". The
training data is used to fit the regression model and the test data is used
for the estimator. The idea is to regress $Y$ on $X$ with the training data,
predict $\hat{Y}$ with the test data and then calculate residuals
$Y_{res}=\hat{Y}-Y_{Test}$. $Y_{res}$ and $X_{Test}$ are then used in the
criterion step, where we use several estimators independently and receive
$\hat{C}_{X\to Y}$, a score for the assumed case $X\to Y$. Similarly, to test
the other case ($Y\to X$), we regress $X$ on $Y$, calculate residuals
$X_{res}=\hat{X}-X_{Test}$ and estimate $\hat{C}_{Y\to X}$. In our test
scenario our generated data always follows $X\to Y$. This verifies the
assumption that only one direction in our data is correct (and not both) and
therefore we can compare both scores directly in order to make a decision on
the cause-effect direction. Thus, for independence tests, we can compare
estimates directly and we do not need to determine the value of $\alpha$ for
the independence tests. Algorithm 1 shows pseudo code of the procedure
explained above.

Algorithm 1 General procedure to decide whether $p(x,y)$ satisfies Additive
Noise Model $X\to Y$ or $Y\to X$ with decoupled estimation.

1:Input:

2: 1) I.i.d. sample data $X$ and $Y$

3: 2) Regression method

4: 3) Score estimator
$\hat{C}:\mathbb{R}^{N}\times\mathbb{R}^{N}\to\mathbb{R}$

5:Output:

6: dir

7:

8:Procedure

9:1) Split data into training and test data:

10: $X_{Train},X_{Test}\leftarrow X$

11: $Y_{Train},Y_{Test}\leftarrow Y$

12:

13:2) Train regression models

14: $reg_{1}\leftarrow$ Regress $Y_{Train}\text{ on }X_{Train}$

15: $reg_{2}\leftarrow$ Regress $X_{Train}\text{ on }Y_{Train}$

16:

17:3) Calculate Residuals:

18: $Y_{res}=reg_{1}.predict(X_{Test})-Y_{Test}$

19: $X_{res}=reg_{2}.predict(Y_{Test})-X_{Test}$

20:

21:4) Calculate Scores:

22: $\hat{C}_{X\to Y}=\hat{C}(X_{Test},Y_{res})$

23: $\hat{C}_{Y\to X}=\hat{C}(Y_{Test},X_{res})$

24:

25:5) Finally, output direction dir:  $\displaystyle dir=\begin{cases}X\to
Y&\text{if }\hat{C}_{X\to Y}<\hat{C}_{Y\to X},\\\ Y\to X&\text{if
}\hat{C}_{X\to Y}>\hat{C}_{Y\to X},\\\ ?&\text{if }\hat{C}_{X\to
Y}=\hat{C}_{Y\to X}.\end{cases}$ (3.1)

(a) Regress $Y$ on $X$.

(b) Regress $X$ on $Y$.

Figure 3.1: Data is generated following $Y=X+N_{y}$ with uniform distributions
only.

Fig. 3.1 explains the concept graphically. The left-hand side of the figure
shows a scatter plot of $X_{Test}\text{ and }Y_{res}$. As one can see, given
any input for $X_{Test}$ the range does not change, and therefore $X_{Test}$
and $Y_{res}$ are independent. The right-hand side of the figure shows the
opposite case. For any given input, the range does change. This means that
$Y_{Test}\text{ and }X_{res}$ are dependent and $Y_{Test}$ contains some
information of $X_{res}$ (which contradicts our assumption of independent
noise!). This allows us to draw conclusions about the causal relationship
direction and therefore we know that $X$ is the cause and not vice-versa. This
illustrates the estimation step Eq. 3.1 in Algorithm 1.

For the regression, Linear Regression is used. Linear regression can also be
used in the non-linear case when an appropriate coordinates transformation is
applied. For the estimator several different scores can be used. In our
experiments we used six different independence tests and six different entropy
measures for the estimation criterion. In general, for the independence tests
we have:

$\hat{C}(X_{Test},Y_{res})=I(X_{Test},Y_{res})$

with $I(\cdot,\cdot)$ being any independence test.

In the case of entropy estimators we have:

$\hat{C}(X_{Test},Y_{res})=H(X_{Test})+H(Y_{res}),$

with $H(\cdot)$ being any entropy measure. The estimator score for entropy is
derived from Lemma 1 in [28]:

Lemma 1 Consider a joint distribution of $X,Y$ with density p(x,y). For
arbitrary functions f,g: $\mathbb{R}\to\mathbb{R}$ we have:

$\displaystyle H(X)+H(Y-f(X))=H(Y)+H(X-g(Y))-I(X-g(Y),Y)+I(Y-f(X),X),$ (3.2)

where $H(\cdot)$ denotes differential Shannon entropy, and $I(\cdot,\cdot)$
denotes differential mutual information ([31]).

The above Lemma can be proven with the chain rule of differential entropy (see
[28]). If the density $p(x,y)$ satisfies an identifiable ANM $X\to Y$, then
there exists a function $f$ with $I(Y-f(X),X)=0$ (e.g., the regression
function $x\mapsto\mathbb{E}(Y|X=x)$), but $I(X-g(Y),Y)>0$ for any function
$g$. Therefore, in the case of $X\to Y$, in 3.2 we have:

$H(X)+H(Y-f(X))<H(Y)+H(X-g(Y))$

which is equivalent to the first row in 3.1:

$\hat{C}_{X\to Y}<\hat{C}_{Y\to X}.$

The following estimators were used in this work. The implementation of
estimators with numbers 2 - 12 was taken from the information theoretical
estimators toolbox [32]:111More details on the estimators is given in the
documentation of the toolbox.

  1. 1.

HSIC: Hilbert-Schmidt Independence Criterion with RBF Kernel 222Source:
https://github.com/amber0309/HSIC

$I_{HSIC}(x,y):=||C_{xy}||^{2}_{HS}$

where $C_{xy}$ is the cross-covariance operator and $HS$ the squared Hilbert-
Schmidt norm.

  2. 2.

HSIC_IC: Hilbert-Schmidt Independence Criterion using incomplete Cholesky
decomposition (low rank decomposition of the Gram matrices, which permits an
accurate approximation to HSIC as long as the kernel has a fast decaying
spectrum) which has $\eta=1*10^{-6}$ precision in the incomplete cholesky
decomposition.

  3. 3.

HSIC_IC2: Same as HSIC_IC but with $\eta=1*10^{-2}$.

  4. 4.

DISTCOV: Distance covariance estimator using pairwise distances. This is
simply the $L^{2}_{w}$ norm of the characteristic functions $\varphi_{12}$ and
$\varphi_{1}\varphi_{2}$ of input $x,y$:

$\varphi_{12}(\boldsymbol{u}^{1},\boldsymbol{u}^{2})=\mathbb{E}[e^{i\langle\boldsymbol{u}^{1},\boldsymbol{x}\rangle+i\langle\boldsymbol{u}^{2},\boldsymbol{y}\rangle}],$
$\varphi_{1}(\boldsymbol{u}^{1})=\mathbb{E}[e^{i\langle\boldsymbol{u}^{1},\boldsymbol{x}\rangle}],$
$\varphi_{2}(\boldsymbol{u}^{2})=\mathbb{E}[e^{i\langle\boldsymbol{u}^{2},\boldsymbol{y}\rangle}].$

With $i=\sqrt{-1}$, $\langle\cdot,\cdot\rangle$ the standard Euclidean inner
product and $\mathbb{E}$ the expectation. Finally, we have:

$I_{dCov}(x,y)=||\varphi_{12}-\varphi_{1}\varphi_{2}||_{L^{2}_{w}}$

  5. 5.

DISTCORR: Distance correlation estimator using pairwise distances. It is
simply the standardized version of the distance covariance:

$I_{dCor}(x,y)=\begin{cases}\frac{I_{dCov}(x,y)}{\sqrt{I_{dVar}(x,x)I_{dVar}(y,y)}},&\text{if
}I_{dVar}(x,x)I_{dVar}(y,y)>0\\\ 0,&\text{otherwise,}\end{cases}$

with

$I_{dVar}(x,x)=||\varphi_{11}-\varphi_{1}\varphi_{1}||_{L^{2}_{w}},\>I_{dVar}(y,y)=||\varphi_{22}-\varphi_{2}\varphi_{2}||_{L^{2}_{w}}$

(see characteristic functions under 3. DISTCOV)

  6. 6.

HOEFFDING: Hoeffding’s Phi

$I_{\Phi}(x,y)=I_{\Phi}(C)=\left(h_{2}(d)\int_{[0,1]^{d}}[C(\boldsymbol{u})-\Pi(\boldsymbol{u})]^{2}d\boldsymbol{u}\right)^{\frac{1}{2}}$

with $C$ standing for the copula of the input and $\Pi$ standing for the
product copula.

  7. 7.

SH_KNN: Shannon differential entropy estimator using kNNs (k-nearest
neighbors)

$H(\boldsymbol{Y}_{1:T})=log(T-1)-\psi(k)+log(V_{d})+\frac{d}{T}\sum^{T}_{t=1}log(\rho_{k}(t))$

with $T$ standing for the number of samples, $\rho_{k}(t)$ \- the Euclidean
distance of the $k^{th}$ nearest neighbour of $\boldsymbol{y}_{t}$ in the
sample $\boldsymbol{Y}_{1:T}\backslash\\{\boldsymbol{y}_{t}\\}$ and
$V\subseteq\mathbb{R}^{d}$ a finite set.

  8. 8.

SH_KNN_2: Shannon differential entropy estimator using kNNs with k=3 and kd-
tree for quick nearest-neighbour lookup

  9. 9.

SH_KNN_3: Shannon differential entropy estimator using kNNs with k=5

  10. 10.

SH_MAXENT1: Maximum entropy distribution-based Shannon entropy estimator

$H(\boldsymbol{Y}_{1:T})=H(n)-\left[k_{1}\left(\frac{1}{T}\sum^{T}_{t=1}G_{1}(y^{\prime}_{t})\right)^{2}+k_{2}\left(\frac{1}{T}\sum^{T}_{t=1}G_{2}(y^{\prime}_{t})-\sqrt{\frac{2}{\pi}}\right)^{2}\right]+log(\hat{\sigma}),$

with

$\hat{\sigma}=\hat{\sigma}(\boldsymbol{Y}_{1:T})=\sqrt{\frac{1}{T-1}\sum^{T}_{t=1}(y_{t})^{2}},$
$y^{\prime}_{t}=\frac{y_{t}}{\hat{\sigma}},(t=1,\dots,T)$
$G_{1}(z)=ze^{\frac{-z^{2}}{2}},$ $G_{2}(z)=|z|,$
$k_{1}=\frac{36}{8\sqrt{3}-9},$ $k_{2}=\frac{1}{2-\frac{6}{\pi}},$

  11. 11.

SH_MAXENT2: Maximum entropy distribution-based Shannon entropy estimator, same
as SH_MAXENT1 with the following changes:

$G_{2}(z)=e^{\frac{-z^{2}}{2}},$ $k_{2}=\frac{24}{16\sqrt{3}-27},$

  12. 12.

SH_SPACING_V: Shannon entropy estimator using Vasicek’s spacing method.

$H(\boldsymbol{Y}_{1:T})=\frac{1}{T}\sum^{T}_{t=1}log\left(\frac{T}{2m}[y_{(t+m)}-y_{(t-m)}]\right)$

with $T$ number of samples, the convention that $y_{(t)}:=y_{(1)}$ if $t<1$
and $y_{(t)}:=y_{(T)}$ if $t>T$ and $m=\lfloor\sqrt{T}\rfloor$.

####  3.1.2 Execution

First, 199 different $i$ factors have been generated:

$i\in\\{0.01,0.02,\dots,1.00\\}\cup\\{1,2,\dots,100\\}$

For each $i$ every linear and non-linear combination of distribution types
have been tested (18 in total). That is, we have the general structures
$Y=X+N_{y}$ and $Y=X^{3}+N_{y}$ where $X\text{ and }N_{y}$ are drawn from the
three different distributions, $\mathcal{N},\mathcal{U}\text{ or
}\mathcal{L}$.

$Y=X\sim\mathcal{N}+N_{y}\sim\mathcal{N},$
$Y=X\sim\mathcal{N}+N_{y}\sim\mathcal{U},$ $\vdots$
$Y=X\sim\mathcal{L}^{3}+N_{y}\sim\mathcal{U}.$

Note that $\mathcal{L}^{3}$ here signifies the non-linear case
$Y=X^{3}+N_{y}$.
Next, for each of the 18 combinations we perform 100 tests each time. For
every test we generate 1000 new samples for the distribution of $X$ and the
distribution of $N_{y}$.
Lastly, we simply count how many tests are successful in these 100 tests and
define this ratio as our accuracy measure.

####  3.1.3 Experimental Results

In the following figures the y-axis shows the accuracy
($\frac{\text{\\#successful tests}}{100}$) and the x-axis shows the range of
the $i$ factor. Each figure contains two subfigures, the left with
$i\in\\{0.01,0.02,\dots,1.00\\}$ and the right figure with
$i\in\\{1,2,\dots,100\\}$. If the values of the estimators shown in the plots
are close to $0.5$, then this means that in 50% of the tests the algorithm
decided the correct direction (and vice versa 50% the wrong direction) and
thus indicates unidentifiability. If plots are closer to accuracy $1$ then we
have very good/consistent identifiability. The legend for both subfigures is
placed on the right side of the right subfigure and for the decoupled
experiments the estimator names are suffixed with "_S", S(plit). However, in
the analysis text the suffix is dropped for better readability. The next
subsection describes each figure individually and the following subsection
thereafter provides a summary and draws conclusions. Note that the first group
of results are obtained with the assumption "$X\to Y\oplus Y\to X$" and using
decoupled estimation. Also the plots for DISTCOV (dark green) and DISTCORR
(medium purple) often overlap (more than in 90% of cases), resulting in a dark
purple line.

##### Individual Analysis

This section can be skipped as we provide a summary table in the next section
(Table 3.1).

Fig. 3.3 shows the only case where we never achieve identifiability. This is
the well known linear Gaussian structural causal model and has only recently
been tackled successfully in [27] and [26].

Fig. 3.3 shows the linear model with $Y=\mathcal{N}+\mathcal{U}$, that is
independent variable $X$ being distributed according to the normal
distribution and noise - according to uniform distribution. SH_SPACING_V
performs the best with the accuracy of 100% for $i\in[0.55;7]$. HSIC_IC and
HSIC_IC2 perform the worst here with only an accuracy above 90% for
$i\in[3;7]$. All other estimators perform mediocre with an accuracy $>$ 80%
for $i\in[0.5;7]$.

Figure 3.2: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}+\mathcal{N}$

Figure 3.3: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}+\mathcal{U}$

Fig. 3.5 shows the linear model $Y=\mathcal{N}+\mathcal{L}$. The best
estimators are SH_MAXENT1 and SH_MAXENT2 with accuracy above 90% for
$i\in[0.35;2]$. HSIC also performs good with accuracy over 90% for
$i\in[0.38;2]$. The worst estimators are the three Shannon differential
entropy estimators using kNNs which never remain consistently above 80%
accuracy. The remaining estimators lie within the range 90% $\pm$ 8% accuracy
with $i\in[0.3;3]$. Fig. 3.5 shows the non-linear model
$Y=\mathcal{N}^{3}+\mathcal{N}$. Here all estimators perform very good with
$i\in[0.4;25]$ having an accuracy of almost 100%. With $i<0.3$ most estimators
drop fast below 90% accuracy. With $i\in[20;100]$ all estimators remain above
90% accuracy, except for HSIC_IC, HSIC_IC2, SH_MAXENT1 and SH_MAXENT2 which
drop below 90% after $i=45$. DISTCOV, SH_SPACING_V and the three Shannon kNN
estimators remain close to 100% in $i\in[0.01;100]$.

Figure 3.4: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}+\mathcal{L}$

Figure 3.5: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{N}$

Fig. 3.7 shows the non-linear model $Y=\mathcal{N}^{3}+\mathcal{U}$. This
result is similar as the previous one. With $i\in[0.45;80]$ we have 90% or
higher accuracy for all estimators. DISTCOV, SH_SPACING_V and the three
Shannon kNN estimators remain close to 100% in $i\in[0.01;100]$. Fig. 3.7
shows the non-linear model $Y=\mathcal{N}^{3}+\mathcal{L}$. Here all
estimators perform very good with $i\in[0.4;30]$ having an accuracy close to
100%. With $i<0.25$ all estimator drop rapidly and with $i>30$ HSIC_IC,
HSIC_IC2, SH_MAXENT1 and SH_MAXENT2 drop below 90% accuracy for higher $i$.
All others remain over 90% accuracy while HSIC remains around 90%.

Figure 3.6: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{U}$

Figure 3.7: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{L}$

Fig. 3.9 shows the linear model $Y=\mathcal{U}+\mathcal{U}$. HSIC_IC and
HSIC_IC2 only reach accuracy above 90% around $i=3$. All other estimator
perform quite good with $i\in[0.3;5]$ while only HOEFFDING and DISTCOV drop
slightly below 90% accuracy for $i=1$. SH_SPACING_V has 100% accuracy for
$i\in[0.12;10]$ and has on the remaining $i$ factors better accuracy than all
other estimators. Fig. 3.9 shows the linear model $Y=\mathcal{U}+\mathcal{N}$.
Here all estimators differ stronger than in the previous experiments. First,
SH_SPACING_V is performing the best with 100% accuracy for $i\in[0.08;2]$.
With $i=1$ all other estimators remain above 90%, expect HOEFFDING ($\sim$88%)
and HSIC_IC and HSIC_IC2 (both $\sim 75\%$). After $i=2$ all estimators drop
drastically towards 50% accuracy except SH_SPACING_V which remains above 70%.
For $i\in[0.2;1]$ some estimators remain between 80% and 95% while HSIC is
above 95%. Only HSIC_IC and HSIC_IC2 perform worse than all other estimators.

Figure 3.8: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}+\mathcal{U}$

Figure 3.9: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}+\mathcal{N}$

Fig. 3.11 shows the linear model $Y=\mathcal{U}+\mathcal{L}$. For
$i\in[0.3;1]$ all estimators perform well with 90% or higher accuracy, except
HSIC_IC and HSIC_IC2 which remain above 90% accuracy only after $i=0.45$.
After $i=1$ each estimator drops drastically and all converge towards 50%
accuracy, except for SH_SPACING_V which remains with a mean of 70% accuracy
higher than other estimators. For $i\in[0.08;1]$ SH_SPACING_V also has
accuracy 100%. For $i<0.3$ all other estimators drop fast towards 50%. Fig.
3.11 shows the non-linear model $Y=\mathcal{U}^{3}+\mathcal{U}$. For
$i\in[0.09;1]$ all estimators (except HSIC_IC and HSIC_IC2) remain above 95%
accuracy, while SH_KNN, SH_KNN_2, and SH_SPACING_V continue to do so for
$i\in[1;100]$. DISTCOV, DISCORR and HOEFFDING remain between 80% and 90%. HSIC
and SH_MAXENT1 drop to $>60$% after $i=20$ and remain above 60% for $i<100$.
SH_MAXENT2, HSIC_IC and HSIC_IC2 drop to 50% for $i\in[20;100]$.

Figure 3.10: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}+\mathcal{L}$

Figure 3.11: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{U}$

Fig. 3.13 shows the non-linear model $Y=\mathcal{U}^{3}+\mathcal{N}$. All
estimators, except HSIC_IC and HSIC_IC2, remain above 95% for $i\in[0.05;1]$
and for $i\in[1;100]$ the estimators converge differently. All three Shannon
differential entropy measures with kNNs and SH_SPACING_V remain above 95%
accuracy. DISTCOV, DISCORR and HOEFFDING keep a mean of $\sim$85% accuracy.
HSIC and SH_MAXENT1 remain above 60% accuracy. SH_MAXENT2 is pretty much
unidentifiable and HSIC_IC and HSIC_IC2 are unidentifiable for
$i\in[0.01;100]$. Fig. 3.13 shows the non-linear model
$Y=\mathcal{U}^{3}+\mathcal{L}$. The behaviour of different estimators is
almost the same as for $Y=\mathcal{U}^{3}+\mathcal{N}$. The only differences
are that HSIC_IC performs slightly better for $i\in[0.2;1]$ and DISTCOV,
DISTCORR, HSIC and SH_MAXENT2 perform worse.

Figure 3.12: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{N}$

Figure 3.13: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{L}$

Fig. 3.15 shows the linear model $Y=\mathcal{L}+\mathcal{L}$. First, for
$i\in[0.4;2]$ SH_MAXENT1 and SH_MAXENT2 have accuracy close to 100%. Next, for
$i\in[0.4;1]$ HSIC, HSIC_IC, HSIC_IC2, SH_SPACING_V, DISTCOV and DISTCORR
remain above 90% accuracy. HOEFFDING, and the three Shannon kNN estimators
never reach an accuracy above 90%. After $i=1$ all estimators drop fast
towards 50%. Fig. 3.15 shows the linear model $Y=\mathcal{L}+\mathcal{N}$.
This has a similar pattern as the previous one. For $i\in[0.3;1]$ HSIC,
SH_MAXENT1 and SH_MAXENT2 have accuracy greater than 90%. SH_SPACING_V,
HSIC_IC, HSIC_IC2, DISTCOV and DISTCORR lie between 85% and 95% accuracy for
$i\in[0.4;1]$ and HOEFFDING remains between 80% and 90%. Again, the three
Shannon kNN estimators never reach an accuracy higher than 80%. After $i=1$
all estimators drop fairly fast towards unidentifiability. Note: in Fig. 3.15
and Fig. 3.15 SH_KNN and SH_KNN_2 are overlapping completely.

Figure 3.14: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}+\mathcal{L}$

Figure 3.15: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}+\mathcal{N}$

Fig. 3.17 shows the linear model $Y=\mathcal{L}+\mathcal{U}$. SH_SPACING_V
performs the best of all estimators and has an accuracy of 100% for
$i\in[0.5;5]$. All other estimators slowly climb towards good identifiability
and for $i\in[0.7;7]$ they remain above 90% accuracy. Afterwards, all other
estimators drop with similar pace towards unidentifiability. Fig. 3.17 shows
the non-linear model $Y=\mathcal{L}^{3}+\mathcal{L}$. This experiments shows
the best results of all. For $i\in[0.1;100]$ all estimators (except SH_MAXENT1
and SH_MAXENT2) have an accuracy 90% or higher SH_SPACING_V and the three
Shannon kNN estimators have an accuracy of 100% for $i\in[0.01;100]$. Only
SH_MAXENT1 and SH_MAXENT2 perform bad at the beginning but still have an
accuracy of 90% or higher for $i\in[0.35;100]$.

Figure 3.16: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}+\mathcal{U}$

Figure 3.17: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{L}$

Fig. 3.19 shows the non-linear model $Y=\mathcal{L}^{3}+\mathcal{N}$. This
shows similar results as Fig. 3.17. For $i\in[0.1;100]$ all estimators (except
SH_MAXENT1 and SH_MAXENT2) have an accuracy 90% or higher SH_SPACING_V and the
three Shannon kNN estimators have an accuracy of 100% for $i\in[0.01;100]$.
Only SH_MAXENT1 and SH_MAXENT2 perform bad at the beginning but still have an
accuracy of 90% or higher for $i\in[0.35;100]$. Fig. 3.19 shows the non-linear
model $Y=\mathcal{L}^{3}+\mathcal{U}$. Similar as the two previous results.
For $i\in[0.15;100]$ all estimators (except SH_MAXENT1 and SH_MAXENT2) have an
accuracy 90% or higher SH_SPACING_V, and the three Shannon kNN estimators have
an accuracy of 100% for $i\in[0.01;100]$. Only SH_MAXENT1 and SH_MAXENT2
perform bad at the beginning but still have an accuracy of 90% or higher for
$i\in[0.7;100]$. For $i\in[1;100]$ all estimators are very close to 100%
accuracy.

Figure 3.18: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{N}$

Figure 3.19: RESIT & different noise levels & decoupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{U}$

##### Summary and Conclusion

Table 3.1 contains an overview of the results. The rows represent each
structural equation and each column represents the estimator used. The values
in the cells describe on what range of $i$ the estimator can reach over 90%.
Estimators have some variance in the results and thus on some intervals they
fall below 90% accuracy. The limits in the cells were chosen as follows: the
lower limit designates where estimators reaches the first time 90% or higher,
and the upper limit designates the last time where it reaches 90% or higher.
In between, most of the time estimators remain above 90% or rarely fall below,
but not more than 10%. Empty cells mean that in that case the estimator never
reached 90% or higher.
As the results show, different noise levels do have an impact on the
identifiability performance in RESIT methods. In general the linear equation
models are more fragile in RESIT than the non-linear equation models because
non-linear relationships tend to break the symmetry between the variables
easier, see [14]. Furthermore, in all cases the test results themselves have a
standard deviation between 0.05 to 0.1 as one can see in the sharp wiggles in
the plots.

Looking now only at the best estimation function and assuming a strong
identifiability equal or greater than 90% accuracy then for all linear cases
factor $i$ is smaller than 10 (even smaller than 5 in some cases, e.g., Fig.
3.5) and bigger than 0.5. In other words, if $i\notin[0.5;5]$ then the
accuracy is below 90%. This looks different for the non-linear cases. For Fig.
3.7 e.g., accuracy is greater than 90% if $i\in[0.05;100]$. Accuracy is even
equal to 100% in the cases where $X\sim\mathcal{L}$ (Fig. 3.17 \- Fig. 3.19)
with $i\in[0.01;100]$.

Some estimators perform differently depending on the setup. For example, in
Fig. 3.11 \- Fig. 3.13, the three Shannon differential entropy estimators with
kNNs always perform above 90% accuracy for $i\in[0.01;100]$ and even with 100%
accuracy for $i\in[0.01;100]$ in Fig. 3.17 \- Fig. 3.19. At the same time, the
HSIC with incomplete cholesky (HSIC_IC and HSIC_IC2) has the worst
performance. There is even one case where HSIC_IC never achieves
identifiability (Fig. 3.13). Overall, SH_SPACING_V performs the best in almost
all cases, and is only outperformed by SH_MAXENT1 and SH_MAXENT2 in Fig. 3.5,
Fig. 3.15 and Fig. 3.15. Some independence tests lose some of the accuracy
while entropy estimators retain an accuracy over 90% (e.g., Fig. 3.11, Fig.
3.13 and Fig. 3.13).
Additionally, it is worth to mention that entropy estimators are less
computationally demanding than independence tests but can be quite sensitive
to discretization effects [30]. However, entropy estimators can only be used
with the prior assumption we made: there is only one causal direction and it
is present in the model. Section 3.3 contains experiments with the same setup
but without this prior assumption.

In the next section, we compare decoupled estimation and coupled estimation.

Equation | HSIC | HISC_IC | HSIC_IC2 | DISTCOV | DISTCORR | HOEFFDING | SH_KNN | SH_KNN_2 | SH_KNN_3 | SH_MAXENT1 | SH_MAXENT2 | SH_SPACING_V
---|---|---|---|---|---|---|---|---|---|---|---|---
$Y=\mathcal{N}+\mathcal{N}$ |  |  |  |  |  |  |  |  |  |  |  |
$Y=\mathcal{N}+\mathcal{U}$ | 0.6 - 4 | 3 - 7 | 3 - 7 | 1 | 1 |  |  |  | 0.85 - 4 | 0.65 - 5 | 0.4 - 8 | 0.2 - 22
$Y=\mathcal{N}+\mathcal{L}$ | 0.38 - 2 | 0.75 - 2 | 0.75 - 2 | 0.8 - 1 | 0.87 - 1 |  |  |  |  | 0.3 - 3 | 0.33 - 3 | 0.82 - 1
$Y=\mathcal{N}^{3}+\mathcal{N}$ | 0.09 - 68 | 0.1 - 40 | 0.1 - 40 | 0.07 - 100 | 0.07 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.2 - 45 | 0.28 - 45 | 0.01 - 100
$Y=\mathcal{N}^{3}+\mathcal{U}$ | 0.16 - 100 | 0.16 - 88 | 0.17 - 88 | 0.14 - 100 | 0.14 - 100 | 0.05 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.33 - 85 | 0.43 - 91 | 0.01 - 100
$Y=\mathcal{N}^{3}+\mathcal{L}$ | 0.11 - 70 | 0.1 - 35 | 0.1 - 35 | 0.08 - 88 | 0.05 - 88 | 0.04 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.18 - 67 | 0.25 - 65 | 0.01 - 100
$Y=\mathcal{U}+\mathcal{U}$ | 0.18 - 6 | 3 | 3 | 0.2 - 3 | 0.2 - 3 | 0.2 - 3 | 0.2 - 3 | 0.2 - 2 | 0.2 - 3 | 0.23 - 3 | 0.12 - 8 | 0.05 - 21
$Y=\mathcal{U}+\mathcal{N}$ | 0.21 - 1 |  |  | 0.25 - 1 | 0.25 - 1 | 1 | 0.26 - 1 | 0.26 - 1 | 0.2 - 1 | 0.2 - 1 | 0.1 - 4 | 0.04 - 9
$Y=\mathcal{U}+\mathcal{L}$ | 0.15 - 1 | 0.35 - 1 | 0.35 - 1 | 0.2 - 1 | 0.2 - 1 | 0.15 - 1 | 0.21 - 1 | 0.21 - 1 | 0.16 - 1 | 0.12 - 3 | 0.1 - 3 | 0.03 - 8
$Y=\mathcal{U}^{3}+\mathcal{U}$ | 0.05 - 6 |  |  | 0.03 - 20 | 0.03 - 20 | 0.03 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.05 - 7 | 0.09 - 5 | 0.01 - 100
$Y=\mathcal{U}^{3}+\mathcal{N}$ | 0.04 - 3 |  |  | 0.01 - 85 | 0.01 - 85 | 0.01 - 95 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.02 - 4 | 0.08 - 2 | 0.01 - 100
$Y=\mathcal{U}^{3}+\mathcal{L}$ | 0.05 - 3 | 0.7 - 1 | 0.7 - 1 | 0.02 - 5 | 0.02 - 5 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.03 - 4 | 0.06 - 3 | 0.01 - 100
$Y=\mathcal{L}+\mathcal{L}$ | 0.3 - 5 | 0.4 - 4 | 0.4 - 4 | 0.23 - 4 | 0.23 - 4 |  |  |  | 0.6 - 1 | 0.2 - 6 | 0.17 - 5 | 0.17 - 4
$Y=\mathcal{L}+\mathcal{N}$ | 0.32 - 3 | 0.36 - 3 | 0.36 - 3 | 0.33 - 2 | 0.33 - 2 | 0.57 - 1 |  |  |  | 0.21 - 4 | 0.21 - 4 | 0.49 - 4
$Y=\mathcal{L}+\mathcal{U}$ | 0.4 - 8 | 0.6 - 7 | 0.6 - 7 | 0.4 - 5 | 0.4 - 5 | 0.4 - 5 | 0.53 - 4 | 0.53 - 4 | 0.51 - 5 | 0.32 - 10 | 0.32 - 10 | 0.16 - 27
$Y=\mathcal{L}^{3}+\mathcal{L}$ | 0.05 - 100 | 0.06 - 100 | 0.1 - 100 | 0.03 - 100 | 0.03 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.26 - 100 | 0.28 - 100 | 0.01 - 100
$Y=\mathcal{L}^{3}+\mathcal{N}$ | 0.04 - 100 | 0.07 - 100 | 0.1 - 100 | 0.03 - 100 | 0.03 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.28 - 100 | 0.3 - 100 | 0.01 - 100
$Y=\mathcal{L}^{3}+\mathcal{U}$ | 0.07 - 100 | 0.12 - 100 | 0.15 - 100 | 0.04 - 100 | 0.04 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.53 - 100 | 0.43 - 100 | 0.01 - 100

Table 3.1: Summary Table for RESIT & different noise levels & Decoupled
estimation. The numbers reflect the ranges of noise that allow identifiability
with accuracy around 90%.

###  3.2 RESIT with Different Noise Levels and coupled estimation

In the previous section, Section 3.1, we discussed decoupled estimation when
data was split into training (80%) an testing (20%) sets. The reason one would
want to split data is of computational nature: smaller sets of samples allow
for faster computation of the estimates, notably the independence estimators,
but might decrease the accuracy of the algorithm. The same data can be
analyzed using the coupled estimation method when the entire set of 1000
samples is used for both the regression step and the estimation step, and in
this case we achieve higher accuracy (or overall performance of the
algorithm). Therefore, choosing between decoupled and coupled estimation is a
trade-off between identifiability performance and computation speed. In this
section we analyze the performance of RESIT for this setup. As the
experimental procedure and the setup stay the same as in Section 3.1, we
directly proceed to the discussion of results.

####  3.2.1 Results

The figures for the results here are to be interpreted in the same way as in
the section before (see Section 3.1.3).

##### Individual Analysis

Again, reading individual results is not necessary as we provide an summary
table in the summary section (Table 3.2).
Fig. 3.21 is the only case which remains unchanged performance wise
($Y=\mathcal{N}+\mathcal{N}$). Fig. 3.21 shows the linear model with
$Y=\mathcal{N}+\mathcal{U}$. All estimators reach now an accuracy close to
100% at some interval $i\in[0.8;5]$. The accuracy in the $[0.01;1]$ interval
climbs now faster for all estimators and drops more slowly for $i\in[1;100]$.

Figure 3.20: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}+\mathcal{N}$

Figure 3.21: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}+\mathcal{U}$

Fig. 3.23 shows the linear model $Y=\mathcal{N}+\mathcal{L}$. All estimators
except for the three Shannon kNN estimators now reach an accuracy close to
100% for $i\in[0.42;5]$. The three kNN estimators achieve better performance
than in the decoupled estimation but remain worse than the other estimators.
Fig. 3.23 shows the non-linear model $Y=\mathcal{N}^{3}+\mathcal{N}$. Here all
estimators perform very good with $i\in[0.15;50]$ having an accuracy close to
100%. With $i<0.1$ most estimators drop fast below 90% accuracy. For $i>50$
most estimators remain above 90% accuracy, except SH_MAXENT2, HSIC_IC and
HSIC_IC2 which still remain above 80% accuracy.

Figure 3.22: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}+\mathcal{L}$

Figure 3.23: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{N}$

Fig. 3.25 shows the non-linear model $Y=\mathcal{N}^{3}+\mathcal{U}$. This
result is similar as the previous one. With $i\in[0.2;100]$ we have accuracy
close to 100% for all estimators. SH_MAXENT1 and SH_MAXENT2 now perform almost
as good as other estimators for $i<0.4$. Fig. 3.25 shows the non-linear model
$Y=\mathcal{N}^{3}+\mathcal{L}$. For $i\in[0.17;40]$ all estimators are close
to or at 100% accuracy. SH_MAXENT1 and SH_MAXENT2 perform as good as other
estimators for $i<0.2$. and $i>60$. The accuracy for HSIC_IC and HSIC_IC2
still drops fairly far after $i>50$ compared to the rest.

Figure 3.24: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{U}$

Figure 3.25: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{N}^{3}+\mathcal{L}$

Fig. 3.27 shows the linear model $Y=\mathcal{U}+\mathcal{U}$. This shows some
interesting changes. In the decoupled estimation HSIC_IC and HSIC_IC2 never
reached a consistent range with accuracy above 90% but in the coupled
estimation they do now. All estimators have an accuracy of 100% for
$i\in[0.25;5]$. Interestingly, SH_SPACING_V does not converge to 50%. HSIC_IC
and HSIC_IC2 do clearly fall below 40% for $i>\sim 85$. (Note that 0% accuracy
means that the wrong direction ($Y\to X$) is consistently identified.) Fig.
3.27 shows the linear model $Y=\mathcal{U}+\mathcal{N}$. For $i\in[0.2;1]$ all
estimators are close to or at 100% accuracy. In the decoupled estimation only
two estimators were close to 100% accuracy in that interval. After $i=1$ most
estimators converge towards $\sim$ 60% accuracy. HSIC_IC and HSIC_IC2
completely drop to unidentifiability. SH_SPACING_V remains above 90% accuracy
for $i\in[1;100]$.

Figure 3.26: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}+\mathcal{U}$

Figure 3.27: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}+\mathcal{N}$

Fig. 3.29 shows the linear model $Y=\mathcal{U}+\mathcal{L}$. This case has
similar results as the previous one. For $i\in[0.1;3]$ all estimators are
close to or at 100% accuracy. SH_SPANNING_V remains again over 90% accuracy
for $i\in[1;100]$. All other converge towards unidentifiability while the
three Shannon kNN estimators converge towards 60% accuracy. Fig. 3.29 shows
the non-linear model $Y=\mathcal{U}^{3}+\mathcal{U}$. HSIC_IC and HSIC_IC2 now
perform better and are close to 100% accuracy for $i\in[0.3;5]$ while others
have highest accuracy for $i\in[0.05;5]$. For $i$ larger than 5, results look
similar to the decoupled estimation but overall the performance increased by
an additive $5-15$% for the various estimators. The exception again is HSIC_IC
and HSIC_IC2 which lie clearly below 50% accuracy now.

Figure 3.28: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}+\mathcal{L}$

Figure 3.29: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{U}$

Fig. 3.31 shows the non-linear model $Y=\mathcal{U}^{3}+\mathcal{N}$. In the
decoupled estimation HSIC_IC and HSIC_IC2 never reached identifiability. This
is now different for the coupled estimation: for $i\in[0.1;0.75]$ they are
both above 90% accuracy but begin to drop already after $i=0.75$ towards
unidentifiability. SH_SPACING_V and the three Shannon kNN estimators have
perfect accuracy for $i\in[0.01;100]$. All other estimators have their
accuracy increased additively by $5-20$% for $i$ larger than 1. Fig. 3.31
shows the non-linear model $Y=\mathcal{U}^{3}+\mathcal{L}$. This is similar to
the previous case. SH_SPACING_V and the three Shannon kNN estimators are at
100% accuracy for $i\in[0.01;100]$. All others reach 100% accuracy for
$i\in[0.25;1]$ and for $i>1$ their accuracy have increased by an additive
$5-20$%.

Figure 3.30: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{N}$

Figure 3.31: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{U}^{3}+\mathcal{L}$

Fig. 3.33 shows the linear model $Y=\mathcal{L}+\mathcal{L}$. For
$i\in[0.01;1]$ estimators have improved significantly. The three Shannon kNN
estimators still perform worse than other estimators, but they do reach now
accuracy above 90% for $i\in[0.45;1]$. All other estimators are close to or at
100% for $i\in[0.2;5]$. Fig. 3.33 shows the linear model
$Y=\mathcal{L}+\mathcal{N}$. This has a similar pattern as the previous one.
The three Shannon kNN estimators still perform worse than other estimators and
they barely reach 90% accuracy for $i=1$. All other estimators are now close
to or at 100% for $i\in[0.25;5]$.

Figure 3.32: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}+\mathcal{L}$

Figure 3.33: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}+\mathcal{N}$

Fig. 3.35 shows the linear model $Y=\mathcal{L}+\mathcal{U}$. While
SH_SPACING_V performed the best of all estimators in the decoupled estimation
most estimators are now close to the performance of SH_SPACING_V or even
better. For $i\in[0.4;5]$ all estimators have an accuracy close to or at 100%.
Fig. 3.35 shows the non-linear model $Y=\mathcal{L}^{3}+\mathcal{L}$. For
$i\in[0.2;100]$ all estimators are now close to or at 100% accuracy.

Figure 3.34: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}+\mathcal{U}$

Figure 3.35: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{L}$

Fig. 3.37 shows the non-linear model $Y=\mathcal{L}^{3}+\mathcal{N}$. Same as
Fig. 3.17. For $i\in[0.2;100]$ all estimators are now close to or at 100%
accuracy. Fig. 3.37 shows the non-linear model
$Y=\mathcal{L}^{3}+\mathcal{U}$. Similar as the two previous results.
SH_MAXENT1 and SH_MAXENT2 do perform better now but are still worse than all
other estimators which, for $i\in[0.1;100]$, are now close to or at 100%
accuracy.

Figure 3.36: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{N}$

Figure 3.37: RESIT & different noise levels & coupled estimation &
$Y=\mathcal{L}^{3}+\mathcal{U}$

##### Summary and Conclusion

Table 3.2 follows the same scheme as Table 3.1. Similar to the previous
section, different noise levels do have an impact on the identifiability
performance of RESIT methods. Overall, all estimators perform better than in
the decoupled estimation case. Some estimators show better performance in
particular structural causal models (e.g., SH_SPACING_V in Fig. 3.27). Others
do perform worse than the rest of the used estimators (e.g., HSIC_IC and
HSIC_IC2 in Fig. 3.21). For all non-linear equation models the SH_SPACING_V
and the three Shannon kNN estimators are always at or at least very close to
100% accuracy for $i\in[0.01;100]$. SH_SPACING_V also keeps its good
performance in the linear equation models. For the independence measures HSIC,
DISTCOV, DISTCORR and HOEFFDING perform quite similar and good overall.
Furthermore, for linear equation models if the standard deviation of the noise
term is smaller than the standard deviation of the independent variable ($X$)
then RESIT is strong for $i\in[\sim 0.2;1]$. In the opposite direction, when
the standard deviation of the noise term is larger, RESIT does not perform
that well. Often it is only robust up to a factor $i=10$ (e.g., Fig. 3.23).
Note again that this results are also based on the assumption that in our
bivariate structure only one and exact one direction of cause and effect is
present (namely $X\to Y$). Therefore, without this assumption we cannot
compare the estimates directly but rather need to compare the estimate to a
derived p-value given some value of alpha (e.g., commonly $\alpha=0.05$). The
next section describes experiments and results without this assumption.

Equation | HSIC | HISC_IC | HSIC_IC2 | DISTCOV | DISTCORR | HOEFFDING | SH_KNN | SH_KNN_2 | SH_KNN_3 | SH_MAXENT1 | SH_MAXENT2 | SH_SPACING_V
---|---|---|---|---|---|---|---|---|---|---|---|---
$Y=\mathcal{N}+\mathcal{N}$ |  |  |  |  |  |  |  |  |  |  |  |
$Y=\mathcal{N}+\mathcal{U}$ | 0.17 - 18 | 0.65 - 26 | 0.7 - 26 | 0.16 - 23 | 0.16 - 23 | 0.16 - 25 | 0.32 - 12 | 0.32 - 12 | 0.24 - 12 | 0.23 - 12 | 0.15 - 22 | 0.13 - 33
$Y=\mathcal{N}+\mathcal{L}$ | 0.13 - 8 | 0.31 - 7 | 0.33 - 7 | 0.13 - 7 | 0.13 - 8 | 0.13 - 8 | 0.76 - 1 | 0.76 - 1 | 0.51 - 1 | 0.12 - 10 | 0.13 - 7 | 0.17 - 5
$Y=\mathcal{N}^{3}+\mathcal{N}$ | 0.04 - 100 | 0.04 - 83 | 0.08 - 83 | 0.02 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.05 - 100 | 0.11 - 98 | 0.01 - 100
$Y=\mathcal{N}^{3}+\mathcal{U}$ | 0.08 - 100 | 0.06 - 100 | 0.08 - 100 | 0.02 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.06 - 100 | 0.16 - 100 | 0.01 - 100
$Y=\mathcal{N}^{3}+\mathcal{L}$ | 0.04 - 100 | 0.04 - 70 | 0.09 - 70 | 0.02 - 100 | 0.02 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.05 - 100 | 0.1 - 100 | 0.01 - 100
$Y=\mathcal{U}+\mathcal{U}$ | 0.06 - 16 | 0.06 - 15 | 0.14 - 15 | 0.05 - 21 | 0.05 - 21 | 0.05 - 21 | 0.07 - 12 | 0.07 - 12 | 0.07 - 14 | 0.1 - 12 | 0.05 - 17 | 0.03 - 40
$Y=\mathcal{U}+\mathcal{N}$ | 0.05 - 6 | 0.04 - 3 | 0.1 - 3 | 0.04 - 7 | 0.04 - 7 | 0.04 - 7 | 0.08 - 4 | 0.08 - 4 | 0.05 - 5 | 0.06 - 4 | 0.03 - 7 | 0.01 - 100
$Y=\mathcal{U}+\mathcal{L}$ | 0.04 - 7 | 0.04 - 5 | 0.11 - 5 | 0.04 - 10 | 0.04 - 10 | 0.04 - 8 | 0.09 - 4 | 0.09 - 4 | 0.05 - 5 | 0.04 - 8 | 0.04 - 8 | 0.01 - 100
$Y=\mathcal{U}^{3}+\mathcal{U}$ | 0.03 - 16 | 0.14 - 13 | 0.17 - 13 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.02 - 90 | 0.04 - 12 | 0.01 - 100
$Y=\mathcal{U}^{3}+\mathcal{N}$ | 0.02 - 6 | 0.1 - 0.92 | 0.12 - 0.91 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.03 - 4 | 0.01 - 100
$Y=\mathcal{U}^{3}+\mathcal{L}$ | 0.03 - 7 | 0.1 - 4 | 0.17 - 4 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 88 | 0.04 - 5 | 0.01 - 100
$Y=\mathcal{L}+\mathcal{L}$ | 0.1 - 13 | 0.1 - 8 | 0.12 - 8 | 0.08 - 15 | 0.08 - 15 | 0.1 - 10 | 0.37 - 3 | 0.37 - 3 | 0.32 - 4 | 0.07 - 17 | 0.1 - 13 | 0.09 - 13
$Y=\mathcal{L}+\mathcal{N}$ | 0.1 - 7 | 0.1 - 4 | 0.1 - 4 | 0.1 - 7 | 0.1 - 7 | 0.1 - 7 | 0.61 - 1 | 0.61 - 1 | 0.37 - 3 | 0.07 - 13 | 0.1 - 7 | 0.14 - 6
$Y=\mathcal{L}+\mathcal{U}$ | 0.12 - 23 | 0.14 - 26 | 0.14 - 26 | 0.1 - 25 | 0.1 - 25 | 0.1 - 25 | 0.27 - 12 | 0.27 - 12 | 0.21 - 15 | 0.11 - 24 | 0.11 - 23 | 0.11 - 33
$Y=\mathcal{L}^{3}+\mathcal{L}$ | 0.02 - 100 | 0.03 - 100 | 0.09 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.1 - 100 | 0.15 - 100 | 0.01 - 100
$Y=\mathcal{L}^{3}+\mathcal{N}$ | 0.02 - 100 | 0.03 - 100 | 0.7 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.1 - 100 | 0.14 - 100 | 0.01 - 100
$Y=\mathcal{L}^{3}+\mathcal{U}$ | 0.04 - 100 | 0.05 - 100 | 0.07 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.01 - 100 | 0.17 - 100 | 0.21 - 100 | 0.01 - 100

Table 3.2: Summary Table for RESIT & different noise levels & Coupled
estimation. The numbers reflect the ranges of noise that allow identifiability
with accuracy around 90%.

###  3.3 RESIT with Different Noise Levels without prior assumption

Without assuming that exactly one direction has to be present in our bivariate
model, we must then account for four different possibilities instead:

  1. 1.

$X\to Y$ and $Y\cancel{\to}X$

  2. 2.

$X\cancel{\to}Y$ and $Y\to X$

  3. 3.

$X\to Y$ and $Y\to X$

  4. 4.

$X\cancel{\to}Y$ and $Y\cancel{\to}X$

This means that we cannot compare estimates directly but have to estimate an
$\alpha$ value for the independence test when testing each direction and we
can no longer use entropy estimators.

####  3.3.1 Setup

This setup is the same as in Section 3.1.1. The only difference is a small
modification in the algorithm we use (Algorithm 1). More precisely, the first
steps remain unchanged but the step number 5 needs to be adjusted. Here, we
perform independence test on the independent variable and the residuals and
compare the independence estimate to the alpha estimate with alpha = 0.05.
Furthermore, since the equation of our data is $Y=X+N_{y}$ only the outcome
$X\to Y$ and $Y\cancel{\to}X$ from our independence tests is the correct one
and thus we can formulate the following output function as the second step in
our modified algorithm:

$\displaystyle dir=\begin{cases}Correct&\text{if
}I(X_{Test},Y_{res})<\alpha(0.05)\land I(Y_{Test},X_{res})>\alpha(0.05),\\\
Incorrect&\text{otherwise }\end{cases}$ (3.3)

As before, Linear Regression is used. However, for simplicity we only used one
independence test: Hilbert-Schmidt Independence Criterion with RBF Kernel
(HSIC) 333Source: https://github.com/amber0309/HSIC. Here, we use the gamma
test to compare the HSIC estimate, that is with the alpha quantile of the
gamma distribution with mean and variance of HSIC under independence
hypothesis. The HSIC estimator [33] is chosen because of its general good
performance in previous results.

####  3.3.2 Execution

The execution is exactly the same as in Section 3.1.2.

####  3.3.3 Results

In the following figures the y-axis shows the accuracy
($\frac{\text{\\#successful tests}}{100}$) and the x-axis shows the range of
the $i$ factor. Each figure contains two subfigures, the left with
$i\in\\{0.01,0.02,\dots,1.00\\}$ and the right figure with
$i\in\\{1,2,\dots,100\\}$.. Differently from Section 3.1, if HSIC is closer to
0 then we have unidentifiability. If plots are closer to accuracy $1$ then we
have very good/consistent identifiability. The next subsection describes each
figure individually and the subsection thereafter provides a summary and a
small conclusion. Additionally, we also directly cover the decoupled
estimation together with the coupled estimation approach. Lastly, the figures
apply some abbreviations which we shortly explain here:

  * •

GAU: Gaussian distribution.

  * •

UNI: Uniform distribution.

  * •

LAP: Laplace distribution.

  * •

NL***: Non-linear variant.

Furthermore, if in a legend of a plot there are two of the above abbreviations
concatenated with a "+", e.g., GAU+UNI, then this signifies that
$X\sim\mathcal{N}$ and $N_{y}\sim\mathcal{U}$ (respecting order). If an
abbreviation stands alone, for example GAU, then both variables ($X,N_{y}$)
are drawn from that distribution.

##### Individual Analysis

Again, individual analysis can be skipped as we provide a summary in Table
3.3. The following paragraphs describe Fig. 3.39 \- Fig. 3.41 which show the
performance of using decoupled estimation with a split of 80%/20%. Fig. 3.39
shows all cases where both $X$ and $N_{y}$ are drawn from the same
distribution. $Y=\mathcal{N}+\mathcal{N}$ is the only case which never
achieves identifiability with $i\in[0.01;100]$. The non-linear cases are
pretty robust for $i\in[0.2;1]$; NLLAP remains above 90% accuracy for
$i\in[0.3;100]$ and NLGAU remains over 90% accuracy for $i\in[0.25;35]$.
However, NLUNI already drops fast for $i>2$ but remains above 90% accuracy for
$i\in[0.12;2]$. UNI and LAP are over 80% accuracy for $i\in[0.5;2]$ but drop
then fast after $i=2$. Fig. 3.39 shows all cases where $X\sim\mathcal{N}$ and
$N_{y}\cancel{\sim}\mathcal{N}$. GAU+UNI only reaches an accuracy of 90% when
$i=2$. GAU+LAP never reaches over 90% accuracy but remains partially over 80%
for $i\in[0.6;1]$. NL_GAU+UNI has an accuracy over 90% for $i\in[0.3;65]$ and
NL_GAU+LAP has also a good accuracy for $i\in[0.2;20]$.

Figure 3.38: $Y=X+N_{y}$. Contains all cases where $X$ and $N_{y}$ are drawn
from the same type of distribution. Dashed lines are non-linear cases.
Decoupled estimation (80% split).

Figure 3.39: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{N}$ and
$N_{y}\cancel{\sim}\mathcal{N}$. Dashed lines are non-linear cases. Decoupled
estimation (80% split).

Fig. 3.41 shows all cases where $X\sim\mathcal{U}$ and
$N_{y}\cancel{\sim}\mathcal{U}$. All cases drop fast to $<$ 10% accuracy for
$i>1$. NL_UNI+GAU and NL_UNI+LAP are over 90% accuracy for $i\in[0.1;1]$.
UNI+LAP hangs around 90% accuracy for $i\in[0.2;1]$ and UNI+GAU never reaches
a consistent accuracy above 90%, only over 80% for $i\in[0.33;0.9]$. Fig. 3.41
shows all cases where $X\sim\mathcal{L}$ and $N_{y}\cancel{\sim}\mathcal{L}$.
NL_LAP+GAU performs the best here with accuracy above 90% for $i\in[0.3;100]$,
followed by NL_LAP+UNI for $i\in[0.5;100]$. LAP_UNI only reaches 90% accuracy
around $i=1$. LAP+GAU never $>$ 90%, only 80% for $i\in[0.95;1]$.

Figure 3.40: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{U}$ and
$N_{y}\cancel{\sim}\mathcal{U}$. Dashed lines are non-linear cases. Decoupled
estimation (80% split).

Figure 3.41: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{L}$ and
$N_{y}\cancel{\sim}\mathcal{L}$. Dashed lines are non-linear cases. Decoupled
estimation (80% split).

The next paragraphs describe Fig. 3.43 \- Fig. 3.45 which show the performance
of using coupled estimation. Only differences to the decoupled estimation
counter part will be described. Fig. 3.43 shows all cases where both $X$ and
$N_{y}$ are drawn from the same distribution. For $i<1$ all cases (except
Gaussian noise only) reach over 90% accuracy earlier, $i\in[0.08;1]$ for the
non-linear equation models and $i\in[0.2;1]$ for the UNI and LAP models.
Furthermore, NLGAU drops more slowly for higher $i$. Fig. 3.43 shows all cases
where $X\sim\mathcal{N}$ and $N_{y}\cancel{\sim}\mathcal{N}$. For $i<1$ all
cases reach over 90% accuracy earlier. Also both linear models now reach an
accuracy over 90%, for $i\in[0.25;5]$ (GAU+LAP) and for $i\in[0.4;7]$
(GAU+UNI). Both non-linear equation models drop more slowly now for higher
$i$.

Figure 3.42: $Y=X+N_{y}$. Contains all cases where $X$ and $N_{y}$ are drawn
from the same type of distribution. Dashed lines are non-linear cases. Coupled
estimation.

Figure 3.43: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{N}$ and
$N_{y}\cancel{\sim}\mathcal{N}$. Dashed lines are non-linear cases. Coupled
estimation.

Fig. 3.45 shows all cases where $X\sim\mathcal{U}$ and
$N_{y}\cancel{\sim}\mathcal{U}$. For $i<1$ all cases reach over 90% accuracy
earlier. However, both linear equation models show the strongest improvement
and have now a consistent accuracy over 90% for $i\in[0.1;1]$. Fig. 3.45 shows
all cases where $X\sim\mathcal{L}$ and $N_{y}\cancel{\sim}\mathcal{L}$. For
$i<1$ all cases reach over 90% accuracy earlier. Again, both linear equation
models show the strongest improvement and have now a consistent accuracy over
90% for $i\in[0.25;1]$.

Figure 3.44: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{U}$ and
$N_{y}\cancel{\sim}\mathcal{U}$. Dashed lines are non-linear cases. Coupled
estimation.

Figure 3.45: $Y=X+N_{y}$. Contains all cases where $X\sim\mathcal{L}$ and
$N_{y}\cancel{\sim}\mathcal{L}$. Dashed lines are non-linear cases. Coupled
estimation.

##### Summary and Conclusion

Table 3.3 follows the same scheme as Table 3.1. However, columns designate
decoupled and coupled estimation instead. In general, non-linear cases are
often better identifiable than the linear cases. Our results show the
differences quite well and further confirm that "nonlinearities in the data-
generating process are in fact a blessing rather than a curse" \- [14].
Without the prior assumption models can be well identifiable but are also
affected by different noise levels. Also note that in the decoupled estimation
we used only 200 samples in the independence test, whereas in the coupled
estimation we used the entire set consisting of 1000 samples. For independence
tests this can make a difference when considering computation time.
Additionally, the changes in terms of improvement between decoupled and
coupled estimations are bigger in linear equation models than in non-linear
structural causal models. Therefore, for linear structural equations we
generally advise against decoupled estimation, and for non-linear structural
equations decoupled estimation can almost always be considered. Again, the
biggest advantage of decoupled estimation is lesser computation time, and if
the difference of computation time between decoupled and coupled estimation is
not significant, then we suggest to go with coupled estimation.

Combinations | Decoupled | Coupled
---|---|---
GAU |  |
UNI | 0.38 - 3 | 0.16 - 6
LAP | 0.57 - 2 | 0.23 - 4
NLGAU | 0.21 - 37 | 0.09 - 74
NLUNI | 0.14 - 4 | 0.05 - 7
NLLAP | 0.13 - 100 | 0.05 - 100
GAU+UNI |  | 0.39 - 7
GAU+LAP |  | 0.24 - 4
NL_GAU+UNI | 0.3 - 74 | 0.15 - 100
NL_GAU+LAP | 0.28 - 33 | 0.12 - 60
UNI+GAU | 0.42 - 0.84 | 0.13 - 3
UNI+LAP | 0.24 - 1 | 0.07 - 4
NL_UNI+GAU | 0.08 - 1 | 0.04 - 4
NL_UNI+LAP | 0.09 - 1 | 0.05 - 3
LAP+GAU |  | 0.23 - 4
LAP+UNI | 0.64 - 4 | 0.24 - 12
NL_LAP+GAU | 0.12 - 100 | 0.05 - 100
NL_LAP+UNI | 0.15 - 100 | 0.07 - 100
Table 3.3: Summary Table for RESIT & different noise levels & Coupled
estimation. The numbers reflect the ranges of noise that allow identifiability
with accuracy around 90%.

###  3.4 RESIT with different means

One question which emerged was how different means for the variables in the
equation of additive noise models have an impact on the outcome of the RESIT
method. In the previous experiments (Section 3.1 and Section 3.3), we always
used linear regression which is computational wise very cheap and showed very
good results in linear and non-linear structural causal models. Many recent
scientific papers state that any regression method can be used and explicitly
or implicitly indicate that means of 0 are assumed in the data (both
artificial and natural data). Although, from a mathematical point of view one
can already guess that for linear equations there should be no difference and
for the non-linear equations identifiability might only be observable around 0
mean for the $X$ variable. Nonetheless, we conducted a range of experiments to
analyze how different means for the noise term have an impact on the outcome
of RESIT. Again, we have two scenarios, the first one with the prior
assumption we made in Section 3.1 (that is only one causal direction must be
present in the bivariate case) and the second scenario without this
assumption.

####  3.4.1 Setup

The setup here is a bit different than in Section 3.1.1. For all empirical
tests we generate data following $X\to Y$ and generate linear as well as non-
linear cases for each test. Thus, we have the following structural causal
model:

$Y=\beta X+N_{y}$ $(Y=\beta X^{3}+N_{y}\text{ for the non-linear case})$

where

$\beta=1,$

and

$X\sim\begin{cases}\mathcal{N}(\mu_{X},1)&\text{or}\\\
\mathcal{U}(\mu_{X}-1,\mu_{X}+1)&\text{or}\\\
\mathcal{L}(\mu_{X},1)\end{cases}$

and

$N_{y}\sim\begin{cases}\mathcal{N}(\mu_{N},1)&\text{or}\\\
\mathcal{U}(\mu_{N}-1,\mu_{N}+1)&\text{or}\\\
\mathcal{L}(\mu_{N},1)\end{cases}$

where $\mu_{X},\mu_{N}\in\\{-100,-90,-80,\dots,100\\}$. Note that standard
deviation (range between lower and upper bound in the uniform case) for both
variables are the same. This equals the case where $i=1$ which showed one of
the best results in Section 3.1 and in Section 3.3 so we focus only on these
cases for this experiment setup. The algorithm used is the same as in Section
3.1 and we also make use of linear regression and the 12 estimators introduced
under Section 3.1.1.

####  3.4.2 Execution

As previously mentioned, the means $\mu_{X}$ and $\mu_{N}$ are drawn from a
set of 21 values:

$\mu_{X},\mu_{N}\in\\{-100,-90,-80,\dots,100\\}$

Thus we have 441 different combinations and for each combination a linear and
a non-linear model of different combinations of distribution types have been
tested (18 in total, these are again the same as in Section 3.1.2).

####  3.4.3 Results

In the following figures the y-axis shows the accuracy
($\frac{\text{\\#successful tests}}{100}$) and the x-axis shows all possible
combinations of the means for the $X$ and $N_{y}$ variable, in the following
form: $\mu_{X}/\mu_{N}$. If plots of the estimators are close to $0.5$
accuracy then this means that in 50% of the tests the algorithm decided the
correct direction (and vice versa 50% the wrong direction) and thus indicates
unidentifiability. If plots are closer to accuracy $1$ then we have very
good/consistent identifiability. This time we will not cover each individual
model because our results confirmed that there is no change at all in terms of
accuracy, and results are showing that if $\mu_{X}$ is significantly different
from 0 then for non-linear cases the models becomes unidentifiable.

##### Individual analysis

Fig. 3.47 contains two linear models: on the left we have
$Y=\mathcal{N}+\mathcal{U}$ and on the right $Y=\mathcal{U}+\mathcal{L}$. Now
if we compare them to Fig. 3.3 and Fig. 3.11, for $i=1$, respectively, we can
see that all estimators remain at the same accuracy. Fig. 3.47 contains two
non-linear models: on the left we have $Y=\mathcal{U}^{3}+\mathcal{U}$ and on
the right $Y=\mathcal{U}^{3}+\mathcal{N}$. For both models if
$\mu_{X}=0\land\mu_{N}\in[-100;100]$ then we have identifiability and
otherwise the models are unidentifiable as expected. However, in the left
model HSIC_IC only reaches $\sim$ 80% accuracy and for the right model HSIC_IC
does not achieve identifiability at all, same as in Fig. 3.11 and Fig. 3.13,
respectively.

Figure 3.46: Left: $X\sim\mathcal{N},N_{y}\sim\mathcal{U}$. Right:
$X\sim\mathcal{U},N_{y}\sim\mathcal{L}$.

Figure 3.47: Left: $X\sim\mathcal{U}^{3},N_{y}\sim\mathcal{U}$. Right:
$X\sim\mathcal{U}^{3},N_{y}\sim\mathcal{N}$.

##### Summary and Conclusion

As expected the causal relationship directions are recoverable for the linear
case on the entire test range for the means. However, for the non-linear case
the relationship directions are only recoverable for $\mu_{X}=0$. Furthermore,
for the non-linear case the estimators HSIC_IC and HSIC_IC2 perform badly,
similar to the previous results (see Tables 3.1, 3.2 and 3.3).

##  Chapter 4 Experiments on Identification using Conditional Variances

As already mentioned in Section 3.1.2 the problem of the identifiability of
the linear Gaussian structural equation model has only been recently tackled
successfully. [23], [34] and [35] proved that Gaussian linear structural
equation models with equal or known error variances are identifiable, and [36]
and [27] proved that Gaussian linear structural equation models with unknown
heterogeneous error variances can be identifiable. In this chapter, we will
discuss and test the algorithm of [26] and [29] on bivariate additive noise
models, which proved the identifiability of Gaussian structural equation
models with both homogeneous and heterogeneous unknown error variances.

###  4.1 Introduction

Let

$G_{1}:Y=\beta X+N_{Y}$ $G_{2}:X=\beta Y+N_{X}$ $G_{3}:X=N_{X}\text{ and
}Y=N_{Y}$

be the three possible structural equation models for any bivariate additive
noise models (acyclic). For $G_{1}$ [29] states that if the error variance
ratio satisfies $\frac{Var(Y)}{Var(X)}>(1-\beta^{2})$ then from the law of
total variance we have the following two conditions:

$(A)\text{
}Var(Y)=\mathbb{E}(Var(Y|X))+Var(\mathbb{E}(Y|X))=Var(Y)+\beta^{2}Var(X)>Var(X)$

$\displaystyle(B)\text{ }\mathbb{E}(Var(X|Y))$
$\displaystyle=Var(X)-Var(\mathbb{E}(X|Y))$
$\displaystyle=Var(X)-\frac{\beta^{2}Var(X)^{2}}{\beta^{2}Var(X)+Var(Y)}<Var(Y)=\mathbb{E}(var(Y|X))$

For $(A)$ the intuition here is that the level of uncertainty of $X$ is lower
than the level of uncertainty of $Y$ since $X$ has only one random source and
$Y$ has two random sources. For $(B)$ the intuition is that after eliminating
the other variable effect the level of uncertainty of $X$ is smaller than the
level of uncertainty of $Y$. This is because after eliminating the effect of
$X$ the remaining part of $Y$ is $N_{Y}$ and when eliminating the effect of
$Y$ on $X$ the remaining part is $N_{X}$ ($=X$) since $Y$ contains some
information of $N_{X}$. Therefore, one can order the observed variables and
one can always recover the true ordering as long as
$\frac{Var(Y)}{Var(X)}>(1-\beta^{2})$, even if the error variances are
different.
In the same fashion, one can find the true ordering for $G_{2}$ as long as
$\frac{Var(X)}{Var(Y)}>(1-\beta^{2})$. For the last structural equation model
$G_{3}$, we have no guarantee which marginal or conditional variance is bigger
and therefore any ordering is considered correct.
Lastly, after determining the orderings for each graph, we simply need to test
for the presence of an edge by verifying the dependence relationships between
variables. For the graphs $G_{1}$ and $G_{2}$, the variables $X$ and $Y$ are
dependent and for the graph $G_{3}$ the variables are independent. Thus one
can recover the true graphs.

####  4.1.1 Algorithm

The algorithm used for the experiments in this work is composed of Algorithm 2
and Algorithm 3 from [29]. As previously hinted, the algorithm is composed of
two parts: 1) ordering and 2) conditional independence testing.

##### 1) Ordering - Backward step wise selection

For the ordering step the paper [29] proposes two algorithms (Algorithm 1 and
2 from [29]):

Theorem 2 (Identifiability Conditions for ANMs) Let $P(X)$ be generated from
an Additive Noise Model with directed acyclic graph G and true ordering $\pi$.
Suppose that causal minimality holds. Then, G is uniquely identifiable if
either of the two following conditions is satisfied: For any node
$j=\pi_{m}\in V,k\in De(j)$, and $l\in An(j)$,

$(A)\text{ Forward step wise selection:
}\sigma^{2}_{j}<\sigma^{2}_{k}+\mathbb{E}(Var(\mathbb{E}(X_{k}|X_{Pa(k)})|X_{\pi_{1}},\dots,X_{\pi_{m-1}}))$

, or

$(B)\text{ Backward step wise selection:
}\sigma^{2}_{j}>\sigma^{2}_{l}-\mathbb{E}(Var(\mathbb{E}(X_{l}|X_{\pi_{1}},\dots,X_{\pi_{m}}\text{\textbackslash}X_{l})|X_{Pa(l)}))$

For $(A)$ an additive noise model is identifiable if the conditional variance
of a node $j$ is smaller than that of its descendant, De($j$), given the non-
descendants, Nd($j$). This can be understood that the variance of $N_{j}$ is
overestimated owing to lack of parents. For $(B)$ an additive noise model is
identifiable if the conditional variance of a node $j$ given its parents,
Pa($j$), is bigger than that of its ancestor, An($j$), given the union of its
parents and any of its descendants.

In this work we used Backward step wise selection (Algorithm 2) as it was more
convenient to implement in Python. First we have our set $S$ which contains
all nodes from the additive noise model. We iterate over this set $S$ and for
each node we calculate its conditional variance given all other remaining
nodes in the set $S$. We select the node with the highest conditional variance
and append it to the ordering $\pi$ and also remove it from the set $S$. With
the updated set $S$ we repeat the process of finding the node with the highest
conditional variance, append it to the ordering $\pi$ and remove it from $S$.
This is repeated until $S$ is empty. Lastly, the reverse of the ordering $\pi$
is returned. (The first node to be appended to the ordering is actually the
last one in the ordering, therefore "backward step wise selection" is used).

##### 2) Uncertainty Scoring

The second step consists of the uncertainty scoring (Algorithm 3). Here, we
iterate over the ordering $\pi$ and for each node $j$ we perform conditional
independence tests with each other node $l$ appearing before node $j$ in the
ordering $\pi$. If any node $l$ is dependent on $j$, then we add node $l$ to
Pa($j$). Here, the first node in the ordering never has a parent and we start
the iteration at the second node in the ordering. For the conditional
independence test, Fisher’s z-transform of the partial correlation is used. In
the following sections we refer to this algorithm as Park algorithm or simply
Park.

Algorithm 2 Backward step-wise selection

1:Input: All variables from an ANM: $X=(x_{1},x_{2},\dots,x_{n})$

2:Output: Estimated ordering $\pi=(\pi_{1},\pi_{2},\dots,\pi_{n})$

3:

4:Procedure

5:Set $S=\\{1,2,\dots,n\\}$

6:List $\pi=[\hskip 5.69054pt]$

7:for $m=1\dots n\textbf{ do}$

8: for $j\in S$ do

9: Estimate the conditional variance $x_{j}$ given
$\\{x_{1},\dots,x_{n}\\}\textbackslash x_{j}$, $\sigma^{2}_{j|S\textbackslash
j}$

10: end

11: Append $\pi_{m}=argmax_{j}\sigma^{2}_{j|S\textbackslash j}$ to $\pi$

12: Update $S=S\textbackslash\pi_{m}$

13:end

14:Reverse list $\pi$

Algorithm 3 Uncertainty Scoring

1:Input: All variables from an ANM: $X=(x_{1},x_{2},\dots,x_{n})$

2:Output: Dictionary with estimated parents for all variables:
$G=\\{Pa(x_{1}):[\dots],Pa(x_{2}):[\dots],\dots,Pa(x_{n}):[\dots]\\}$

3:

4:Procedure

5:1) Get ordering from backward step-wise selection:

6:$\pi=(\pi_{1},\pi_{2},\dots,\pi_{n})$

7:

8:2) Parents estimation

9:$G=\\{\\}$

10:for $m=2\dots n\textbf{ do}$

11: $Pa(\pi_{m})=[\text{ }]$

12: for $j=1\dots m-1$ do

13: Conditional independence test between $\pi_{m}$ and $\pi_{j}$ given
$\\{\pi_{1},\dots,\pi_{m-1}\\}\textbackslash\pi_{j}$

14: If dependent, include $\pi_{j}$ into $Pa(\pi_{m})$

15: end

16: Insert $Pa(\pi_{m})$ into $G$

17:end

###  4.2 Setup

The setup is the same as in Section 3.1.1.

####  4.2.1 Execution

The execution is the same as in Section 3.1.2 except that we use the Park
algorithm instead. Here, when using the Park algorithm the output will be a
list of parent sets for all nodes in the additive noise model:

$\text{Output: }[(Pa(X),Pa(Y)].$

A test is then successful if and only if

$(Pa(X)=\\{\\})\land(Pa(Y)=\\{X\\}).$

####  4.2.2 Results

In the following figures the y-axis shows the accuracy
($\frac{\text{\\#successful tests}}{100}$) and the x-axis shows the range of
the $i$ factor. Each figure contains two subfigures, the left with
$i\in\\{0.01,0.02,\dots,1.00\\}$ and the right figure with
$i\in\\{1,2,\dots,100\\}$. Differently from Section 3.1, if HSIC is closer to
0 then we have unidentifiability. If plots are closer to accuracy $1$ then we
have very good/consistent identifiability. The next subsection describes each
figure individually and the following subsection thereafter provides a summary
and a small conclusion.

##### Individual Analysis

Fig. 4.2 shows the results for all linear cases. For $i\in[0.5;3]$ all linear
structural equation models are always identifiable. For $i<0.5$ all models
slowly drop towards 50% accuracy. For $i>3$ all cases start dropping in terms
of accuracy. Some do drop faster (UNIxLAP, UNIxGAU, GAUxLAP, GAU, UNI, LAP)
and others more slowly (LAPxGAU, GAUxUNI, LAPxUNI). For $i>5$ the first case
(UNIxLAP) is below 90% accuracy and for $i>25$ all cases are below 90%
accuracy.

Fig. 4.2 shows the results for non-linear cases. Here, the results are quite
different and varying. NL_UNIxLAP,NL_UNIxGAU and NLUNI have accuracy 100% for
$i\in[0.12;1]$. After $i=1$ they drop fast below 50% accuracy. NLGAU and
NL_GAUxLAP reach an accuracy over 90% for $i\in[0.35;25]$ and slowly drop
outside this range in terms of accuracy. NL_GAUxUNI has accuracy over 90% for
$i\in[0.77;58]$ and slowly drops outside this range. The last three, NLLAP,
NL_LAPxGAU and NL_LAPxUNI only reach an accuracy higher than 90% for
$i\in[3;100]$.

Figure 4.1: $Y=X+N_{y}$. Contains all linear cases.

Figure 4.2: $Y=X^{3}+N_{y}$. Contains all non-linear cases.

##### Summary and Conclusion

Table 4.1 follows the same scheme as Table 3.3. Interestingly, for
$i\in[0.01;1]$ the linear cases perform better than the non-linear cases and
for $i\in[1;100]$ several non-linear cases perform better than most linear
cases, which wasn’t the case for RESIT. Furthermore, as the results show,
after specific noise levels the identifiability from the Park algorithm starts
to suffer and thus not all graphs can be successfully recovered. The reason
for this lies within step 2, the conditional independence test. If noise
levels are significantly different then the independence test fails to capture
the correlation between the two nodes and therefore concludes that the nodes
are independent (Type II Error). However, for any given $i$, the ordering step
always performs correctly.111A quick test in python shell, with $i=57$,
$X\sim\mathcal{L}$ and $N_{y}\sim\mathcal{U}$ and 100 repetitions showed that
in these runs the ordering was always correct but only in 35 runs (from the
100 repetitions) the independence tests were correct.

Combinations | Coupled
---|---
GAU | 0.08 - 10
UNI | 0.1 - 8
LAP | 0.1 - 10
GAU+UNI | 0.16 - 10
GAU+LAP | 0.05 - 6
UNI+GAU | 0.04 - 5
UNI+LAP | 0.03 - 3
LAP+GAU | 0.14 - 13
LAP+UNI | 0.19 - 26
NLGAU | 0.33 - 37
NLUNI | 0.05 - 6
NLLAP | 2 - 100
NL_GAU+UNI | 0.52 - 67
NL_GAU+LAP | 0.23 - 25
NL_UNI+GAU | 0.04 - 4
NL_UNI+LAP | 0.03 - 3
NL_LAP+GAU | 4 - 100
NL_LAP+UNI | 5 - 100
Table 4.1: Summary Table for Uncertainty Scoring & different noise levels &
Coupled estimation. The numbers reflect the ranges of noise that allow
identifiability with accuracy around 90%.

##  Chapter 5 Conclusion and Future Work

The results from the experiments showed that two analyzed causal discovery
methods (RESIT and Uncertainty Scoring) are affected by different noise
scales. For significantly small noise levels in the disturbance term $N_{y}$
(almost deterministic data), or significantly high noise levels, these causal
discovery methods fail to capture the true causal relationship of the given
structural equation model. Recall that significantly here depends on the
model. For example, on some models if the noise level was already twice larger
then the methods failed to determine the causal direction consistently. Other
models remained identifiable with 100 times higher noise levels (The range of
different noise levels in the experiments was quite exhaustive, and
realistically speaking having noise levels 100 times higher then the potential
cause variable is very rare111Additionally, with very high noise levels the
effect of the cause variable is very likely negligible anyways). This also
shows different behavior for different distribution types (e.g., Gaussian or
Uniform). Therefore, if observed data differs significantly in terms of
variances scientists need to be careful when analyzing causal relationships
and in drawing conclusions.

In Chapter 3 we tested RESIT in the bi-variate case using several estimators
(12 in total). We compared decoupled estimation (splitting data into training
and test data) and coupled estimation (no splitting) in Section 3.1 and
Section 3.2. In Section 3.3 we omitted a particular assumption which allowed
us to use entropy estimators and therefore restricted our setting in using
only independence tests. All results from these sections are quite similar. In
general, if the standard deviation of the noise term is smaller than the
standard deviation of the cause then models remained identifiable as opposed
to the case when the standard deviation of the noise term was bigger. For
example, often when the standard deviation of the noise term was only half of
that of the cause, the model was still identifiable. However, in several
cases, if the standard deviation of the noise term was already twice larger
than the standard deviation of the cause, then the model became
unidentifiable. We also tested linear and non-linear models and our results
show that non-linear models were still identifiable in situations where a
linear model were unidentifiable. For example, some non-linear models, where
the noise term’s standard deviation was 100 times higher than that of the
cause, were still perfectly identifiable while their linear counterpart
wasn’t. Similar results can be observed in the Uncertainty Scoring method
(Chapter 4).
Specific to RESIT, we also tested whether splitting data into training and
test data was beneficial. The reason to do this is to use only the test data
in the estimators and therefore reduce the overall computation time of the
algorithm. Hence, we analysed the identifiability performance between
decoupled estimation and coupled estimation. Our results show that there is a
performance drop from coupled to decoupled estimation. However, in linear
cases, this drop is bigger than in non-linear cases. Also, our results show
that decoupled estimation for non-linear cases has almost no significant
impact on the performance of RESIT. For linear models we generally advise to
use coupled estimation. In the end this is a trade-off between performance and
computation speed and should be carefully selected depending on the individual
experiment. Lastly, for RESIT we used several estimators (6 independence
estimators and 6 entropy estimators). Our results show differences in terms of
performance in these estimators. In our experiments Hilbert-Schmidt
Independence Criterion with RBF Kernel was the best independence estimator,
and Shannon entropy estimator using Vasicek’s spacing method was the best
entropy estimator.

Future work. In our experiments, we tested only two particular methods and
three different distribution types. Similar results are expected with other
methods for causal discovery for additive noise models, as the failing point
is the independence estimation (or entropy estimation) in the two used
methods. Therefore, methods relying on these estimations are generally prone
to errors for some levels of noise. Furthermore, this work does not formalize
the effect of different noise levels in ANM causal discovery methods but could
be done in a future work.

In reality, observed data does not always strictly follow a certain
distribution type. As there are many different combinations possible, it would
be interesting to generalize the impact of different noise levels on any
distribution by using the different properties an observed distribution
exhibits.

### Acknowledgement

I want to express my many thanks to Dr. Marharyta Aleksandrova. I enjoyed
working on this topic a lot, and I am grateful to have had Dr. Aleksandrova as
my advisor. She is a very proactive and professional researcher, and I learned
a lot from her. Furthermore, she gave me a lot of feedback on my thesis text,
which helped me improve the quality of the thesis text a lot.

## References

  * [1] Sewall Wright  “Correlation and causation”  In _J. agric. Res._ 20, 1921, pp. 557–580
  * [2] George Rebane and Judea Pearl  “The Recovery of Causal Poly-Trees from Statistical Data”  In _CoRR_ abs/1304.2736, 2013  arXiv: http://arxiv.org/abs/1304.2736
  * [3] Peter Spirtes, Clark Glymour and Richard Scheines  “Causation, Prediction, and Search”  In _Causation, Prediction, and Search_ 81, 1993  DOI: 10.1007/978-1-4612-2748-9
  * [4] Peter Spirtes, Clark N Glymour, Richard Scheines and David Heckerman  “Causation, prediction, and search”  MIT press, 2000
  * [5] Pearl Judea  “Causality: models, reasoning, and inference”  In _Cambridge University Press. ISBN 0_ 521.77362, 2000, pp. 8
  * [6] Dominik Janzing et al.  “Information-geometric approach to inferring causal directions”  In _Artificial Intelligence_ 182 Elsevier, 2012, pp. 1–31
  * [7] Eleni Sgouritsa, Dominik Janzing, Philipp Hennig and Bernhard Schölkopf  “Inference of cause and effect with unsupervised inverse regression”  In _Artificial intelligence and statistics_ , 2015, pp. 847–855  PMLR
  * [8] Nir Friedman and Iftach Nachman  “Gaussian Process Networks”  In _CoRR_ abs/1301.3857, 2013  arXiv: http://arxiv.org/abs/1301.3857
  * [9] Yutaka Kano and Shohei Shimizu  “Causal inference using nonnormality”  In _Proceedings of the international symposium on science of modeling, the 30th anniversary of the information criterion_ , 2003, pp. 261–270
  * [10] Shohei Shimizu et al.  “A linear non-Gaussian acyclic model for causal discovery.”  In _Journal of Machine Learning Research_ 7.10, 2006
  * [11] Shohei Shimizu, Aapo Hyvarinen and Yoshinobu Kawahara  “A direct method for estimating a causal ordering in a linear non-Gaussian acyclic model”, 2014  arXiv:1408.2038 [cs.LG]
  * [12] Xiaohai Sun, Dominik Janzing and Bernhard Schölkopf  “Causal inference by choosing graphs with most plausible Markov kernels”  In _Ninth International Symposium on Artificial Intelligence and Mathematics (AIMath 2006)_ , 2006, pp. 1–11
  * [13] Xiaohai Sun, Dominik Janzing and Bernhard Schölkopf  “Causal reasoning by evaluating the complexity of conditional densities with kernel methods”  In _Neurocomputing_ 71.7-9 Elsevier, 2008, pp. 1248–1256
  * [14] Patrik Hoyer et al.  “Nonlinear causal discovery with additive noise models”  In _Advances in neural information processing systems_ 21, 2009, pp. 689–696
  * [15] Dominik Janzing, Patrik O Hoyer and Bernhard Schölkopf  “Telling cause from effect based on high-dimensional observations”  In _arXiv preprint arXiv:0909.4386_ , 2009
  * [16] Joris Mooij, Dominik Janzing, Jonas Peters and Bernhard Schölkopf  “Regression by dependence minimization and its application to causal inference in additive noise models”  In _Proceedings of the 26th annual international conference on machine learning_ , 2009, pp. 745–752
  * [17] Oliver Stegle et al.  “Probabilistic latent variable models for distinguishing between cause and effect”  In _Advances in neural information processing systems_ 23 Citeseer, 2010, pp. 1687–1695
  * [18] Joris M Mooij, Dominik Janzing, Tom Heskes and Bernhard Schölkopf  “On causal discovery with cyclic additive noise model”  [Sl]: NIPS, 2011
  * [19] Kun Zhang and Aapo Hyvarinen  “On the Identifiability of the Post-Nonlinear Causal Model”, 2012  arXiv:1205.2599 [stat.ML]
  * [20] Povilas Daniusis et al.  “Inferring deterministic causal relations”  In _arXiv preprint arXiv:1203.3475_ , 2012
  * [21] Aapo Hyvärinen and Stephen M Smith  “Pairwise likelihood ratios for estimation of non-Gaussian structural equation models”  In _Journal of Machine Learning Research_ 14.Jan, 2013, pp. 111–152
  * [22] Shohei Shimizu  “LiNGAM: Non-Gaussian methods for estimating causal structures”  In _Behaviormetrika_ 41.1 Springer, 2014, pp. 65–98
  * [23] J. Peters and P. Bühlmann  “Identifiability of Gaussian structural equation models with equal error variances”  In _Biometrika_ 101.1, 2013, pp. 219–228  DOI: 10.1093/biomet/ast043
  * [24] Jonas Peters, Joris M Mooij, Dominik Janzing and Bernhard Schölkopf  “Causal discovery with continuous additive noise models”, 2014
  * [25] Christopher Nowzohour and Peter Bühlmann  “Score-based causal learning in additive noise models”  In _Statistics_ 50.3 Taylor & Francis, 2016, pp. 471–485
  * [26] Gunwoong Park and Younghwan Kim  “Identifiability of Gaussian Structural Equation Models with Homogeneous and Heterogeneous Error Variances”, 2019  arXiv:1901.10134 [stat.ML]
  * [27] Wenyu Chen, Mathias Drton and Y Samuel Wang  “On causal discovery with an equal-variance assumption”  In _Biometrika_ 106.4 Oxford University Press, 2019, pp. 973–980
  * [28] Samory Kpotufe, Eleni Sgouritsa, Dominik Janzing and Bernhard Schölkopf  “Consistency of causal inference under the additive noise model”  In _International Conference on Machine Learning_ , 2014, pp. 478–486  PMLR
  * [29] Gunwoong Park  “Identifiability of Additive Noise Models Using Conditional Variances.”  In _Journal of Machine Learning Research_ 21.75, 2020, pp. 1–34
  * [30] Joris M Mooij et al.  “Distinguishing cause from effect using observational data: methods and benchmarks”  In _The Journal of Machine Learning Research_ 17.1 JMLR. org, 2016, pp. 1103–1204
  * [31] Thomas M Cover  “Elements of information theory”  John Wiley & Sons, 1999
  * [32] Zoltán Szabó  “Information Theoretical Estimators Toolbox”  In _Journal of Machine Learning Research_ 15, 2014, pp. 283–287
  * [33] Arthur Gretton, Olivier Bousquet, Alex Smola and Bernhard Schölkopf  “Measuring statistical dependence with Hilbert-Schmidt norms”  In _International conference on algorithmic learning theory_ , 2005, pp. 63–77  Springer
  * [34] Po-Ling Loh and Peter Bühlmann  “High-dimensional learning of linear causal networks via inverse covariance estimation”, 2013  arXiv:1311.3492 [stat.ML]
  * [35] Asish Ghoshal and Jean Honorio  “Learning identifiable gaussian bayesian networks in polynomial time and sample complexity”  In _arXiv preprint arXiv:1703.01196_ , 2017
  * [36] Asish Ghoshal and Jean Honorio  “Learning linear structural equation models in polynomial time and sample complexity”  In _International Conference on Artificial Intelligence and Statistics_ , 2018, pp. 1466–1475  PMLR