[1]Masatoshi Kobayashi

1] Graduate School of Information Science and Technology, The University of
Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan

2] IBM Research – Tokyo, Tokyo, Japan

# Detection of Unobserved Common Causes based on NML Code in Discrete, Mixed,
and Continuous Variables

<EMAIL_ADDRESS>Kohei Miyaguchi<EMAIL_ADDRESS>Shin Matsushima  https://orcid.org/0000-0002-8160-4310 [  [

###### Abstract

Causal discovery in the presence of unobserved common causes from
observational data only is a crucial but challenging problem. We categorize
all possible causal relationships between two random variables into the
following four categories and aim to identify one from observed data: two
cases in which either of the direct causality exists, a case that variables
are independent, and a case that variables are confounded by latent
confounders. Although existing methods have been proposed to tackle this
problem, they require unobserved variables to satisfy assumptions on the form
of their equation models. In our previous study [10], the first causal
discovery method without such assumptions is proposed for discrete data and
named CLOUD. Using Normalized Maximum Likelihood (NML) Code, CLOUD selects a
model that yields the minimum codelength of the observed data from a set of
model candidates. This paper extends CLOUD to apply for various data types
across discrete, mixed, and continuous. We not only performed theoretical
analysis to show the consistency of CLOUD in terms of the model selection, but
also demonstrated that CLOUD is more effective than existing methods in
inferring causal relationships by extensive experiments on both synthetic and
real-world data.

###### keywords:

Causal Discovery, Unobserved Common Causes, Discrete, Mixed, Continuous Data,
SCM, MDL Principle, Model Selection, NML Code

##  1 Introduction

Intelligent systems that utilize accurate prediction based on data has enjoyed
remarkable success by the development of machine learning methodology. It is
expected that an accurate predictor learned from data possesses a certain form
of information on the data. This expectation motivates us to render
information extracted from data by a learning algorithm into a form in which
humans can understand. For this sake, it is considered that the study of
causal inference, aiming at extracting the underlying causal mechanism, has
gained prominence in the context of machine learning as well as other various
fields [25, 31].

As it is impractical to perform randomized control trials in many cases, many
studies have focus on inferring causal structures based soly on observational
data [20]. In the context of classical causal discovery, the decision whether
$X$ is a direct cause of $Y$ or $Y$ is a direct cause of $X$ is already a hard
problem so that standard methods assume there is no unknown common cause
(causal sufficiency). However, in practice, this assumption is often violated,
which can results in the methods producing unreliable results. Therefore,
causal discovery that allows for the presence of unobserved common causes
becomes crucial.

We consider a new problem setting by revisiting Reichenbach’s common cause
principle quoted as follows [25, Principle 1.1]:

> If two random variables $X$ and $Y$ are statistically dependent, then there
> exists a third variable $C$ that causally influences both. (As a special
> case, $C$ may coincide with either $X$ or $Y$.) Furthermore, this variable
> $C$ screens $X$ and $Y$ from each other in the sense that given $C$, they
> become independent.

As a logical conclusion of the statement above, we can always categorize the
relationship between $X$ and $Y$ into 4 cases: 1) $X$ causally influences ($X$
and) $Y$, 2) $Y$ causally influences $X$ (and $Y$), 3) there exists a third
variable $C$ that causally influences $X$ and $Y$, and 4) $X$ and $Y$ are
statistically independent. We call a problem to decide which among those cases
from data _Reichenbach problem_ and deal with it. The solution for this
problem is advantageous over the traditional methods in the sense that it does
not require the prior knowledge on the nonexistence of unobserved confounder
and then it is widely applicable.

There has been many existing studies on the detection of unobserved common
cause as well. As we discuss in later sections, those methods require a type
of assumptions such that unobserved confounding variables and the observed
variables can be described by a specific formulation. However, we consider
that such assumptions on unobserved variables are hard to guarantee, even if
the domain knowledge on the dataset is available. Therefore, we aim to deal
with the Reichenbach problem without assuming such a specificity on the
possible relations between unobserved confounder and observed variables.

In our previous study [10], we proposed CLOUD (CodeLength-based methOd for
Unobserved common causes between Discrete data) to address the Reichenbach
problem for discrete data. In this paper, we extend CLOUD to accomodate all
types of data including discrete, mixed, and continuous. Unlike all existing
methods, CLOUD does not specify a form of unobserved confounders. We take a
strategy in which we select a causal model which yields the minimum codelength
among all candidates, known as the minimum description length (MDL) principle
[26]. The key for our method is to employ normalized maximum likelihood (NML)
code to compute the codelength in models of different capacities. We show this
method exhibits high performance both in theoretically and experimentally.

The rest of this paper is organized as follows. The next section introduces
the MDL principle and Structural Causal Models (SCMs), and then defines the
Reichenbach problem formally. In section 3, we review existing methods from
the viewpoint of two approaches and then see how it is hard to deal with the
Reichenbach problem without assumptions in both approaches. In section 4, we
describe models for which we consider the NML code and proposed method. In
section 5 and section 6, we conduct theoretical analysis and extensive
experiments, respectively. Finally, the conclusion is given in section 7.

##  2 Preliminaries

We introduce three theoretical frameworks upon which we construct our method
for causal discovery, namely the minimum description length (MDL) principle,
Structural Causal Models (SCMs) and the Reichenbach problem.

###  2.1 The MDL principle and the NML Codelength

The Minimum Description Length (MDL) principle is a model selection principle
grounded in the concept of data compression. It asserts that the optimal model
is the one that most succinctly describes both the data $x^{n}$ and the model
$M$, where $x^{n}=(x_{i})_{i=1,\ldots,n}$ is a data sequence of length $n$. In
this principle, we compute the codelength of data with a universal code for
each model, as an information criteria.

In this context, we introduce the Normalized Maximum Likelihood (NML) code as
a universal code. The NML code is justified by the fact that, it is optimal in
terms of the minimax regret criterion [29] and that it exhibits consistency in
model selection [28]. The NML codelength, also known as stochastic complexity,
is derived from the NML distribution. The NML distribution for statistical
model $M$ with respect to data $x^{n}$ is defined as follows:

$P_{\mathrm{NML}}\left(x^{n};M\right)=\frac{P\left(x^{n};M,\hat{\bm{\theta}}(x^{n})\right)}{\sum_{X^{n}}P\left(X^{n};M,\hat{\bm{\theta}}(X^{n})\right)},$
(1)

where $\hat{\bm{\theta}}(x^{n})$ denotes the maximum likelihood estimator of
the parameters of the model $M$ given the data $x^{n}$ and the summation
$\sum_{X^{n}}$ is taken over the space of all the possible values of the data
$X^{n}\in\mathcal{X}^{n}$.

The stochastic complexity, or the NML codelength, of $x^{n}$ is the negative
logarithmic likelihood of the NML distribution:

$\displaystyle\mathcal{SC}(x^{n};M)\coloneqq$
$\displaystyle-\log{P_{\mathrm{NML}}\left(x^{n};M\right)}$ $\displaystyle=$
$\displaystyle-\log
P\left(x^{n};M,\hat{\bm{\theta}}\left(x^{n}\right)\right)+\log\sum\limits_{X^{n}}P\left(X^{n};M,\hat{\bm{\theta}}\left(X^{n}\right)\right)$
(2)

The first term of Eq. (2) is the negative maximum log-likelihood, which is
efficiently computable for many models. The second term, referred to as
parametric complexity, is expressed as follows:

$\displaystyle\log\mathcal{C}_{n}(M)\coloneqq\log\sum\limits_{X^{n}}P\left(X^{n};M,\hat{\bm{\theta}}\left(X^{n}\right)\right)$
(3)

The parametric complexity is not always analytically tractable, and its
computation is one of the major focuses of the NML-based model selection.
Techniques for its computation include deriving asymptotically consistent
approximations [29] and employing the g-function [5].

As we see in below, the parametric complexity of a categorical distribution
model with $K$ categories, denoted by $\log\mathcal{C}_{\mathsf{CAT}}(K,n)$,
is represented by the following equation:

$\displaystyle\mathcal{C}_{\mathsf{CAT}}(K,n)=$
$\displaystyle\sum_{X^{n}\in\left\\{1,\ldots,K\right\\}^{n}}\prod_{k=1}^{K}\left(\frac{n(X=k)}{n}\right)^{n(X=k)}.$

Here, $n(X=k)$ is the frequency of occurrence of value $k$ in the sequence
$x^{n}$, defined as $n(X=k)=\sum_{i=1}^{n}\mathbb{I}(x_{i}=k)$, where
$\mathbb{I}$ is the indicator function. Kontkanen and Myllymäki developed an
efficient recurrence formula for this model with a linear time complexity of
$\mathcal{O}(n+K)$ [12]:

$\displaystyle\mathcal{C}_{\mathsf{CAT}}(K=1,n)$ $\displaystyle=1,$
$\displaystyle\mathcal{C}_{\mathsf{CAT}}(K=2,n)$
$\displaystyle=\sum_{h_{1}+h_{2}=n}\frac{n!}{h_{1}!h_{2}!}\left(\frac{h_{1}}{n}\right)^{h_{1}}\left(\frac{h_{2}}{n}\right)^{h_{2}},$
$\displaystyle\mathcal{C}_{\mathsf{CAT}}(K+2,n)$
$\displaystyle=\mathcal{C}_{\mathsf{CAT}}(K+1,n)+\frac{n}{K}\mathcal{C}_{\mathsf{CAT}}(K,n).$

In dealing with continuous variables, the summation in Eq. (3) is replaced
with an integral. Lastly, we illustrate how to compute the NML codelength for
some statistical models. These examples cover data sequences of either
discrete-type or continuous-type.

######  Example 1 (The NML Codelength for a Discrete Data).

We compute the NML codelength for a data sequence $x^{n}$ under a categorical
distribution model ${\mathsf{CAT}}^{m_{X}}$ with $m_{X}$ categories:

$\displaystyle\mathsf{CAT}^{m_{X}}=\left\\{P(X;\bm{\theta})\ \middle|\
\bm{\theta}=(\theta_{0},\dots,\theta_{m_{X}-1}),\theta_{k}\geq
0,\sum_{k=0}^{m_{X}-1}\theta_{k}=1\right\\}.$

For given data $x^{n}$, the maximum likelihood estimator of parameter
$\hat{\bm{\theta}}=(\hat{\theta}_{0},\dots,\hat{\theta}_{m_{X}-1})$ is
$\hat{\theta}_{k}(x^{n})=\frac{n(X=k)}{n}$ for $k=0,\cdots,m_{X}-1$.
Consequently, the first term in Eq. (2) is computed as

$\displaystyle-\log
P(x^{n};\mathsf{CAT}^{m_{X}},\hat{\bm{\theta}}(x^{n}))=-\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}.$

The second term in Eq. (2), the parametric complexity of
$\mathsf{CAT}^{m_{X}}$, is given by
$\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)$. Thus, the NML codelength for the
discrete data $x^{n}$ under the $\mathsf{CAT}^{m_{X}}$ is given by

$\displaystyle\mathcal{SC}(x^{n};\mathsf{CAT}^{m_{X}})=-\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n).$
(4)

######  Example 2 (The NML codelength for a Continuous Data).

Let the domain of continuous data sequence $x^{n}$ is $\mathcal{X}=[0,1)$. We
divide $\mathcal{X}$ into $m_{X}$ cells $\\{I^{X}_{k}\\}$ of equal length
$\frac{1}{m_{X}}$:

$\displaystyle\mathcal{X}=\cup_{k=0}^{m_{X}-1}I^{X}_{k},\
I^{X}_{k}=\left[\frac{k}{m_{X}},\frac{k+1}{m_{X}}\right)(k=0,\ldots,m_{X}-1)$

We define the histogram density function model $\mathsf{HIS}^{m_{X}}$ as
follows:

$\displaystyle\mathsf{HIS}^{m_{X}}\\!=\Bigg{\\{}p(X;\bm{\theta})=\\!\\!\sum_{k=0}^{m_{X}-1}\\!\theta_{k}\mathbb{I}[X\in
I^{X}_{k}]\,\Bigg{|}\,\bm{\theta}=(\theta_{0},\dots,\theta_{m_{X}-1}),\theta_{k}\geq
0,\\!\sum_{k=0}^{m_{X}-1}\\!\frac{\theta_{k}}{m_{X}}=1\Bigg{\\}}.$

The maximum likelihood estimator for this model results in
$\hat{\theta}_{k}(x^{n})=\frac{n(X\in I^{X}_{k})}{n}m_{X}$. Consequently, the
maximum log-likelihood of data is calculated as:

$\displaystyle\log p(x^{n};\hat{\bm{\theta}}(x^{n}))$
$\displaystyle=\sum_{k=0}^{m_{X}-1}n(X\in I^{X}_{k})\log\hat{\theta}_{k}$
$\displaystyle=\sum_{k=0}^{m_{X}-1}n(X\in I^{X}_{k})\log\frac{n(X\in
I^{X}_{k})}{n}+n\log m_{X},$

where $n(X\in I^{X}_{k})$ is the frequency of observations in interval
$I^{X}_{k}$ in the data sequence $x^{n}$, formally defined as $n(X\in
I^{X}_{k})=\sum_{i=1}^{n}\mathbb{I}(x_{i}\in I^{X}_{k})$. Let
$\tilde{x}^{n}=(\tilde{x}_{1},\ldots\tilde{x}_{n})\in\\{0,\ldots,m_{X}-1\\}^{n}$
be a sequence of $n$ bin labels and $I_{\tilde{x}^{n}}^{X}\coloneqq
I_{\tilde{x}_{1}}^{X}\times\cdots\times I_{\tilde{x}_{n}}^{X}\subset[0,1)^{n}$
be the $n$-dimensional hyper-bin associated with $\tilde{x}^{n}$. The
parametric complexity of $\mathsf{HIS}^{m_{X}}$ is then given by:

$\displaystyle\log\mathcal{C}_{n}(\mathsf{HIS}^{m_{X}})$
$\displaystyle=\log\int_{[0,1)^{n}}p(x^{n};\hat{\bm{\theta}}(x^{n}))dx^{n}$
$\displaystyle=\log\sum_{\tilde{x}^{n}\in\\{0,\ldots,m_{X}-1\\}^{n}}\int_{I_{\tilde{x}^{n}}^{X}}p(x^{n};\hat{\bm{\theta}}(x^{n}))dx^{n}$
$\displaystyle=\log\sum_{\tilde{x}^{n}\in\\{0,\ldots,m_{X}-1\\}^{n}}\prod_{k=0}^{m_{X}-1}\left(\frac{n(\tilde{X}=k)}{n}\right)^{n(\tilde{X}=k)}$
$\displaystyle=\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n),$

which results in the same value as the parametric complexity of the
$m_{X}$-valued categorical model $\mathsf{CAT}^{m_{X}}$.

The NML codelength thus becomes

$\displaystyle\mathcal{SC}(x^{n};\mathsf{HIS}^{m_{X}})$ $\displaystyle=$
$\displaystyle-\sum_{k=0}^{m_{X}-1}n(X\in I^{X}_{k})\log\frac{n(X\in
I^{X}_{k})}{n}-n\log m_{X}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n),$ (5)

for continuous data $x^{n}$ based on $\mathsf{HIS}^{m_{X}}$.

###  2.2 Structural Causal Model

A Structural Causal Model (SCM) [20] represents the data-generating process
through a set of structural assignments. In an SCM, variables are expressed as
functions of their parent variables (direct causes) and exogenous variables.
In particular, when we consider only two variables $X$ and $Y$, which are both
one-dimentional, with the causal graphs $X\to Y$ or $X\leftarrow Y$, SCMs can
be represented as follows:

$\displaystyle M_{X\to Y}:\begin{cases}X=E_{X}\\\ Y=f(X,E_{Y})\end{cases}$
$\displaystyle M_{X\leftarrow Y}:\begin{cases}X=g(Y,E_{X})\\\
Y=E_{Y},\end{cases}$

where $f$ and $g$ are functions, and $E_{X},E_{Y}$ are one-dimensional
exogenous variables such that $E_{X}\mathop{\\!\perp\\!\\!\\!\perp\\!}E_{Y}$.
Here, the statistical models $M_{X\to Y},M_{X\leftarrow Y}$, derived from each
SCM, are referred to as causal models. In general, without constraints on the
distributions of the exogenous variables and/or on the forms of functions $f$
and $g$, it is not identifiable whether samples from the joint distribution
$P(X,Y)$ are induced by the causal relationship of $M_{X\to Y}$ or
$M_{X\leftarrow Y}$ [25, Proposition 4.1]. In other words, causal discovery
from observational data usually requires making specific assumptions on the
functional forms and/or the distributions of exogenous variables, which limits
the scope of the joint distributions led by SCMs. Some notable SCMs in the
context of causal discovery are listed in below.

######  Example 3 (Additive Noise Model (ANM)).

ANM[6, 22, 24] assumes a data generating process as per the following
equations where effects are the nonlinear functions of their causes with
additive noise:

$\displaystyle M_{X\to Y}:\begin{cases}X=E_{X}\\\ Y=f(X)+E_{Y},\end{cases}$
$\displaystyle M_{X\leftarrow Y}:\begin{cases}X=g(Y)+E_{X}\\\
Y=E_{Y},\end{cases}$

where additive noises $E_{X},E_{Y}$ satisfy
$E_{X}\mathop{\\!\perp\\!\\!\\!\perp\\!}Y,E_{Y}\mathop{\\!\perp\\!\\!\\!\perp\\!}X$,
respectively.

######  Example 4 (Linear NonGaussian Acyclic Model (LiNGAM)).

LiNGAM[32, 33], which is a special case of ANM, assumes that causal
relationships are linear and that exogenous variables follow non-Gaussian
distributions:

$\displaystyle M_{X\to Y}:\begin{cases}X=E_{X}\\\ Y=b_{X\to
Y}X+E_{Y},\end{cases}$ $\displaystyle M_{X\leftarrow
Y}:\begin{cases}X=b_{X\leftarrow Y}Y+E_{X}\\\ Y=E_{Y},\end{cases}$

where $b_{X\to Y},b_{X\leftarrow Y}\in\mathbb{R}$ are the linear coefficients,
and $E_{X},E_{Y}$ follow Non-Gaussian distributions.

######  Example 5 (Linear Mixed causal model (LiM)).

LiM [39] is an extension of LiNGAM to accommodate mixed data types, including
both continuous and discrete variables. In the case of continuous variables,
LiM assumes the same SCMs as LiNGAM. When $X$ is a continuous variable and $Y$
is a binary variable, LiM assumes the following SCM:

$\displaystyle M_{X\to Y}:\begin{cases}X=E_{X}\\\ Y=\begin{cases}1\ (b_{X\to
Y}X+E_{Y}>0),\\\ 0\ (\text{otherwise}),\\\ \end{cases}\\\ \end{cases}$

where $E_{Y}$ follows a Logistic distribution.

######  Example 6 (LiNGAM with latent confounder (lvLiNGAM)).

lvLiNGAM extends the basic LiNGAM model to incorporate hidden common causes
[7]. Besides the linear and non-Gaussian assumptions of LiNGAM, lvLiNGAM
explicitly models unobserved common causes $C$. For a one-dimensional $C$, the
model can be represented with the following SCM:

$\displaystyle M_{X\leftarrow C\to Y}:\begin{cases}X=\lambda_{X}C+E_{X}\\\
Y=\lambda_{Y}C+E_{Y},\\\ \end{cases}$

where $\lambda_{X},\lambda_{Y}\in\mathbb{R}$ denote the direct causal effects
from the unobserved common cause $C$ to each observed variable to the observed
variables $X$ and $Y$, respectively. In this SCM, the latent variable $C$ is
assumed to be non-Gaussian and independent of $E_{X}$ and $E_{Y}$, and it is
assumed to have linear effects on the observed variables.

###  2.3 Reichenbach Problem

This section describes the Reichenbach problem formally, which is central to
our study. Suppose we have i.i.d. observational data
$z^{n}=(x^{n},y^{n})\in\mathcal{X}^{n}\times\mathcal{Y}^{n}$ generated from
joint distribution $P(X,Y)$. Here, $X$ and $Y$ can be either discrete or
continuous variables.

Based on Reichenbach’s common cause principle, we can categorize the causal
relationship between $X$ and $Y$ into four cases. The goal is to infer one of
these causal models $M$ that best explains the underlying causal relationship:

  * •

$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$: $X$ and $Y$ are independent, with
no direct causal link.

  * •

$M_{X\leftarrow C\to Y}$: There exist common causes $C$ that causally
influence both $X$ and $Y$.

  * •

$M_{X\to Y}$: $X$ causes $Y$, but not vice versa.

  * •

$M_{X\leftarrow Y}$: $Y$ causes $X$, but not vice versa.

In solving the Reichenbach problem, it is desirable to use methods that do not
make any assumption about unobserved variables $C$. This is crucial because it
is almost impossible to have prior knowledge of all potential unobserved
variables. Therefore, we propose a method capable of detecting the presence of
unobserved common causes, even in situations where the candidates for latent
variables are unknown or when unexpected unobserved variables are present,
without relying on assumptions about $C$.

To employ existing causal discovery methods that rely on modeling unobserved
variables, one needs to know the nature of the unobserved causes beforehand.
However, it is challenging to acquire complete knowledge about all possible
unobserved variables solely from the domain knowledge of the observed
variables. Thus, when applying these methods to real-world data, a heightened
level of caution is required. Nevertheless, existing approaches to the
Reichenbach problem typically depend on models that make assumptions about the
relationships between unobserved variables $C$ and the observed variables
$X,Y$.

##  3 Existing Work

Existing causal discovery methods from a joint distribution of two variables
can be categorized into two approaches: one employs the identifiability of the
model and the other is based on the principle of algorithmic independence of
conditionals. Furthermore, in the context of the Reichenbach problem, there
are methods focused on solving the sub-problem of choosing between $X\to Y$
and $X\leftarrow Y$, and those that attempt to solve the Reichenbach problem
by making assumptions about unobserved common causes. In this section, we
describe these approaches and discuss why it is difficult to detect unobserved
common factors without making assumptions about these causes.

###  3.1 Identifiable models

In this approach, when we formulate causal relationships using SCMs, we
restrict the functional forms and the distributions of the exogenous
variables. The causal models induced by SCMs are said to be identifiable if
different causal structures always lead to different joint distributions of
the observed variables.

In order to infer a causal structure using identifiable models, we assume that
the observed data have been generated by a distribution belonging to one of
these models. Then, we infer $X$ is the cause of $Y$ when the corresponding
model explains the data best of all models. This reasoning applies similarly
to other causal relationships.

To determine the causal direction, various type of indentifiable models are
studied in existing work so that various types of data can be applied. Shimizu
et al. [32] showed that causal models become identifiable if function is
linear and the distribution of exogenous variables is non-Gaussian as in
Example 4 (LiNGAM). A general case of LiNGAM is Additive Noise Models (ANMs,
[6]), as detailed in Example 3, where we assume that the noise is additive and
independent of the cause. In general, ANM is identifiable if functional form
is non-linear even without imposing any restrictions on the noise
distributions [6, 24]. [22] proposed DR, which is an extension of ANMs to
discrete variables and showed that ANM is generally identifiable in discrete
case. For mixed-type data, [39] formulated Linear Mixed causal model (LiM) as
shown in Example 5, and Li et al. [14] proposed an algorithm, HCM, which
formulates nonlinear causal relationships as mixed-SCMs.

In this approach, we can recover the true causal relationship based on
observed data as long as there is no latent confounder, i.e., the true
distribution belongs to one of the models we adopt. Therefore, this framework
requires practitioners to make use of domain knowledge such that causal
relationship, if any, can be formulated in a certain type of SCM. However,
causal discovery methods under the assumption of causal sufficiency, which
assume no unobserved common causes, can lead to incorrect conclusions when
applied to data where unobserved common causes actually exist. To avoid such
issues, it becomes crucial to consider unobserved common causes in causal
discovery.

As for models that allow for the presence of unobserved common causes, Hoyer
et al. proposed a model called lvLiNGAM [7], which extends LiNGAM to the
models with latent counfounding variables by explicitly modeling them.
lvLiNGAM, shown in Example 6, assumes that the latent confounder $C$ follow
non-Ganssian distributions and relationships between $C$ and observed
variables are linear. Parce LiNGAM (BUPL, [37]) and Repetitive Causal
Discovery (RCD, [16]) make the same assumptions on SCMs. Maeda et al. extended
lvLiNGAM to its non-linear variant, and proposed CAMUV algorithm [17] .

For models to be identifiable including $M_{X\leftarrow C\to Y}$, one must
make an assumption on the model of unobserved common causes. This means that
one must have known the nature of the unobserved common causes beforehand to
successfully perform causal discovery. Otherwise, conclusions led by this
framework will be unreliable when unobserved common causes do not follow the
assumptions. Even with domain knowledge of $X$ and $Y$, it remains challenging
to accurately determine the form of SCM for an unobserved common cause, which
a practitioner is not certain if exists, will satisfy.

###  3.2 Algorithmic Independence of Conditionals

This section describes the approach based on the principle of algorithmic
independence of conditionals, which is described as follows: if true causality
is $X\to Y$, then mechanism $P^{*}(Y|X)$ is independent of the cause
$P^{*}(X)$ [8], where $P^{*}$ denotes true distributions we assume under the
corresponding causal relationship. By denoting the Kolmogorov complexity as
$K$, it leads to the following inequality:

$K(P^{*}(X))+K(P^{*}(Y|X))<K(P^{*}(Y))+K(P^{*}(X|Y)),$

if true causality is $X\to Y$ [35]. This inequality can not be evaluated due
to the following two reasons: Kolmogorov complexity is not computable and true
distribution is unknown. Therefore, Marx and Vreeken [18] has justified that
this principle leads to the approximation build upon two-part MDL as follows:

$\displaystyle\mathcal{L}(z^{n};M_{X\to Y})<\mathcal{L}(z^{n};M_{X\leftarrow
Y}),$

where $\mathcal{L}(z^{n};M_{X\to Y})$ is the description length of data
$z^{n}$ under statistical model $M_{X\to Y}=M_{X}\times M_{Y|X}$ in which we
assume $P^{*}(X)\in M_{X}$ and $P^{*}(Y|X)\in M_{Y|X}$, formally defined as
follows:

$\displaystyle\mathcal{L}(z^{n};M_{X\to
Y})=L(x^{n};M_{X})+L(y^{n};x^{n},M_{Y|X}).$

We define $\mathcal{L}(z^{n};M_{X\leftarrow Y})$ analogously. Methods such as
Causal Inference by Stochastic Complexity (CISC, [2]), Accurate Causal
Inference on Discrete data (ACID, [3]) and Distance Correlation (DC, [15])
model $\mathcal{L}(z^{n};M)$ using refined MDL, Shannon Entoropy and distance
correlation, respectively.

If one considers including confounded model $M_{X\leftarrow C\to Y}$ with
unobserved common causes $C$, the description length under the joint
distributions of that model, $\mathcal{L}(z^{n};M_{X\leftarrow C\to Y})$ which
is an approximation of $K(P^{*}(X,Y,C))$, must be evaluated and compared with
directed cases of $M_{X\to Y}$ and $M_{X\leftarrow Y}$. A naive approximation
approach based solely on likelihood invariably leads to the selection of the
confounded model $M_{X\leftarrow C\to Y}$, due to its inherently minimized
complexity. To address this issue, Confounded-or-Causal (COCA) method was
developed [9]. COCA selects between $X\to Y$ and $X\leftarrow C\to Y$ under
the assumption that not only the observed variables but also unobserved common
causes $C$ follow specific-dimensional Gaussian distributions. It employs
Bayesian coding to approximate the description length of data, considering not
only the likelihood but also the complexity of the model class, including $C$.
This approach enables to comparison between different sizes of statistical
models, specifically $M_{X\to Y}$ and $M_{X\leftarrow C\to Y}$. However, it is
important to note that this approach relies on certain assumptions about $C$.

In Summary, in challenging the Reichenbach problem, all existing methods face
a common limitation: they require additional assumptions about unobserved
common causes. Our previous work has already shown that CLOUD successfully
overcomes this limitation in the context of discrete variables [10]. We
claimed that this is achievable by comparing models with different capacities,
as quantified using the NML codelength. In this paper, our objective is to
expand the applicability of CLOUD to encompass continuous and mixed data
types, thereby enhancing its effectiveness in solving the Reichenbach problem
across a wider range of data types.

##  4 Proposed Method

In this section, we extend CLOUD to cover all data types, which was originally
designed for the Reichenbach problem in discrete data. In CLOUD, we formulate
causal models for four causal relationships in the Reichenbach problem
(Section 4.1). Then, based on MDL principle, we calculate the codelength of
the observed data $z^{n}$ based on NML coding and select a causal model $M$
which achieves the shortest codelength (Section 4.2).

We formulate confounded model $M_{X\leftarrow C\to Y}$ to represent any joint
distribution, thus avoiding assumptions about $C$. Since this model has the
highest complexity compared to the other three, the NML-based codelength on
this model is not necessarily the shortest although the negative loglikilhood
is the smallest. Thus, by considering both model complexity and data
likelihood, we can select an appropriate causal model among models of
different complexities.

###  4.1 Model

In this section, we describe causal models for cases where both $X$ and $Y$
are discrete or continuous variables, represented as statistical models
derived from the assumed Structural Causal Models (SCMs) for each causal
relationship.

First, we formulate SCMs for each causal model $M$ to describe the causal
relationships between $X\in\mathcal{X}$ and $Y\in\mathcal{Y}$:

$\displaystyle
M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}:\begin{cases}X=E_{X}\\\
Y=E_{Y}\end{cases}$ $\displaystyle M_{X\leftarrow C\to
Y}:\begin{cases}X=f(C,E_{X})\\\ Y=g(C,E_{Y})\end{cases}$ $\displaystyle
M_{X\to Y}:\begin{cases}X=E_{X}\\\ Y=f(X)+E_{Y}\end{cases}$ $\displaystyle
M_{X\leftarrow Y}:\begin{cases}X=g(Y)+E_{X}\\\ Y=E_{Y}\end{cases}$

Here, the exogenous variables $E_{X}\in\mathcal{X},E_{Y}\in\mathcal{Y}$ are
independent of each other. For $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$, we
assume faithfullness [34] in the sense that we regard $X$ and $Y$ is causally
independent when they are statistically independent.

For $M_{X\leftarrow C\to Y}$, any probability density functions on
$(\mathbb{Z}/m_{X}\mathbb{Z})\times(\mathbb{Z}/m_{Y}\mathbb{Z})$ and
probability density functions on $(\mathbb{R}/\mathbb{Z})^{2}$ belong to the
$M_{X\leftarrow C\to Y}$ by considering appropriate choice of $f$, $g$ and
$C$. In this sense, we do not impose assumptions on the unobserved common
cause $C$ for confounded case.

For $M_{X\to Y}$ and $M_{X\leftarrow Y}$, we assume additive noise models
(ANMs) [22], which employs functions from the function sets
$\mathcal{F}=\\{f:\mathcal{X}\to\mathcal{Y}\mid f\text{ is not constant}\\}$
and $\mathcal{G}=\\{g:\mathcal{Y}\to\mathcal{X}\mid g\text{ is not
constant}\\}$. In case $X$ is discrete, we set
$\mathcal{X}=\\{0,1,\dots,m_{X}-1\\}$ and addition is taken over
$\mathbb{Z}/m_{X}\mathbb{Z}$ as in [25, 23]. In case $X$ is continuous, we set
$\mathcal{X}=[0,1)$ and addition is taken over $\mathbb{R}/\mathbb{Z}$. The
same goes for cases $Y$ is either discrete or continuous. In continuous cases,
the addition can be regarded as addition over $\mathbb{R}$ for data scaled
with a sufficiently small constant $\epsilon$. This implies that, in practical
applications, models with periodic boundary conditions can be considered as
including non-periodic ANMs.

Second, we identify the causal models
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$, $M_{X\leftarrow C\to Y}$, $M_{X\to
Y}$ and $M_{X\leftarrow Y}$ with a set of joint probability distributions on
$(X,Y)$ that models imply. It can be justified in case in which only
observations from the joint distribution are available.

##### Discrete Case:

If both $X$ and $Y$ are discrete, then
$\mathcal{X}=\mathbb{Z}/m_{X}\mathbb{Z}$ and
$\mathcal{Y}=\mathbb{Z}/m_{Y}\mathbb{Z}$. Any discrete probability
distribution on $(X,Y)$ can be identified with a parameter
$\bm{\theta}\in\bm{\Theta}$ by the following relation:

$\displaystyle
P(X,Y;\bm{\theta})=\prod_{k,k^{\prime}}\theta_{k,k^{\prime}}^{\mathbb{I}[X=k,Y=k^{\prime}]},$

where we define
$\bm{\Theta}=\left\\{\bm{\theta}=(\theta_{k,k^{\prime}})_{k,k^{\prime}}\in\mathbb{R}^{m_{X}\times
m_{Y}}\middle|\theta_{k,k^{\prime}}\geq
0,\sum\theta_{k,k^{\prime}}=1\right\\}$. Based on this parametrization, we
represent causal models by characterizing the respective subset of the
parameter space as follows:

$\displaystyle M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$
$\displaystyle=\left\\{P\left(X,Y;\bm{\theta}\right)=\prod_{k,k^{\prime}}\theta_{k,k^{\prime}}^{\mathbb{I}[X=k,Y=k^{\prime}]}\,\middle|\,\bm{\theta}\in\bm{\Theta}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\right\\},$
$\displaystyle M_{X\leftarrow C\to Y}$
$\displaystyle=\left\\{P\left(X,Y;\bm{\theta}\right)=\prod_{k,k^{\prime}}\theta_{k,k^{\prime}}^{\mathbb{I}[X=k,Y=k^{\prime}]}\,\middle|\,\bm{\theta}\in\bm{\Theta}_{X\leftarrow
C\to Y}\right\\},$ $\displaystyle M_{X\to Y}$
$\displaystyle=\left\\{P\left(X,Y;\bm{\theta}\right)=\prod_{k,k^{\prime}}\theta_{k,k^{\prime}}^{\mathbb{I}[X=k,Y=k^{\prime}]}\,\middle|\,\bm{\theta}\in\bm{\Theta}_{X\to
Y}\right\\},$ $\displaystyle M_{X\leftarrow Y}$
$\displaystyle=\left\\{P\left(X,Y;\bm{\theta}\right)=\prod_{k,k^{\prime}}\theta_{k,k^{\prime}}^{\mathbb{I}[X=k,Y=k^{\prime}]}\;\middle|\,\bm{\theta}\in\bm{\Theta}_{X\leftarrow
Y}\right\\},$

where

$\displaystyle\bm{\Theta}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$
$\displaystyle=\left\\{(\theta^{X}_{k}\,\theta^{Y}_{k^{\prime}})_{k,k^{\prime}}\in\bm{\Theta}\,\middle|\,\sum_{k}\theta^{X}_{k}=1,\theta^{X}_{k}\geq
0,\sum_{k^{\prime}}\theta^{Y}_{k^{\prime}}=1,\theta^{Y}_{k^{\prime}}\geq
0\right\\},$ (6a) $\displaystyle\bm{\Theta}_{X\leftarrow C\to Y}$
$\displaystyle=\bm{\Theta},$ (6b) $\displaystyle\bm{\Theta}_{X\to Y}$
$\displaystyle=\left\\{(\theta^{X}_{k}\,\theta^{Y}_{f(k)+k^{\prime}})_{k,k^{\prime}}\in\bm{\Theta}\,\middle|\,\sum_{k}\theta^{X}_{k}=1,\theta^{X}_{k}\geq
0,\sum_{k^{\prime}}\theta^{Y}_{k^{\prime}}=1,\theta^{Y}_{k^{\prime}}\geq
0,f\in\mathcal{F}\right\\},$ (6c) $\displaystyle\bm{\Theta}_{X\leftarrow Y}$
$\displaystyle=\left\\{(\theta^{X}_{g(k^{\prime})+k}\,\theta^{Y}_{k^{\prime}})_{k,k^{\prime}}\in\bm{\Theta}\,\middle|\,\sum_{k}\theta^{X}_{k}=1,\theta^{X}_{k}\geq
0,\sum_{k^{\prime}}\theta^{Y}_{k^{\prime}}=1,\theta^{Y}_{k^{\prime}}\geq
0,g\in\mathcal{G}\right\\}.$ (6d)

Note that the addition in subscripts is taken over each respective finite
space. For discrete case, we encode the data $z^{n}$ based on these discrete
causal models.

##### Continuous Case:

In case the domains of $X$ and $Y$ are both continuous, we represent their
joint density function by the infinite union of histogram densities.

Firstly, we consider partitioning $\mathcal{X}$ into $m_{X}$ equally-sized
cells $\\{I^{X}_{k}\\}\ (k=0,\ldots,m_{X}-1)$ and $\mathcal{Y}$ into $m_{Y}$
equally-sized cells $\\{I^{Y}_{k^{\prime}}\\}\ (k^{\prime}=0,\ldots,m_{Y}-1)$.

We then define the two-dimensional density function model
$\mathsf{HIS}^{m_{X},m_{Y}}$ as follows:

$\displaystyle\mathsf{HIS}^{m_{X},m_{Y}}\\!=\left\\{p(X,Y;\bm{\theta})=\\!\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}\\!\theta_{k,k^{\prime}}\mathbb{I}[X\in
I^{X}_{k},Y\in
I^{Y}_{k^{\prime}}]\,\middle|\,\bm{\theta}\in\bm{\Theta}^{m_{X},m_{Y}}\right\\}$

where

$\displaystyle\bm{\Theta}^{m_{X},m_{Y}}=\left\\{\bm{\theta}=(\theta_{k,k^{\prime}})_{k,k^{\prime}}\in\mathbb{R}^{m_{X}\times
m_{Y}}\,\middle|\,\theta_{k,k^{\prime}}\geq
0,\\!\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}\\!\frac{\theta_{k,k^{\prime}}}{m_{X}m_{Y}}=1\right\\}.$

Then, the infinite union $\bigcup_{m_{X},m_{Y}}\mathsf{HIS}^{m_{X},m_{Y}}$ can
represent any joint density function n an arbitrary precision level. Due to
the universality of $M_{X\leftarrow C\to Y}$, we represent

$\displaystyle M_{X\leftarrow C\to
Y}=\bigcup_{m_{X},m_{Y}}M^{m_{X},m_{Y}}_{X\leftarrow C\to Y},\quad
M^{m_{X},m_{Y}}_{X\leftarrow C\to Y}=\mathsf{HIS}^{m_{X},m_{Y}}.$

It means that we can identify any joint distribution on $(X,Y)$ by an element
of $M_{X\leftarrow C\to Y}$, in the sense that it gives the same distribution
of codelength for a given code when the approximation level is fixed.

By defining
$\bm{\Theta}^{m_{X},m_{Y}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},\bm{\Theta}^{m_{X},m_{Y}}_{X\to
Y},\bm{\Theta}^{m_{X},m_{Y}}_{X\leftarrow Y}\subset\bm{\Theta}^{m_{X},m_{Y}}$
similarly as equations (6), we can represent the other models as follows:

$\displaystyle
M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}=\bigcup_{m_{X},m_{Y}}M^{m_{X},m_{Y}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},$
$\displaystyle\quad
M^{m_{X},m_{Y}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}=\left\\{p(X,Y;\bm{\theta})\,\middle|\,\bm{\theta}\in\bm{\Theta}^{m_{X},m_{Y}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\right\\}$
$\displaystyle M_{X\to Y}=\bigcup_{m_{X},m_{Y}}M^{m_{X},m_{Y}}_{X\to Y},$
$\displaystyle\quad M^{m_{X},m_{Y}}_{X\to
Y}=\left\\{p(X,Y;\bm{\theta})\,\middle|\,\bm{\theta}\in\bm{\Theta}^{m_{X},m_{Y}}_{X\to
Y}\right\\}$ $\displaystyle M_{X\leftarrow
Y}=\bigcup_{m_{X},m_{Y}}M^{m_{X},m_{Y}}_{X\leftarrow Y},$ $\displaystyle\quad
M^{m_{X},m_{Y}}_{X\leftarrow
Y}=\left\\{p(X,Y;\bm{\theta})\,\middle|\,\bm{\theta}\in\bm{\Theta}^{m_{X},m_{Y}}_{X\leftarrow
Y}\right\\}.$

We can regard $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ as a set of any
density function that is decomposable as $p(X,Y)=p(X)p(Y)$. As for $M_{X\to
Y}$, we can see it as a special form of decomposition of the density function
such as $p(X,Y)=p(X)p(Y-f(X))$ using a function $f$.

###  4.2 Algorithm

We regard the Reichenbach problem as a problem of model selection and conduct
a selection under the MDL criterion. That is, among the causal models defined
in Section 4.1, we infer the underlying causal relationship of the data
$z^{n}$ is such a causal model $M$ that yields the shortest description length
of the data.

##### Discrete Case:

Algorithm 1 Main function in Discrete Case

1:Data $z^{n}$, A set of model candidates $\mathcal{M}$

2:Best model $\hat{M}$

3:for $M$ in $\mathcal{M}$ do

4: Compute $\mathcal{L}^{d}(z^{n};M)$ using (9)

5:end for

6:$\hat{M}=\mathrm{argmin}_{M\in\mathcal{M}}\mathcal{L}^{d}(z^{n};M)$

7:return $\hat{M}$

We infer a causal relationship by selecting a discrete causal model according
to the following equation:

$\displaystyle\hat{M}(z^{n})=\mathop{\mathrm{argmin}}_{M\in\mathcal{M}}\mathcal{L}^{d}(z^{n};M),$
(7)

where $\mathcal{L}^{d}(z^{n};M)$ is a universal codelength of the discrete
data $z^{n}$ for the discrete causal model $M$, and we employ NML code to
compute the codelength. $\mathcal{M}$ is a set of model candidates,
$\mathcal{M}=\left\\{M_{X\to Y},M_{X\leftarrow
Y},M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},M_{X\leftarrow C\to Y}\right\\}$.
We can also set $\mathcal{M}$ to be its subset, such as $\left\\{M_{X\to
Y},M_{X\leftarrow Y}\right\\}$, based on the prior knowledge. The algorithm
for discrete case is shown in Algorithm 1.

We calculate $\mathcal{L}^{d}(z^{n};M)$ for
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ and $M_{X\leftarrow C\to Y}$ using
exact NML codes. For $M_{X\to Y}$ and $M_{X\leftarrow Y}$, we employ two-stage
coding based on the NML code since it is hard to exactly calculate the
codelength of the NML code because functions $f$ and $g$ are not fixed. For
$M_{X\to Y}$, we compute the codelength based on two-stage coding with respect
to function $f$ as follows:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\to Y})=L(f;M_{X\to
Y})+L(z^{n};M_{X\to Y},f),$ (8)

The first term on the right-hand side is a codelength required to encode a
function $f$, and the second term represents the NML codelength for the model
$M_{X\to Y}$ with function $f$ fixed. The same applies to $M_{X\leftarrow Y}$.
The forms of each codelength are provided in Proposition 1. We provide its
proof in Appendix B.

######  Proposition 1 (NML-based codelength for discrete data).

For a given discrete data $z^{n}$ and the discrete causal models
$M\in\left\\{M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},M_{X\to
Y},M_{X\leftarrow Y},M_{X\leftarrow C\to Y}\right\\}$, the codelengths defined
as above have the following expressions:

$\displaystyle\mathcal{L}^{d}(z^{n};M)$
$\displaystyle=\begin{dcases}\ell_{X}+\ell_{Y}+\log\left(\mathcal{C}_{\mathsf{CAT}}(m_{X},n)\cdot\mathcal{C}_{\mathsf{CAT}}(m_{Y},n)\right)&\textnormal{if
}M=M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\text{,}\\\
\ell_{X,Y}+\log\mathcal{C}_{\mathsf{CAT}}(m_{X}m_{Y},n)&\textnormal{if
}M=M_{X\leftarrow C\to Y}\text{,}\\\
\ell_{X}+\ell_{Y|X}(\hat{f})+\log\left(\mathcal{C}_{\mathsf{CAT}}(m_{X},n)\cdot\mathcal{C}_{\mathsf{CAT}}(m_{Y},n)\right)\\\
\quad+\log(m_{Y}^{m_{X}-1}-1)&\textnormal{if }M=M_{X\to Y}\text{,}\\\
\ell_{Y}+\ell_{X|Y}(\hat{g})+\log\left(\mathcal{C}_{\mathsf{CAT}}(m_{X},n)\cdot\mathcal{C}_{\mathsf{CAT}}(m_{Y},n)\right)\\\
\quad+\log(m_{X}^{m_{Y}-1}-1)&\textnormal{if }M=M_{X\leftarrow
Y},\end{dcases}$ (9)

where

$\displaystyle\ell_{X}$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n},$
$\displaystyle\ell_{Y|X}(f)$
$\displaystyle=-\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y=f(X)+k^{\prime})\log\frac{n(Y=f(X)+k^{\prime})}{n},$
$\displaystyle\ell_{X,Y}$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n\left(X=k,Y=k^{\prime}\right)\log\frac{n\left(X=k,Y=k^{\prime}\right)}{n},$

and similarity for $\ell_{Y}$ and $\ell_{X|Y}$. Here, $\hat{f}$ and $\hat{g}$
are functions derived through maximum likelihood estimation.

Proposition 1 shows that while $M_{X\leftarrow C\to Y}$ is the most expressive
model and the negative log-likelihood of the data is always minimized in
$M_{X\leftarrow C\to Y}$, its parametric complexity is the largest among all
models. Therefore, by employing NML-based codelength as shown in Eq. (9), we
can compare between models with varies capacities as per Eq. (7) under the
trade-off between data likelihood and model complexity.

The algorithm for estimation of function $\hat{f}$ or $\hat{g}$ is shown in
Algorithm 4 in Appendix A.

##### Continuous Case:

Algorithm 2 Main function in Continuous Case

1:Data $z^{n}$, a set of model candidates $\mathcal{M}$, and candidates
$\mathcal{P}$ for bin numbers $m_{X}$ and $m_{Y}$

2:Best model $\hat{M}$

3:for $(m_{X},m_{Y})$ in $\mathcal{P}$ do

4: for $M$ in $\mathcal{M}$ do

5: Compute $\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$ using (13)

6: end for

7:end for

8:return
$\hat{M}=\mathop{\mathrm{argmin}}_{M\in\mathcal{M}}\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$

We infer a causal relationship by selecting a continuous causal model $M$
according to the following equation:

$\displaystyle\hat{M}(z^{n})$
$\displaystyle=\mathop{\mathrm{argmin}}_{M\in\mathcal{M}}\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M).$
(10)

Here, $\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$ is a codelength of data $z^{n}$
required to encode $z^{n}$ up to an arbitrarily precision $\delta>0$. In order
to construct a universal code of $z^{n}$ under causal model $M$, we employ
two-stage coding for $m_{X},m_{Y}$ as follows:

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$
$\displaystyle=L(m_{X},m_{Y};M)+L(z^{n};M,m_{X},m_{Y}),$ (11)

where $L(m_{X},m_{Y};M)$ is the codelength required to encode the numbers of
bins $m_{X}$ and $m_{Y}$. The second term, $L(z^{n};M,m_{X},m_{Y})$,
represents codelengths for encoding $z^{n}$ based on continuous models of
corresponding bin sizes,
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}^{m_{X},m_{Y}},M_{X\leftarrow C\to
Y}^{m_{X},m_{Y}},M_{X\to Y}^{m_{X},m_{Y}}$, or $M_{X\leftarrow
Y}^{m_{X},m_{Y}}$.

In order to encode $z^{n}=(x^{n},y^{n})$ with the precision $\delta$, we again
employ two-part coding through $\textnormal{disc}(x^{n};m_{X})$ and
$\textnormal{disc}(y^{n};m_{Y})$. Here, we define
$\textnormal{disc}:\mathcal{X}^{n}\to\\{0,\ldots,m_{X}-1\\}^{n}$ as a function
that discretizes continuous data $x^{n}\in\mathcal{X}$ into $m_{X}$ equal
categories. Note that the discretized data is encoded by the strategy
mentioned above as in discrete case. After the discretized data is encoded,
additional codelength are required to attain the prespecified precision level.
This additional codelength is denoted by
$L(x^{n};\textnormal{disc}(x^{n};m_{X}))$ or
$L(y^{n};\textnormal{disc}(y^{n};m_{Y}))$. Thus, the second term in Eq. (11)
has the following expression.

$\displaystyle L(z^{n};M,m_{X},m_{Y})=$
$\displaystyle\mathcal{L}^{d}(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y});\textnormal{DISC}(M;m_{X},m_{Y}))$
$\displaystyle+L(x^{n};\textnormal{disc}(x^{n};m_{X}))+L(y^{n};\textnormal{disc}(y^{n};m_{Y})),$
(12)

where $\textnormal{DISC}(M;m_{X},m_{Y})$ denotes the discrete causal model of
causality $M\in\mathcal{M}$ with the category numbers $m_{X}$ and $m_{Y}$.
Consequently, the form of each codelength have an expression as provided in
Proposition 2. We provide its proof in Appendix C.

######  Proposition 2 (Codelength in Continuous Case).

The codelength of data $z^{n}$ required to encode $z^{n}$ as described above
has the following expression:

$\displaystyle\mathcal{L}^{c}(z^{n};M)$
$\displaystyle=\mathcal{L}^{d}(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y});\textnormal{DISC}(M;m_{X},m_{Y}))$
(13) $\displaystyle\quad+L^{c\to d}(m_{X},n)+L^{c\to
d}(m_{Y},n)+\textnormal{const.},$

where the first term on the right-hand side in Eq. (13) can be calculated
using Eq. (9) in Proposition 1 and $L^{c\to d}(m,n)$ for $m\in\mathbb{N}^{+}$
is defined as:

$\displaystyle L^{c\to d}(m,n)=-n\log m+\log^{*}m,$ (14)

where $\log^{*}m$ is given by Rissanen’s universal integer coding [27]

$\displaystyle\log^{*}m=\log c+\log m+\log\log m+\cdots\quad(c\approx 2.865),$
(15)

where the summation is only taken over nonnegative terms. The constant term
only depends on the precision level $\delta$.

The algorithm for the continuous case is presented in Algorithm 2. For
computational efficiency, a practical algorithm can limit search range for
$m_{X}$ and $m_{Y}$ to a subset $\mathcal{P}\subset(\mathbb{N}^{+})^{2}$.

##### Mixed Case:

We can consider the both of mixed cases as a special case of the continuous
case. If $X$ is a continuous variable and $Y$ is a discrete variable, we
regard $y^{n}=(y_{i})_{i}$ as continuous values by mapping $y^{n}$ to
$\text{cont}(y^{n};m_{Y})=(y_{i}/m_{Y})_{i}$ and calculate the description
length using continuous causal models by
$\mathcal{L}^{c}(x^{n},\text{cont}(y^{n}),m_{X},m_{Y};M)$. The causal
discovery algorithm for the mixed case is presented in Algorithm 3.

Algorithm 3 Main function in Mixed Case in which $X$ is continuous

1:Data $z^{n}=(x^{n},y^{n})$, a set of model candidates $\mathcal{M}$, search
range $\mathcal{P}\subset\mathbb{N}^{+}$ for $m_{X}$, and the number of
categories $m_{Y}$

2:Best model $\hat{M}$

3:for $m_{X}$ in $\mathcal{P}$ do

4: for $M$ in $\mathcal{M}$ do

5: Compute $\mathcal{L}^{c}(x^{n},\text{cont}(y^{n};m_{Y}),m_{X},m_{Y};M)$
using (13)

6: end for

7:end for

8:return
$\hat{M}=\mathop{\mathrm{argmin}}_{M\in\mathcal{M}}\min_{m_{X}\in\mathcal{P}}\mathcal{L}^{c}(x^{n},\text{cont}(y^{n};m_{Y}),m_{X},m_{Y};M)$

##  5 Theoretical analysis of the statistical consistency

In this section, we provide the theoretical analysis on the consistency of our
method. By consistency, we mean that the probability that our method select
the true model converges to 1 at the limit of large $n$. Noting the inclusion
relation in our case, $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},M_{X\to
Y},M_{X\leftarrow Y}\subset M_{X\leftarrow C\to Y}$, we consider the true
model for a given probability distribution is the minimal model that contains
it.

###  5.1 Discrete Case

######  Theorem 1.

Define true model $M^{*}(P^{*})$ given as follows:

$\displaystyle
M^{*}(P^{*})=\begin{cases}M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}&P^{*}\in
M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\\\ M_{X\to Y}&P^{*}\in M_{X\to Y}\\\
M_{X\leftarrow Y}&P^{*}\in M_{X\leftarrow Y}\\\ M_{X\leftarrow C\to
Y}&P^{*}\in M_{X\leftarrow C\to
Y}\setminus\left(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\cup M_{X\leftarrow
Y}\cup M_{X\leftarrow Y}\right).\end{cases}$

Then, the probability that CLOUD outputs $M^{*}(P^{*})$ given $n$ i.i.d.
samples from $P^{*}$, converges to 1 in the limit of $n\to\infty$, provided
that the maximum likelihood estimation of $\hat{f}$ and $\hat{g}$ is
successful.

###### Proof.

##### In case of $P^{*}\in M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$:

The asymptotic expansion of log-likelihood [1] implies

$\displaystyle-\log
P\left(z^{n};M,\hat{\bm{\theta}}\left(z^{n}\right)\right)=nH(P^{*})+n\textnormal{KL}(P^{*}||P)+O_{P}\left(1\right),$
(16)

for all $M\in\mathcal{M}_{\rm all}$. Here, $H(P)$ denotes the entropy of $P$
defined as

$H(P)=-\sum_{k,k^{\prime}}P(X=k,Y=k^{\prime})\log P(X=k,Y=k^{\prime}),$

$\textnormal{KL}(P^{*}||P)$ denotes the KL-divergence defined as

$K(P^{*}||P)=\sum_{k,k^{\prime}}P^{*}(X=k,Y=k^{\prime})\frac{\log
P^{*}(X=k,Y=k^{\prime})}{\log P(X=k,Y=k^{\prime})},$

and $O_{P}(\cdot)$ denotes the asymptotic order with respect to $n$ in
probability. Specifically, in case of
$M\in\left\\{M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},M_{X\leftarrow C\to
Y}\right\\}$, the asymptotic expansion of Eq. (16) becomes

$\displaystyle-\log
P\left(z^{n};M,\hat{\bm{\theta}}\left(z^{n}\right)\right)=nH(P^{*})+O_{P}\left(1\right),$
(17)

since $P^{*}$ belongs to both $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ and
$M_{X\leftarrow C\to Y}$. We see
$\log\mathcal{C}_{n}(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})<\log\mathcal{C}_{n}(M_{X\leftarrow
C\to Y})$, which leads to

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\leftarrow C\to
Y})-\mathcal{L}^{d}(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=\log\mathcal{C}_{n}(M_{X\leftarrow C\to
Y})-\log\mathcal{C}_{n}(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})+O_{P}\left(1\right)$
$\displaystyle=\Omega(\log n)+O_{P}(1).$

Thus, the probability that $M_{X\leftarrow C\to Y}$ achieves the smallest
codelength converges to 0. As for $M\in\\{M_{X\to Y},M_{X\leftarrow Y}\\}$,
the negative log-likelihood function in the first term in Eq. (2) divided by
$n$ converges to

$\displaystyle-\frac{1}{n}\log P(z^{n};\bm{\theta})$
$\displaystyle=-\sum_{k,k^{\prime}}\frac{n(X=k,Y=k^{\prime})}{n}\log
P(X=k,Y=k^{\prime};\,{\bm{\theta}})$
$\displaystyle=-\sum_{k,k^{\prime}}\theta_{k,k^{\prime}}^{*}\log
P(X=k,Y=k^{\prime};\,{\bm{\theta}})+o_{P}\left(1\right)$ (18)

as $n\to\infty$ since
$\frac{n(X=k,Y=k^{\prime})}{n}\to\theta_{k,k^{\prime}}^{*}$. Since
$P^{*}\notin M$ implies $\bm{\theta}\neq\bm{\theta}^{*}$ for those models,
from the Gibbs inequality, this value is strictly larger than $H(P^{*})$. The
difference between the first terms gets dominant since the parametric
complexity as well as the codelength to encode functions divided by $n$
converges to 0 for each model [13]. It then follows that the probability of
having
$\mathcal{L}^{d}(z^{n};M)>\mathcal{L}^{d}(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
tends to one, which implies the consistency of CLOUD.

##### In case of $P^{*}\in M_{X\to Y}$ or $P^{*}\in M_{X\leftarrow Y}$:

By the symmetry, we restrict ourselves to the case of $P^{*}\in M_{X\to Y}$
without loss of generality. Let $\bm{\theta}^{*}$ be a parameter such that
$P^{*}(X,Y)=P(X,Y;\bm{\theta}^{*})$. The first term of Eq. (2) for both
$M_{X\to Y}$ and $M_{X\leftarrow C\to Y}$ converges to $-\log
P(z^{n};\bm{\theta}^{*})+O_{P}(1)$. As for the second term, we see
$\log\mathcal{C}_{n}(M_{X\to Y})<\log\mathcal{C}_{n}(M_{X\leftarrow C\to Y})$,
which leads to

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\leftarrow C\to
Y})-\mathcal{L}^{d}(z^{n};M_{X\to Y})$
$\displaystyle=\log\mathcal{C}_{n}(M_{X\leftarrow C\to
Y})-\log\mathcal{C}_{n}(M_{X\to Y})+O_{P}\left(1\right)$
$\displaystyle=\Omega(\log n)+O_{P}(1).$

Therefore, the probability that $M_{X\leftarrow C\to Y}$ achieves the smallest
codelength converges to 0. As discussed in the case above, the negative log-
likelihood function in the first term in Eq. (2) divided by $n$ converges to
strictly larger value than $H(P^{*})$ if
$M\in\left\\{M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},M_{X\leftarrow
Y}\right\\}$. The difference between the first terms gets dominant as
mentioned above. Hence, the probablity that $\mathcal{L}^{d}(z^{n};M_{X\to
Y})$ is shortest converges to 1.

##### In case of $P^{*}\in M_{X\leftarrow C\to Y}\setminus\left(M_{X\leftarrow
Y}\cup M_{X\leftarrow Y}\cup\cup
M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\right)$:

As discussed in the case above, the negative log-likelihood function in the
first term in Eq. (2) divided by $n$ converges to strictly larger value than
$H(P^{*})$ if $M\neq M_{X\leftarrow C\to Y}$. Since the first term is dominant
as mentioned in above, the probability that $M_{X\leftarrow C\to Y}$ will be
selected converges to 1 as $n\to\infty$. ∎

###  5.2 Continuous Case

######  Theorem 2.

Let $\mathcal{P}$ be a finite set so that
$\cup_{(m_{X},m_{Y})\in\mathcal{P}}\textsf{HIS}^{m_{X},m_{Y}}=\textsf{HIS}^{\bar{m}_{X},\bar{m}_{Y}}$
holds for some $(\bar{m}_{X},\bar{m}_{Y})\in\mathcal{P}$. For true
distribution $p^{*}$, we define
$p^{*}_{m_{X},m_{Y}}\in\textsf{HIS}^{m_{X},m_{Y}}$ as follows:

$p^{*}_{m_{X},m_{Y}}(x,y)=\frac{\iint_{(x^{\prime},y^{\prime})\in
I(x,y)}p^{*}(x^{\prime},y^{\prime})dx^{\prime}dy^{\prime}}{m_{X}m_{Y}},$

where $I(x,y)=I_{k}^{X}\times I_{k^{\prime}}^{Y}$ such that $x\in I_{k}^{X}$
and $y\in I_{k^{\prime}}^{Y}$ holds. Then we define true model $M^{*}(p^{*})$
as follows:

$\displaystyle
M^{*}(p^{*})=\begin{cases}M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}&p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in
M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\\\ M_{X\to
Y}&p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M_{X\to Y}\\\ M_{X\leftarrow
Y}&p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M_{X\leftarrow Y}\\\ M_{X\leftarrow C\to
Y}&p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M_{X\leftarrow C\to
Y}\setminus\left(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\cup M_{X\leftarrow
Y}\cup M_{X\leftarrow Y}\right).\end{cases}$

Then, the probability that CLOUD outputs $M^{*}(p^{*})$ using $\mathcal{P}$
and $n$ i.i.d. samples from $p^{*}$ converges to 1 in the limit of
$n\to\infty$, provided that the maximum likelihood estimation of $\hat{f}$ and
$\hat{g}$ is successful.

###### Proof.

From the definition of $\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$ in Eq. (13), we
see

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$
$\displaystyle=\mathcal{L}^{d}(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y});\textnormal{DISC}(M;m_{X},m_{Y}))-n\log
m_{X}m_{Y}+\textnormal{const.}$
$\displaystyle=\min_{P\in\textrm{DISC}(M;m_{X},m_{Y})}-\log
P(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y}))+\mathcal{C}_{n}(\textrm{DISC}(M;m_{X},m_{Y}))$
$\displaystyle\phantom{\min_{P\in\textrm{DISC}(M;m_{X},m_{Y})}x}-n\log
m_{X}m_{Y}+\textnormal{const.}$ (19)

For fixed $\hat{P}\in\textnormal{DISC}(M,m_{X},m_{Y})$, there is a
corresponding density function $\hat{p}\in
M^{m_{X},m_{Y}}\subset\textsf{HIS}^{m_{X},m_{Y}}$ such that
$m_{X}m_{Y}\hat{P}(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y}))=\hat{p}(x^{n},y^{n})$
for all $x^{n}$ and $y^{n}$. Using this correspondence, we see the following
expansion in probability:

$\displaystyle-\log\hat{P}(\textnormal{disc}(x^{n};m_{X}),\textnormal{disc}(y^{n};m_{Y}))-n\log
m_{X}m_{Y}$ $\displaystyle=-\log\hat{p}(z^{n})$
$\displaystyle=nH(p^{*})+n\textnormal{KL}(p^{*}\|\hat{p})+O_{P}(1),$ (20)

where $H(p)=\iint p(x,y)\log\frac{1}{p(x,y)}dxdy$ and
$\mathrm{KL}(p||p^{\prime})=\iint
p(x,y)(\log{p(x,y)}-\log{p^{\prime}(x,y)})dxdy$. We thus obtain

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$
$\displaystyle=nH(p^{*})+n\min_{\hat{p}\in
M^{m_{X},m_{Y}}}\mathrm{KL}(p^{*}\mathop{\|}\hat{p})$
$\displaystyle+\mathcal{C}_{n}(\mathrm{DISC}(M;m_{X},m_{Y}))+O_{P}(1).$

From the assumption on $\mathcal{P}$, we further obtain

$\displaystyle\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$
$\displaystyle=nH(p^{*})+n\min_{(m_{X},m_{Y})\in\mathcal{P}}\min_{\hat{p}\in
M^{m_{X},m_{Y}}}\mathrm{KL}(p^{*}\mathop{\|}\hat{p})+\mathcal{C}_{n}(\mathrm{DISC}(M;m_{X},m_{Y}))+O_{P}(1)$
$\displaystyle=nH(p^{*})+n\min_{\hat{p}\in
M^{\bar{m}_{X},\bar{m}_{Y}}}\mathrm{KL}(p^{*}\mathop{\|}\hat{p})+\mathcal{C}_{n}(\mathrm{DISC}(M;\bar{m}_{X},\bar{m}_{Y}))+O_{P}(1).$

For any $\hat{p}\in\mathsf{HIS}^{m_{X},m_{Y}}$, we see

$\displaystyle\mathrm{KL}(p^{*}\mathop{\|}\hat{p})$
$\displaystyle=\sum_{k,k^{\prime}}\iint_{(x,y)\in I^{X}_{k}\times
I^{Y}_{k^{\prime}}}p^{*}(x,y)\left(\log\frac{p^{*}(x,y)}{p^{*}_{m_{X},m_{Y}}(x,y)}+\log\frac{p^{*}_{m_{X},m_{Y}}(x,y)}{\hat{p}(x,y)}\right)dxdy$
$\displaystyle=\mathrm{KL}(p^{*}\mathop{\|}p^{*}_{m_{X},m_{Y}})+\sum_{k,k^{\prime}}\iint_{(x,y)\in
I^{X}_{k}\times
I^{Y}_{k^{\prime}}}p^{*}(x,y)\left(\log\frac{p^{*}_{m_{X},m_{Y}}(x,y)}{\hat{p}(x,y)}\right)dxdy$
$\displaystyle=\mathrm{KL}(p^{*}\mathop{\|}p^{*}_{m_{X},m_{Y}})+\sum_{k,k^{\prime}}\iint_{(x,y)\in
I^{X}_{k}\times
I^{Y}_{k^{\prime}}}p^{*}_{m_{X},m_{Y}}(x,y)\left(\log\frac{p^{*}_{m_{X},m_{Y}}(x,y)}{\hat{p}(x,y)}\right)dxdy$
$\displaystyle=\mathrm{KL}(p^{*}\mathop{\|}p^{*}_{m_{X},m_{Y}})+\mathrm{KL}(p^{*}_{m_{X},m_{Y}}\mathop{\|}\hat{p}).$

The third equality holds since $\left(\log
p^{*}_{m_{X},m_{Y}}(x,y)-\log\hat{p}(x,y)\right)$ is constant for all
$(x,y)\in I^{X}_{k}\times I^{Y}_{k^{\prime}}$ and $\iint_{(x,y)\in
I^{X}_{k}\times I^{Y}_{k^{\prime}}}p^{*}(x,y)dxdy=\iint_{(x,y)\in
I^{X}_{k}\times I^{Y}_{k^{\prime}}}p^{*}_{m_{X},m_{Y}}(x,y)dxdy$. Therefore,
$\hat{p}^{*}_{m_{X},m_{Y}}$ is a unique minimizer of
$\mathrm{KL}(p^{*}\mathop{\|}\hat{p})$, which implies that

$\displaystyle\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)$
$\displaystyle=nH(p^{*})+n\mathrm{KL}(p^{*}\mathop{\|}p^{*}_{\bar{m}_{X},\bar{m}_{Y}})$
$\displaystyle+\mathcal{C}_{n}(\mathrm{DISC}(M;\bar{m}_{X},\bar{m}_{Y}))+O_{P}(1).$

for $M$ such that $p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M$. If
$p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\notin M$ otherwise, we see

$\displaystyle\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)-\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M_{X\leftarrow
C\to Y})$
$\displaystyle=\min_{(m_{X},m_{Y})\in\mathcal{P}}\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M)-nH(p^{*})-n\mathrm{KL}(p^{*}\mathop{\|}p^{*}_{\bar{m}_{X},\bar{m}_{Y}})+o_{P}(n)$
$\displaystyle=\min_{\hat{p}\in
M^{\bar{m}_{X},\bar{m}_{Y}}}n\mathrm{KL}(p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\mathop{\|}\hat{p})+o_{P}(n),$

since it always holds $p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M_{X\leftarrow C\to
Y}$. This implies that the probability that $M$ is chosen converges to 0. The
consistency among models in which $p^{*}_{\bar{m}_{X},\bar{m}_{Y}}\in M$
follows similarly as that of CLOUD in the discrete case. ∎

##  6 Experiment

In this section, we demonstrate that 1) our proposed method CLOUD can solve
the Reichenbach problem in situations where it is difficult to make
assumptions about unobserved common causes, and 2) CLOUD shows higher
inference accuracy in identifying causal relationships compared to existing
methods, even when the true data-generating mechanism violates the assumptions
of our method.

In the first experiment, we designed the Reichenbach problem scenarios with
all set of synthetic data types — discrete, mixed, and continuous, and
verified that our proposed method exhibits high accuracy and consistency in
solving the problem. In particular, it effectively detects the presence of
unobserved common causes $C$ even in situations where the observed variables
$X,Y$ are generated from complex mechanisms $f(X,C),g(Y,C)$.

In the second experiment, we compared the performance of our proposed method
against existing methods in identifying causal relationships in synthetic data
generated from either $M_{X\to Y}$ or $M_{X\leftarrow C\to Y}$, and
demonstrated its effectiveness.

The third experiment tested the ability of our proposed method to determine
the directions of causality and detect unobserved common causes in real-world
data generated from unknown and complex data-generating process.

We implemented CLOUD in Python and provide the source code at
https://github.com/Matsushima-lab/CLOUD.

###  6.1 Consistency of CLOUD on the Reichenbach problem

In the first experiment, we verified the performance of CLOUD on the
Reichenbach problem with synthetic data and confirmed its consistency in model
selection. We randomly selected each SCM corresponding to the four causal
relationships, and then generated data $z^{n}$ with sample size
$n=10^{2},10^{3},10^{4}$. The data-generating processes for each combination
were defined as follows:

##### Discrete Case:

$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$:
$X$ and $Y$ were independently generated from categorical distributions.
$M_{X\leftarrow C\to Y}$:
$C\in\left\\{0,1,\ldots,99\right\\}$ was independently generated from a
categorical distribution, and $X$ and $Y$ were set to the quotient and
remainder of $C$ divided by 10, respectively. $M_{X\to Y}$:
$X$ and $E_{Y}$ were independently generated from categorical distributions,
$f$ was generated uniformly randomly from all non-constant functions, and
subsequently $Y$ was set to $Y=f(X)+E_{Y}(\bmod\ 10)$. The same applied to the
case of $M_{X\leftarrow Y}$.

##### Mixed Case ($X$ is continuous and $Y$ is discrete):

$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$:
$X$ was independently generated from a Gaussian distribution, and $Y$ from a
categorical distribution. $M_{X\leftarrow C\to Y}$:
$C\in\left\\{0,1,\ldots,99\right\\}$ was independently generated from a
categorical distribution. $X$ and $Y$ were then generated as follows:
$X=b\sin{C}+E_{X},Y=\lfloor\frac{C}{10}\rfloor$, where $b$ was sampled from a
uniform distribution $\mathcal{U}(2,4)$ and $E_{X}$ follows a Gaussian
distribution $\mathcal{N}(0,0.1^{2})$. $M_{X\to Y}$:
$X$ was generated from a mixture of Gaussian distributions with three clusters
as:

$p(X)=0.6\cdot\mathcal{N}(X;-5,2^{2})+0.2\cdot\mathcal{N}(X;0,1^{2})+0.2\cdot\mathcal{N}(X;5,2^{2}).$

Then, $X$ was divided into $m_{X}\sim\text{Uniform}\\{2,3,4\\}$ equal
intervals, with $f(X)\sim\text{Uniform}\\{0,1,\ldots,10\\}$ assigned to each
interval. Additive noise $E_{Y}\sim\text{Uniform}\\{-1,0,1\\}$ was added to
generate $Y$, with addition over $\mathbb{R}/10\mathbb{Z}$. The correlation
coefficient was ensured to be greater than 0.2. $M_{X\leftarrow Y}$:
$Y$ was generated from an $m_{Y}$-valued categorical distribution, and then
$X$ was set to $X=2Y+3\sin{Y}+E_{Y}$, where $m_{Y}$ was generated from uniform
distribution $\text{Uniform}\\{2,3,\ldots,8\\}$ and
$E_{Y}\sim\mathcal{N}(0,1)$, with addition taken over
$\mathbb{R}/20\mathbb{Z}$.

##### Continuous Case:

$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$:
$X$ and $Y$ were independently generated from Gaussian distributions.
$M_{X\leftarrow C\to Y}$:
$X$ and $Y$ were generated based on an ellipse equation with eccentricity $e$
and semi-major axis $a$: $r=\frac{a(1-e^{2})}{1+e\cos(\phi)}$ with
$0\leq\phi<2\pi$ In this process, we sampled $e$ from $\mathcal{U}(0.5,0.9)$
and $a$ from $\text{Uniform}\\{1,2,3\\}$. $X$ and $Y$ were then set to
$X=r\cos(\phi+\eta)+E_{X},Y=r\sin(\phi+\eta)+E_{Y}(0\leq\phi\leq 2\pi)$ ,
where $\eta$ were sampled from $\mathcal{U}(\pi/4,\pi/3)$, and $E_{X}$ and
$E_{Y}$ follow $\mathcal{N}(0,(0.1a)^{2})$. $M_{X\to Y}$:
$X$ was generated from a probability density function of a mixture Gaussian
distribution with three clusters. Then $Y$ was set to
$Y=a*\text{disc}(X,m_{X})+b+E_{Y}$, where
$m_{X}\sim\text{Uniform}\\{2,3,4\\},a\sim\mathcal{U}(4,7),b\sim\mathcal{U}(1,5),E_{Y}\sim\mathcal{N}(0,1)$,
with addition taken over $\mathbb{R}/20\mathbb{Z}$. The procedure for the
$M_{X\leftarrow Y}$ was analogous.

Table 1: Results on experiment 1 $n$ | $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ | $M_{X\leftarrow C\to Y}$ | $M_{X\to Y}$ | $M_{X\leftarrow Y}$
---|---|---|---|---
2-13 | disc. | mix. | cont. | disc. | mix. | cont. | disc. | mix. | cont. | disc. | mix. | cont.
$10^{2}$ | 95.1 | 91.9 | 90.4 | 85.0 | 96.1 | 95.2 | 24.6 | 91.4 | 74.8 | 26.8 | 77.4 | 75.4
$10^{3}$ | 100 | 99.8 | 96.3 | 87.6 | 100 | 100 | 99.6 | 100 | 98.4 | 99.9 | 99.4 | 98.9
$10^{4}$ | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 100 | 99.7 | 100 | 100 | 99.9

Table 1 shows the transition of accuracy as a fraction of correct inference
with sample size $n$. Figures 1, 2, and 3 visualize the inference results of
CLOUD as confusion matrices for each data type. These indicate that accuracy
improves as $n$ increases. In all cases, the accuracy reaches $\sim$100% at
$n=10000$. We thus empirically observed the consistency of CLOUD.

Figure 1: Confusion matrices in the Discrete Case of Experiment 1 Figure 2:
Confusion matrices in the Mixed Case of Experiment 1 Figure 3: Confusion
matrices in the Continuous Case of Experiment 1

CLOUD calculates the codelength of the observed data for each causal model and
selects the one that achieves the shortest codelength. Therefore, we can
expect that our method is more confident in its inference when the difference
in the codelengths between the shortest one and the rest is larger. The next
experiment examined whether the difference in the codelengths per sample size
of the shortest and the second shortest model, denoted as $\Delta$, can be
interpreted as the confidence of CLOUD. We generated 1000 synthetic datasets
from discrete causal models and calculated the accuracy at each decision rate
$d$. The accuracy at decision rate $d$ is defined as the accuracy at the upper
$d$% of datasets when datasets are sorted in descending order of $\Delta$. The
result is shown in Fig 4. For each model, the accuracy is higher when the
decision rate is smaller, i.e., $\Delta$s are larger. We thus conclude that
$\Delta$ is interpreted as the confidence of the inference in CLOUD.

Figure 4: Accuracy vs. decision rate of CLOUD on synthetic data

###  6.2 Comparison to existing methods in case of $X\to Y$ and $X\leftarrow
C\to Y$

In the second experiment, we compared the accuracy of out proposed method
CLOUD with existing methods in two scenarios where the ground truth of the
causal relationship is either $X\to Y$ or $X\leftarrow C\to Y$. This
comparison aimed to evaluate CLOUD’s performance in identifying the direction
of causality and in detecting unobserved common causes. We set the sample size
at $n=500$. For each SCM, synthetic datasets $z^{n=500}=x^{n=500}\times
y^{n=500}$ were generated 100 times. The accuracy was determined by
calculating the proportion of correctly identified causal relationships across
these iterations.

To verify the cases where the true data-generating process violates the
assumptions of our models, for $M_{X\to Y}$, we employ a non-cyclic ANM in
which addition is taken without modulo operation: Discrete Case
$X$ and $E_{Y}$ were randomly and independently generated from categorical
distributions, and $Y$ was set to $Y=f(X)+E_{Y}$ with a mapping function $f$
that was also randomly set. Mixed Case [$X$ is continuous and $Y$ is discrete]
$X$ was generated from $\mathcal{N}(0,10^{2})$, and then $X$ was divided into
$m_{X}\sim\text{Uniform}\\{2,3,4\\}$ equal intervals. For each interval, a
value $f(X)\sim\text{Uniform}\\{1,2,\ldots,24\\}$ was randomly assigned, and
additive noise $E_{Y}\sim\text{Uniform}\\{-1,0,1\\}$ was added to generate
$Y$. Continuous Case
$X$ was generated from a three-class mixture Gaussian distribution as in
Experiment 1, and $Y$ was generated from $Y=af(X)+b\sin{2\pi X}+E_{Y}$, where
$a,b\sim\mathcal{U}([-2,-0.5]\cup[0.5,2])$, and $E_{Y}\sim\mathcal{N}(0,1)$.
For the function $f(X)$, we considered two cases: one linear and the other a
cubic function $x^{3}$.

In the case where $M_{X\leftarrow C\to Y}$ is the ground truth, the same
process used in Experiment 1 was used.

For existing methods, we employed ECI [11], DR [23], CISC [2], ACID [3], LiM
[39], HCM [14], LiNGAM [33], ANM [24], RCD [16], CAMUV [17], BUPL [37], and
COCA [9]. These methods are categorized into those that can detect
$M_{X\leftarrow C\to Y}$ and those that cannot. While ECI, DR, CISC, ACID,
LiM, HCM, LiNGAM, and ANM distinguish between $M_{X\to Y}$ and $M_{X\leftarrow
Y}$, RCD, CAMUV, and BUPL infer a causal model from a set of model candidates
including $M_{X\leftarrow C\to Y}$. Moreover, COCA selects a model only from
$\mathcal{M}=\\{M_{X\to Y},M_{X\leftarrow C\to Y}\\}$. We note that our
proposed method CLOUD as well as LiM and HCM can accept all type of data,
wheres others specialize for either discrete or continuous data.

We utilized the implementations of HCM by Li et al. (2022), COCA by Kaltenpoth
et al. (2019), and others by Ikeuchi et al. (2023). Default hyper parameter
values were used.

Results are shown in Table 2. Unlike existing methods, CLOUD is applicable to
all experimental conditions. In particular, CLOUD is the first method capable
of detecting unobserved common causes in discrete and mixed cases. As Table 2
demonstrates, CLOUD showed consistently high inference accuracy across all
cases, regardless of the data type, even though the number of model candidates
of CLOUD is as many as 4 models. Especially, CLOUD outperformed previous
methods in the discrete case.

Table 2: Results on experiment 2

Methods | Direct Case | Confounded Case | $|\mathcal{M}|$ | $C$
---|---|---|---|---
2-8 | disc. | mix. | cont. | disc. | mix. | cont.
4-5 |  |  | linear | cubic |  |  |  |  |
ECI | 89% | - | - | - | - | - | - | 2 | -
DR | 85% | - | - | - | - | - | - | 2 | -
CISC | 96% | - | - | - | - | - | - | 2 | -
ACID | 92% | - | - | - | - | - | - | 2 | -
LiM | 86% | 12% | 82% | 87% | - | - | - | 3 | -
HCM | 90% | 100% | 85% | 100% | - | - | - | 3 | -
LiNGAM | - | - | 95% | 61% | - | - | - | 3 | -
ANM | - | - | 91% | 100% | - | - | - | 2 | -
RCD | - | - | 48% | 5% | - | - | 96% | 4 | $\ast$
CAMUV | - | - | 92% | 99% | - | - | 100% | 4 | $\ast$
BUPL | - | - | 37% | 8% | - | - | 100% | 4 | $\ast$
COCA | - | - | 2% | 34% | - | - | 13% | 2 | $\ast$
CLOUD | 98% | 99% | 96% | 99% | 88% | 100% | 100% | 4 | $\ast$

Performance comparison of CLOUD against existing methods w.r.t. accuracy in
the discrete case (disc.), mixed case (mix.) and continuous case (cont.) based
on synthetic data generated either from direct case or confounded case.
$|\mathcal{M}|$ column denotes the number of model candidates each method
considers, and $C$ one represents whether each method allows for the existence
of unobserved common causes or not ($\ast$: Yes, -: No)

###  6.3 Real World Data

####  6.3.1 Direct case: Tübingen Benchmark Pairs

In this section, we examined the effectiveness of CLOUD of inferring direct
causality in data generated from complex causal mechanisms by real-world
datasets.

We employed datasets in various application fields from the Tübingen Cause-
Effect Pairs Database [19], which provides a collection of datasets for
testing causal discovery methods. The database contains datasets with known
ground truth to distinguish between cause and effect variables. We note that
the ground truth does not mean there are no unobserved counfounders in
general, except for dataset No.101 which was explicitly generated in an
unconfounded experimental environment. The statistical information of the
datasets used in the experiments is shown in Table 3, and scatter plots for
each data pair are presented in Figure 5. Descriptions for each data pair are
given in Appendix D.

Table 3: Characteristics of Tübingen Cause-Effect-Pairs Dataset | Data-type | $n$ | Ground Truth | $X$ | $Y$
---|---|---|---|---|---
No. 47 | disc. | 254 | $X\leftarrow Y$ | number of cars | working days or not
No. 68 | disc. | 498 | $X\leftarrow Y$ | bytes sent at minute | open http connections
No. 107 | disc. | 240 | $X\to Y$ | contrast | answer correct or not
No. 85 | mix. | 994 | $X\to Y$ | day | protein content of the milk
No. 95 | mix. | 9504 | $X\to Y$ | hour of the day | total electricity consumption
No. 99 | mix. | 2287 | $X\leftarrow Y$ | language test score | social-economic status
No. 23 | cont. | 452 | $X\to Y$ | age | weight
No. 77 | cont. | 8401 | $X\leftarrow Y$ | daily average temperature | solar radiation
No. 101 | cont. | 300 | $X\to Y$ | grey value of a pixel | light intensity
Figure 5: Scatter plots of the Tübingen Cause-Effect Pairs. The horizontal
axis represents $X$, while the vertical axis represents $Y$. Each plot
corresponds to a dataset pair.

We determined the data type for each data pair based on the information of the
variables and run applicable causal discovery methods. Results are shown in
Table 4. CLOUD demonstrated an excellent ability to determine the causal
directions across various data types. Notable, it correctly identified the
causal directions in every case both in mixed and continuous cases. While LiM
and HCM are applicable across all cases, CLOUD achieved the highest number of
correct answers.

In discrete case, CLOUD showed performance comparable to CISC, a state-of-the-
art causal discovery method for discrete data. CLOUD also correctly inferred
the right directions with high confidence (large $\Delta$) for cases No.47 and
No.68, where LiM and HCM were incorrect. Notably, in case No.107, CLOUD
detected the presence of a confounding factor (indicated as ‘conf’) rather
than a clear causal direction. This is a crucial feature in real-world data
analysis involving potential confounding variables.

In continuous case, CLOUD successfully identified the causal directions across
all pairs, despite having four model candidates. In contrast, RCD, CAMUV, and
BUPL, which also have $|\mathcal{M}|=4$, frequently resulted in the ‘conf’
(confounded) classification. This tendency suggests a bias in these methods
towards indicating confoundedness rather than directly identifying causality.
Particularly, it would be incorrect to conclude the presence of confounders in
dataset No.101 mentioned above.

Overall, the results from the Tübingen Benchmark Pairs suggests that CLOUD is
a reliable causal discovery method to identify the true causal direction,
particularly in the situations involving complex data-generating process
across all data types.

Table 4: Results on Tübingen Benchmark Pairs Dataset

Methods | Discrete Case | Mixed Case | Continuous Case
---|---|---|---
2-10 |  No.47 |  No.68 |  No.107 |  No.85 |  No.95 |  No.99 |  No.23 |  No.77 |  No.101
ECI |  $\checkmark$ |  $\times$ |  $\times$ |  - |  - |  - |  - |  - |  -
DR |  $\thickapprox$ |  $\thickapprox$ |  $\thickapprox$ |  - |  - |  - |  - |  - |  -
CISC |  $\checkmark$ |  $\checkmark$ |  $\times$ |  - |  - |  - |  - |  - |  -
ACID |  $\times$ |  $\times$ |  $\thickapprox$ |  - |  - |  - |  - |  - |  -
LiM |  $\times$ |  $\times$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\times$ |  $\times$
HCM |  $\times$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$ |  $\times$ |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
LiNGAM |  - |  - |  - |  - |  - |  - |  $\times$ |  $\checkmark$ |  $\checkmark$
ANM |  - |  - |  - |  - |  - |  - |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
RCD |  - |  - |  - |  - |  - |  - |  conf |  conf |  conf
CAMUV |  - |  - |  - |  - |  - |  - |  $\checkmark$ |  conf |  conf
BUPL |  - |  - |  - |  - |  - |  - |  conf |  conf |  conf
COCA |  - |  - |  - |  - |  - |  - |  $\checkmark$ |  $\checkmark$ |  $\checkmark$
CLOUD |  $\boldsymbol{\checkmark}$
$(\Delta=0.16)$ |  $\boldsymbol{\checkmark}$
$(\Delta=1.5)$ |  conf
$(\Delta=0.01)$ |  $\boldsymbol{\checkmark}$
$(\Delta=0.03)$ |  $\boldsymbol{\checkmark}$
$(\Delta=0.03)$ |  $\boldsymbol{\checkmark}$
$(\Delta=0.09)$ |  $\boldsymbol{\checkmark}$
$(\Delta=0.10)$ |  $\boldsymbol{\checkmark}$
$(\Delta=0.21)$ |  $\ \boldsymbol{\checkmark}$
$(\Delta=0.12)$

$\checkmark$ indicates that a method inferred the true causal direction.
$\times$ indicates that the output of a method was wrong direction.
$\thickapprox$ indicates that a method drew undisicive conclusion. conf
indicates that a method inferred $M_{X\leftarrow C\to Y}$.

####  6.3.2 Confounded case: SOS DNA Repair Network Dataset

Finally, we tested CLOUD’s ability of detecting latent confounding variables
in complex causal relationships using SOS DNA repair network in E.coli [30].
This dataset describes the causal relationships at the protein level between
genes, consisting of measurements for eight different genes under four
distinct ultraviolet radiation conditions, with a total sample size of
$n=200$. A ground truth network, as established by [21], is depicted in Fig.
6. The gene lexA has causal influence on all other genes, creating a situation
where at least lexA is an unobserved common cause among variables downstream
(children) of lexA. Therefore, we randomly selected pairs of child nodes of
lexA and conducted experiments to detect the presence of an unobserved common
cause (lexA) between each pair. We also note that while the experimental setup
correctly represent confounded case, the correct directed cases could not be
extracted, such as the arrow from lexA to umuDC. This is because there might
be common causes between the two.

lexArecAumuDCuvrDuvrAuvrYruvApolB Figure 6: Ground truth graph of SOS DNA
repair network Figure 7: Scatter plots of variable pairs from the SOS DNA
repair network dataset

For comparison, we employed RCD, CAMUV, BUPL, and COCA, which are capable of
inferring $M_{X\leftarrow C\to Y}$. Since COCA considers asymmetry between $X$
and $Y$, we conducted experiments in two settings: one with
$\mathcal{M}=\\{M_{X\to Y},M_{X\leftarrow C\to Y}\\}$ (referred to as
COCA($X\to Y$)) and another with $\mathcal{M}=\\{M_{X\leftarrow
Y},M_{X\leftarrow C\to Y}\\}$ (referred to as COCA($X\leftarrow Y$)).

The results are shown in Table 5. For the pairs (polB, umuDC) and (uvrD,
uvrA), the methods CLOUD, RCD, and CAMUV successfully identified the presence
of unobserved common causes. However, for the pair (uvrY, ruvA), only COCA
detected the unobserved common cause. We observed that CLOUD inferred causal
independence for this pair, likely due to ruvA exhibiting zero-inflation and
the statistical independence of the pair.

CLOUD demonstrated performance comparable to existing state-of-the-art methods
that allow for unobserved common causes.

Given the prevalence of unobserved confounders in real-world applications,
which can lead to incorrect causal conclusions if not properly addressed, we
can conclude that CLOUD is equipped with a crucial feature for handling such
scenarios. Moreover, considering the results from the Tübingen Benchmark Pairs
as well, as presented in Table 4, we further affirm that CLOUD is a reliable
method for real-world data analysis.

Table 5: Results on SOS DNA repair network

Ground Truth |  CLOUD | RCD | CAMUV | BUPL | COCA($X\to Y$) | COCA($X\leftarrow Y$)
---|---|---|---|---|---|---
polB $\leftarrow C\to$ umuDC |  ✓ $(\Delta=0.21)$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | $\times$
uvrD $\leftarrow C\to$ uvrA |  ✓ $(\Delta=0.013)$ | $\checkmark$ | $\checkmark$ | $\checkmark$ | $\times$ | $\times$
uvrY $\leftarrow C\to$ ruvA |  $\times$ $(\Delta=0.25)$ | $\times$ | $\times$ | $\times$ | $\checkmark$ | $\checkmark$

$\checkmark$ indicates that a method inferred $M_{X\leftarrow C\to Y}$, while
$\times$ indicates that a method did not.

##  7 Conclusion

This paper proposed CLOUD, a novel causal discovery method for causal
relationships between two variables with unobserved common causes across
discrete, mixed, and continuous data types.

Based on the Reichenbach’s common cause principle, we defined the Reichenbach
problem as a problem to statistically infer the causal relationships among
four models: one independent model, one with unobserved common cause, and two
models with direct causality.

By employing the NML code, CLOUD offers a model selection approach to solve
the Reichenbach problem without relying on assumptions about unobserved common
causes. CLOUD formulates four models for each causal relationship and data-
type. In particular, CLOUD expresses all joint distributions of $X$ and $Y$
under $M_{X\leftarrow C\to Y}$, which enables us to avoid explicitly modeling
$C$. CLOUD calculates the NML-based codelength of the observational data under
those four models and then infers a causal relationship by selecting the
corresponding causal model that achieves the shortest codelength. We
successfully extended the CLOUD from discrete to continuous data, through
discretization. CLOUD has a consistency with respect to selected models and it
is theoretically proven.

Through both synthetic and real-world data experiments, CLOUD has proven its
effectiveness in solving the Reichenbach problem with high accuracy and
consistency. It stands out in its ability to identify causal directions with
greater precision than existing methods, across a variety of data types and
under complex data-generating conditions. Additionally, CLOUD has demonstrated
a strong performance in detecting latent variables, showcasing its robustness
and reliability in causal discovery.

However, challenges remain, as evidenced by its performance on zero-inflated
data in our final experiment where CLOUD mistakenly determined they are
causally independent. This implies the applicable range of CLOUD is still
restricted, despite its ability to detect unobserved common causes without
additional assumptions on it. Moreover, Additive Noise Model which we assume
has known to be vulnerable in handling data with heteroscedastic noise, i.e
noise variances are dependent of observed variables unlike ANM’s assumption.
Recent research is actively addressing these challenges [36, 38, 4].

Future work aims to broaden CLOUD’s scope, addressing its current SCM
assumptions to enhance robustness and applicability across diverse data
scenarios.

## References

  *   * Akaike [1998] Akaike H (1998) Information theory and an extension of the maximum likelihood principle. In: Selected papers of hirotugu akaike. Springer, p 199–213
  * Budhathoki and Vreeken [2017] Budhathoki K, Vreeken J (2017) MDL for causal inference on discrete data. In: 2017 IEEE International Conference on Data Mining (ICDM), pp 751–756
  * Budhathoki and Vreeken [2018] Budhathoki K, Vreeken J (2018) Accurate causal inference on discrete data. In: 2018 IEEE International Conference on Data Mining (ICDM), pp 881–886
  * Choi and Ni [2023] Choi J, Ni Y (2023) Model-based causal discovery for zero-inflated count data. Journal of Machine Learning Research 24(200):1–32
  * Hirai and Yamanishi [2013] Hirai S, Yamanishi K (2013) Efficient computation of normalized maximum likelihood codes for gaussian mixture models with its applications to clustering. IEEE Transactions on Information Theory 59(11):7718–7727
  * Hoyer et al [2008a] Hoyer P, Janzing D, Mooij JM, et al (2008a) Nonlinear causal discovery with additive noise models. Advances in neural information processing systems 21
  * Hoyer et al [2008b] Hoyer PO, Shimizu S, Kerminen AJ, et al (2008b) Estimation of causal effects using linear non-gaussian causal models with hidden variables. International Journal of Approximate Reasoning 49(2):362–378
  * Janzing and Schölkopf [2010] Janzing D, Schölkopf B (2010) Causal inference using the algorithmic Markov condition. IEEE Transactions on Information Theory 56(10):5168–5194
  * Kaltenpoth and Vreeken [2019] Kaltenpoth D, Vreeken J (2019) We are not your real parents: Telling causal from confounded using mdl. In: Proceedings of the 2019 SIAM International Conference on Data Mining, SIAM, pp 199–207, URL https://github.com/davidkwca/CoCa
  * Kobayashi et al [2022] Kobayashi M, Miyaguchi K, Matsushima S (2022) Detection of unobserved common cause in discrete data based on the mdl principle. In: 2022 IEEE International Conference on Big Data (Big Data), IEEE, pp 45–54
  * Kocaoglu et al [2017] Kocaoglu M, Dimakis AG, Vishwanath S, et al (2017) Entropic causal inference. In: Thirty-First AAAI Conference on Artificial Intelligence
  * Kontkanen and Myllymäki [2007] Kontkanen P, Myllymäki P (2007) A linear-time algorithm for computing the multinomial stochastic complexity. Information Processing Letters 103(6):227–233
  * Kontkanen et al [2008] Kontkanen P, Wettig H, Myllymäki P (2008) NML computation algorithms for tree-structured multinomial Bayesian networks. EURASIP Journal on Bioinformatics and Systems Biology 2007:1–11
  * Li et al [2022] Li Y, Xia R, Liu C, et al (2022) A hybrid causal structure learning algorithm for mixed-type data. In: Proceedings of the AAAI Conference on Artificial Intelligence, pp 7435–7443, URL https://github.com/DAMO-DI-ML/AAAI2022-HCM
  * Liu and Chan [2016] Liu F, Chan L (2016) Causal inference on discrete data via estimating distance correlations. Neural computation 28(5):801–814
  * Maeda and Shimizu [2020] Maeda TN, Shimizu S (2020) Rcd: Repetitive causal discovery of linear non-gaussian acyclic models with latent confounders. In: International Conference on Artificial Intelligence and Statistics, PMLR, pp 735–745
  * Maeda and Shimizu [2021] Maeda TN, Shimizu S (2021) Causal additive models with unobserved variables. In: Uncertainty in Artificial Intelligence, PMLR, pp 97–106
  * Marx and Vreeken [2021] Marx A, Vreeken J (2021) Formally justifying mdl-based inference of cause and effect. arXiv preprint arXiv:210501902
  * Mooij et al [2016] Mooij JM, Peters J, Janzing D, et al (2016) Distinguishing cause from effect using observational data: methods and benchmarks. The Journal of Machine Learning Research 17(1):1103–1204
  * Pearl [2009] Pearl J (2009) Causality. Cambridge university press
  * Perrin et al [2003] Perrin BE, Ralaivola L, Mazurie A, et al (2003) Gene networks inference using dynamic bayesian networks. Bioinformatics-Oxford 19(2):138–148
  * Peters et al [2010] Peters J, Janzing D, Schölkopf B (2010) Identifying cause and effect on discrete data using additive noise models. In: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, pp 597–604
  * Peters et al [2011] Peters J, Janzing D, Scholkopf B (2011) Causal inference on discrete data using additive noise models. IEEE Transactions on Pattern Analysis and Machine Intelligence 33(12):2436–2450
  * Peters et al [2014] Peters J, Mooij JM, Janzing D, et al (2014) Causal discovery with continuous additive noise models. Journal of Machine Learning Research 15(58):2009–2053
  * Peters et al [2017] Peters J, Janzing D, Schölkopf B (2017) Elements of causal inference: foundations and learning algorithms. The MIT Press
  * Rissanen [1978] Rissanen J (1978) Modeling by shortest data description. Automatica 14(5):465–471
  * Rissanen [1983] Rissanen J (1983) A universal prior for integers and estimation by minimum description length. The Annals of statistics 11(2):416–431
  * Rissanen [1989] Rissanen J (1989) Stochastic complexity in statistical inquiry. World Scientific
  * Rissanen [2012] Rissanen J (2012) Optimal Parameter Estimation. Cambridge University Press
  * Ronen et al [2002] Ronen M, Rosenberg R, Shraiman BI, et al (2002) Assigning numbers to the arrows: parameterizing a gene regulation network by using accurate expression kinetics. Proceedings of the national academy of sciences 99(16):10555–10560
  * Schölkopf [2022] Schölkopf B (2022) Causality for machine learning. In: Probabilistic and Causal Inference: The Works of Judea Pearl. p 765–804
  * Shimizu et al [2006] Shimizu S, Hoyer PO, Hyvärinen A, et al (2006) A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research 7(10)
  * Shimizu et al [2011] Shimizu S, Inazumi T, Sogawa Y, et al (2011) Directlingam: A direct method for learning a linear non-gaussian structural equation model. The Journal of Machine Learning Research 12:1225–1248
  * Spirtes et al [2000] Spirtes P, Glymour CN, Scheines R (2000) Causation, prediction, and search. MIT press
  * Stegle et al [2010] Stegle O, Janzing D, Zhang K, et al (2010) Probabilistic latent variable models for distinguishing between cause and effect. Advances in neural information processing systems 23
  * Tagasovska et al [2020] Tagasovska N, Chavez-Demoulin V, Vatter T (2020) Distinguishing cause from effect using quantiles: Bivariate quantile causal discovery. In: International Conference on Machine Learning, PMLR, pp 9311–9323
  * Tashiro et al [2014] Tashiro T, Shimizu S, Hyvärinen A, et al (2014) Parcelingam: A causal ordering method robust against latent confounders. Neural computation 26(1):57–83
  * Xu et al [2022] Xu S, Mian OA, Marx A, et al (2022) Inferring cause and effect in the presence of heteroscedastic noise. In: International Conference on Machine Learning, PMLR, pp 24615–24630
  * Zeng et al [2022] Zeng Y, Shimizu S, Matsui H, et al (2022) Causal discovery for linear mixed data. In: Conference on Causal Learning and Reasoning, PMLR, pp 994–1009

##  Appendix A Optimization of $f$ and $g$ through Likelihood Maximization

The algorithm for estimating the function $\hat{f}$ is shown in the Algorithm
4. In order to compute both terms in Eq. (8), we must estimate
$\hat{f}\colon\mathcal{X}\to\mathcal{Y}$ that achieves a higher likelihood to
calculate a shorter codelength of data under $M_{X\to Y}$. First, we
initialize the function $\hat{f}$ that returns the most frequent $y$ for each
$x\in\mathcal{X}$ (lines 1-3). Subsequently, we iteratively update the
function value for $x\in\mathcal{X}$. At $j$-th step, we update the function
value $f(x)$ while fixing all other mapping $f(x^{\prime})$ for
$x^{\prime}\neq x$ and check if the likelihood increases. If it increases, we
change the value $f(x)$ to new $y$. This update is done for all
$x\in\mathcal{X}$ and repeated until the likelihood no longer increases or for
at most $J$ times. Similarly, the function $g:\mathcal{Y}\to\mathcal{X}$ can
be estimated by replacing $X$ and $Y$ in algorithm 4.

Algorithm 4 Optimize Regression Function $f:\mathcal{X}\to\mathcal{Y}$ with
Likelihood Maximization

1:$z^{n}\in\mathcal{X}^{n}\times\mathcal{Y}^{n}\,$: Dataset with sample size
of $n$

2:$\quad\ J$ : Maximum number of update iterations

3:$\hat{f}\colon\mathcal{X}\to\mathcal{Y}$ : Estimated Regression Function

4:for $x\in\mathcal{X}$ do

5:
$f^{(0)}(x)\leftarrow\mathop{\mathrm{argmax}}\limits_{y\in\mathcal{Y}}n(X=x,Y=y)$

6:end for

7:$j\leftarrow 0$

8:$r_{\mathrm{max}}\leftarrow\max\limits_{\bm{\theta}\in\bm{\Theta}_{X\to
Y}}P\\!\left(z^{n};M_{X\to Y},f^{(0)},\bm{\theta}\right)$

9:while $\mathit{converged}=\mathrm{False}$ or $j<J$ do

10: $j\leftarrow j+1$

11: $\mathit{converged}\leftarrow\mathrm{True}$

12: for $x\in\mathcal{X}$ do

13:
$r\leftarrow\max\limits_{f^{(j-1)}(x)\in\mathcal{Y}}\max\limits_{\bm{\theta}\in\bm{\Theta}_{X\to
Y}}P\\!\left(z^{n};M_{X\to Y},f^{(j-1)},{\bm{\theta}}\right)$ $\triangleright$
subject to $f^{(j-1)}$ is not constant

14: if $r>r_{\mathrm{max}}$ then

15: $r_{\mathrm{max}}\leftarrow r$

16:
$f^{(j)}(x)\leftarrow\mathop{\mathrm{argmax}}\limits_{f^{(j-1)}(x)\in\mathcal{Y}}\max\limits_{\bm{\theta}\in\bm{\Theta}_{X\to
Y}}P\\!\left(z^{n};M_{X\to Y},f^{(j-1)},\bm{\theta}\right)$ $\triangleright$
subject to $f^{(j-1)}$ is not constant

17: $\mathit{converged}\leftarrow\mathrm{False}$

18: end if

19: end for

20:end while

21:return $f^{(j)}$

##  Appendix B Proof of Proposition 1

In this section, we provide the proof for Proposition 1.

###  B.1 Independent Case

We exactly calculate the stochastic complexity defined in Eq. (2) for causal
models $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=\mathcal{SC}(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=-\log
P(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},{\hat{\bm{\theta}}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\left(z^{n}\right))+\log\mathcal{C}_{n}(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$

Using the result of likelihood estimation
$P(X=k,Y=k^{\prime};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},{\hat{\bm{\theta}}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}\left(z^{n}\right))=\frac{n(X=k)}{n}\frac{n(Y=k^{\prime})}{n}$,
the maximum likelihood which the first term on the right-hand side is
represented as:

$\displaystyle\log
P(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},{\hat{\bm{\theta}}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})=\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}+\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y=k^{\prime})\log\frac{n(Y=k^{\prime})}{n},$

where
${\hat{\bm{\theta}}}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}=(\hat{\theta}_{k,k^{\prime}})$
is the maximum likelihood estimator in
$\bm{\Theta}_{{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}}$. The parametric
complexity for $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ as the second term
is calculated as

$\displaystyle\log\mathcal{C}_{n}(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=\log\sum_{Z^{n}\in\mathcal{X}^{n}\times\mathcal{Y}^{n}}\max_{\bm{\theta}\in\bm{\Theta}_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}}P(Z^{n};\,M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},\bm{\theta})$
$\displaystyle=\log\left\\{\sum_{Z^{n}}\prod_{k=0}^{m_{X}-1}\left(\frac{n(X=k)}{n}\right)^{n(X=k)}\prod_{k^{\prime}=0}^{m_{Y}-1}\left(\frac{n(Y=k^{\prime})}{n}\right)^{n(Y=k^{\prime})}\right\\}$
$\displaystyle=\log\left\\{\sum_{X^{n}}\prod_{k=0}^{m_{X}-1}\left(\frac{n(X=k)}{n}\right)^{n(X=k)}\sum_{Y^{n}}\prod_{k^{\prime}=0}^{m_{Y}-1}\left(\frac{n(Y=k^{\prime})}{n}\right)^{n(Y=k^{\prime})}\right\\}$
$\displaystyle=\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n).$

Therefore, the NML codelength of the data $z^{n}$ for
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ is calculated as:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)$
$\displaystyle\phantom{=}+\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y=k^{\prime})\log\frac{n(Y=k^{\prime})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n),$

which is equal to $\mathcal{SC}(z^{n};\mathsf{CAT}^{m_{X}})$.

###  B.2 Confounded Case

We exactly calculate the stochastic complexity defined in Eq. (2) for causal
model $M_{X\leftarrow C\to Y}$:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\leftarrow C\to Y})$
$\displaystyle=\mathcal{SC}(z^{n};M_{X\leftarrow C\to Y})$
$\displaystyle=-\log P(z^{n};M_{X\leftarrow C\to
Y},{\hat{\bm{\theta}}}_{X\leftarrow C\to
Y})+\log\mathcal{C}_{n}(M_{X\leftarrow C\to Y}),$

where ${\hat{\bm{\theta}}}_{X\leftarrow C\to Y}=(\hat{\theta}_{k,k^{\prime}})$
is the maximum likelihood estimator in $\bm{\Theta}_{{X\leftarrow C\to Y}}$.
Each element is
$\hat{\theta}_{k,k^{\prime}}=\frac{n\left(X=k,Y=k^{\prime}\right)}{n}$, where
$n\left(X=k,Y=k^{\prime}\right)$ counts the frequency of data satisfying $X=k$
and $Y=k^{\prime}$ in $z^{n}$. Subsequently, the maximum likelihood as the
first term on the right-hand side is represented as

$\displaystyle\log P(z^{n};M_{X\leftarrow C\to
Y},{\hat{\bm{\theta}}}_{X\leftarrow C\to
Y})=\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(X=k,Y=k^{\prime})\log\frac{n(X=k,Y=k^{\prime})}{n}.$

Since the causal model of $M_{X\leftarrow C\to Y}$ is the model of
$(m_{X}m_{Y})$-categorical distributions, the parametric complexity as the
second term is calculated as

$\displaystyle\log\mathcal{C}_{n}(M_{X\leftarrow C\to
Y})=\log\mathcal{C}_{\mathsf{CAT}}(m_{X}m_{Y},n).$

Thus, the NML codelength of the data $z^{n}$ for $M_{X\leftarrow C\to Y}$ is
calculated as:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\leftarrow C\to Y})$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n\left(X=k,Y=k^{\prime}\right)\log\frac{n\left(X=k,Y=k^{\prime}\right)}{n}+\log\mathcal{C}_{\mathsf{CAT}}(m_{X}m_{Y},n),$

which is equal to
$\mathcal{SC}(z^{n};\mathsf{CAT}^{m_{X}m_{Y}})=\mathcal{SC}(z^{n};\mathsf{HIS}^{m_{X},m_{Y}})$.

###  B.3 Direct Case

First, we consider the first term on the right-hand side of Eq. (8). The first
term, $L(f;M_{X\to Y})$, represents the codelength required to select one
function from a finite set of possible functions $f\in\mathcal{F}$, and the
following holds true:

######  Theorem 1.

Let $\mathcal{F}$ be a set of non-constant functions from $\mathcal{X}$ to
$\mathcal{Y}$. Then, we see

$\displaystyle L(f;M_{X\to Y})=\log|\mathcal{F}|=\log(m_{Y}^{m_{X}-1}-1)$ (21)

###### Proof.

Naively, the number of the possible functions $f$ amounts to $m_{Y}^{m_{X}}$,
but one can remove redundant functions to shorten the resulting codelength.

First, one can remove constant functions since they are associated with the
independence model $M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ and not $M_{X\to
Y}$. $m_{Y}$ such functions exist in total.

Next, distinct functions are associated with the same NML codelength. For a
function $f_{1}$, consider $f_{2}$ given by a constant shift,

$\displaystyle f_{2}(x)=f_{1}(x)+k^{\prime}\pmod{m_{Y}},$

where $k^{\prime}\in\\{1,\cdots,m_{Y}-1\\}$ and $x\in\mathcal{X}$. Now, for
all $k^{\prime}$, the NML codelength with $f_{1}$ is the same as that with
$f_{2}$. Since $m_{Y}$ such different but equivalent functions exist for any
$f_{1}$ including itself, one can further reduce the number of functions by a
factor of $m_{Y}$.

Summing up, there remain
$|\mathcal{F}|=\left(m_{Y}^{m_{X}}-m_{Y}\right)/m_{Y}=m_{Y}^{m_{X}-1}-1$
functions to encode. ∎

As for the second term in Eq. (8), $L(z^{n};M_{X\to Y},f)$, the following
statement holds:

######  Theorem 2.

We define $n(Y=f(X)+k^{\prime})$ as the frequency of data in $z^{n}$ that
satisfies $Y=f(X)+k^{\prime}$. Then, for any $f$, it holds that

$\displaystyle L(z^{n};M_{X\to Y},f)$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y=f(X)+k^{\prime})\log\frac{n(Y=f(X)+k^{\prime})}{n}$
$\displaystyle+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n).$

###### Proof.

Let us denote the probability mass functions of $E_{X}$ and $E_{Y}$ by
$P(E_{X};\bm{\pi}_{X})$ and $P(E_{Y};\bm{\pi}_{Y})$, respectively, where
$\bm{\pi}_{X},\bm{\pi}_{Y}$ are the corresponding parameters.

Now, the observable pair $(X,Y)$ is one-to-one with the exogenous variable
pair $(E_{X},E_{Y})$ when the function $f$ is fixed. Thus, under an
appropriate transformation of data, the joint probability mass function of
$(X,Y)$ is equivalent to that of $(E_{X},E_{Y})$,

$\displaystyle P(X=x,Y=y;M_{X\to Y},f,\bm{\pi}_{X},\bm{\pi}_{Y})$
$\displaystyle=P(E_{X}=x,E_{Y}=y-f(x);M_{X\to Y},f,\bm{\pi}_{X},\bm{\pi}_{Y})$
$\displaystyle=P(E_{X}=x^{\prime};\bm{\pi}_{X})P(E_{Y}=y^{\prime};\bm{\pi}_{Y}),$

which implies the equivalence of $M_{X\to Y}$ with fixed $f$ and
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$

$\displaystyle P(X=x,Y=y;M_{X\to Y},f,\bm{\pi}_{X},\bm{\pi}_{Y})$
$\displaystyle=P(X=x^{\prime},Y=y^{\prime};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},\bm{\pi}_{X},\bm{\pi}_{Y}),$

where $z^{\prime}=(x^{\prime},y^{\prime})=(x,y-f(x))$ is the transformation of
a datum $z=(x,y)$ with fixed $f$.

The equivalence in terms of the probability mass functions immediately extends
to the equivalence in terms of the NML codelengths. Particularly, the NML
codelength of $M_{X\to Y}$ with fixed $f$ is the same as that of
$M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$ with the appropriate
transformation,

$\displaystyle L(z^{n};M_{X\to Y},f)=\mathcal{L}^{d}(z^{\prime
n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}),$

where $z^{\prime n}=(z^{\prime}_{1},...,z^{\prime}_{n})$ and
$z^{\prime}_{i}=(x_{i},y_{i}-f(x_{i}))$ for all $1\leq i\leq n$. This
completes the proof. ∎

Therefore, the codelength of
$z^{n}\in\mathcal{X}^{n}\times\mathcal{Y}^{n}=\\{0,1,\cdots
m_{X}-1\\}^{n}\times\\{0,1,\cdots m_{Y}-1\\}^{n}$ for the causal model
$M_{X\to Y}$ is calculated as:

$\displaystyle\mathcal{L}^{d}(z^{n};M_{X\to Y})$
$\displaystyle=L(z^{n};M_{X\to Y},\hat{f})+L(\hat{f};M_{X\to Y})$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(X=k)\log\frac{n(X=k)}{n}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y=\hat{f}(X)+k^{\prime})\log\frac{n(Y=\hat{f}(X)+k^{\prime})}{n}$
$\displaystyle\phantom{=}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n)+\log(m_{Y}^{m_{X}-1}-1),$

which is equal to
$\mathcal{SC}(x^{n};\mathsf{CAT}^{m_{X}})+\mathcal{SC}((y-\hat{f}(x))^{n};\mathsf{CAT}^{m_{Y}})+L(\hat{f};M_{X\to
Y})$. The subtraction is taken over $\mathbb{Z}/m_{Y}\mathbb{Z}$.

##  Appendix C Proof of Proposition 2

In this section, we provide the proof for Proposition 2.

We first derive the first term on the right-hand side of Eq. (11),
$L(m_{X},m_{Y};M)$, which is the codelength required to encode $(m_{X},m_{Y})$
under causal model $M$. By the Rissanen’s integer coding (Eq. (15)) we have:

$\displaystyle L(m_{X},m_{Y};M)=\log^{*}m_{X}+\log^{*}m_{Y}.$ (22)

For the second term on the right-hand side of Eq. (11),
$L(z^{n};M,m_{X},m_{Y})$, we calculate it under each causal model in the
following sections, and we complete the proof of (13) in Proposition 2 for
each causal relationship:

###  C.1 Independent Case

For given data $z^{n}=(x^{n},y^{n})$, codelength
$L(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},m_{X},m_{Y})$ is calculated
by two NML codes with respect to the histogram models $\mathsf{HIS}^{m_{X}}$
and $\mathsf{HIS}^{m_{Y}}$ for $x^{n}$ and $y^{n}$, respectively. That is, we
have:

$\displaystyle
L(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},m_{X},m_{Y})=\mathcal{SC}(x^{n};\mathsf{HIS}^{m_{X}})+\mathcal{SC}(y^{n};\mathsf{HIS}^{m_{X}}).$

Then, the total codelength
$\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
is expressed as follows:

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=L(m_{X},m_{Y};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})+L(z^{n};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y},m_{X},m_{Y})$
$\displaystyle=\log^{*}m_{X}+\log^{*}m_{Y}+\mathcal{SC}(x^{n};\mathsf{HIS}^{m_{X}})+\mathcal{SC}(y^{n};\mathsf{HIS}^{m_{X}}).$

By the results of Example 1 and 2, the above is written as follows:

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y})$
$\displaystyle=\log^{*}m_{X}+\log^{*}m_{Y}+\mathcal{SC}(x^{n};\mathsf{HIS}^{m_{X}})+\mathcal{SC}(y^{n};\mathsf{HIS}^{m_{X}})$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(X\in I^{X}_{k})\log\frac{n(X\in
I^{X}_{k})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)-n\log
m_{X}+\log^{*}m_{X}$
$\displaystyle\phantom{=}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y\in
I^{Y}_{k^{\prime}})\log\frac{n(Y\in
I^{Y}_{k^{\prime}})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n)-n\log
m_{Y}+\log^{*}m_{Y}$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(\text{disc}(X;m_{X})=k)\log\frac{n(\text{disc}(X;m_{X})=k)}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)+L^{c\to
d}(m_{X},n)$
$\displaystyle\phantom{=}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(\text{disc}(Y;m_{Y})=k^{\prime})\log\frac{n(\text{disc}(Y;m_{Y})=k^{\prime})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n)+L^{c\to
d}(m_{Y},n)$
$\displaystyle=\mathcal{L}^{d}(\text{disc}(x^{n};m_{X}),\text{disc}(y^{n};m_{Y});\text{DISC}(M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y};m_{X},m_{Y}))+L^{c\to
d}(m_{X},n)+L^{c\to d}(m_{Y},n).$

Thus, our claim holds in case $M=M_{X\mathop{\\!\perp\\!\\!\\!\perp\\!}Y}$.

###  C.2 Confounded Case

The codelength $L(z^{n};M_{X\leftarrow C\to Y},m_{X},m_{Y})$ is calculated by
NML code with respect to the histogram model $\mathsf{HIS}^{m_{X},m_{Y}}$ for
$z^{n}$. That is, we have:

$\displaystyle L(z^{n};M_{X\leftarrow C\to Y},m_{X},m_{Y})$
$\displaystyle=\mathcal{SC}(z^{n};\mathsf{HIS}^{m_{X},m_{Y}})$
$\displaystyle=-\max_{p\in\mathsf{HIS}^{m_{X},m_{Y}}}\log
p(z^{n};\bm{\hat{\theta}}(z^{n}))+\log\mathcal{C}_{n}(\mathsf{HIS}^{m_{X},m_{Y}})$

The maximum likelihood estimator for $\mathsf{HIS}^{m_{X},m_{Y}}$ results in
$\hat{\theta}_{k,k^{\prime}}(z^{n})=\frac{n(X\in I^{X}_{k},Y\in
I^{Y}_{k^{\prime}})}{n}m_{X}m_{Y}$. Thus, the maximum log-likelihood of data
is calculated as

$\displaystyle\max_{p\in\mathsf{HIS}^{m_{X},m_{Y}}}\log
p(z^{n};\bm{\hat{\theta}}(z^{n}))$
$\displaystyle=\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(X\in
I^{X}_{k},Y\in I^{Y}_{k^{\prime}})\log\hat{\theta}_{k,k^{\prime}}(z^{n})$
$\displaystyle=\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(X\in
I^{X}_{k},Y\in I^{Y}_{k^{\prime}})\log\frac{n(X\in I^{X}_{k},Y\in
I^{Y}_{k^{\prime}})}{n}+n\log(m_{X}m_{Y}).$

The parametric complexity of $\mathsf{HIS}^{m_{X},m_{Y}}$ is given by

$\displaystyle\log\mathcal{C}_{n}(\mathsf{HIS}^{m_{X},m_{Y}})$
$\displaystyle=\log\int p(z^{n};\bm{\hat{\theta}}(z^{n}))dz^{n}$
$\displaystyle=\log\sum_{Z^{n}\in\\{0,\ldots,m_{X}-1\\}^{n}\times\\{0,\ldots,m_{Y}-1\\}^{n}}\int_{\Delta(Z^{n})}p(z^{n};\bm{\hat{\theta}}(z^{n}))dz^{n}$
$\displaystyle=\log\sum_{Z^{n}\in\\{0,\ldots,m_{X}-1\\}^{n}\times\\{0,\ldots,m_{Y}-1\\}^{n}}\left(\frac{n(X\in
I^{X}_{k},Y\in I^{Y}_{k^{\prime}})}{n}\right)^{n}$
$\displaystyle=\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X}m_{Y},n)$

Consequently, the NML codelength of $z^{n}$ for $\mathsf{HIS}^{m_{X},m_{Y}}$
becomes

$\displaystyle\mathcal{SC}(z^{n};\mathsf{HIS}^{m_{X},m_{Y}})$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(X\in
I^{X}_{k},Y\in I^{Y}_{k^{\prime}})\log\frac{n(X\in I^{X}_{k},Y\in
I^{Y}_{k^{\prime}})}{n}-n\log(m_{X}m_{Y})$
$\displaystyle\phantom{=}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X}m_{Y},n).$
(23)

By comparing the result with Example 1, the total codelength
$\mathcal{L}^{c}(z^{n},m_{X},m_{Y};M_{X\leftarrow C\to Y})$ obtains the
following representation:

$\displaystyle\mathcal{L}^{c}(z^{n},m_{X},m_{Y},M_{X\leftarrow C\to Y})$
$\displaystyle=\log^{*}m_{X}+\log^{*}m_{Y}+\mathcal{SC}(z^{n};\mathsf{HIS}^{m_{X},m_{Y}})$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(X\in
I^{X}_{k},Y\in I^{Y}_{k^{\prime}})\log\frac{n(X\in I^{X}_{k},Y\in
I^{Y}_{k^{\prime}})}{n}-n\log(m_{X}m_{Y})$
$\displaystyle\phantom{=}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X}m_{Y},n)+\log^{*}m_{X}+\log^{*}m_{Y}$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}\sum_{k^{\prime}=0}^{m_{Y}-1}n(\text{disc}(X;m_{X})=k,\text{disc}(Y;m_{Y})=k^{\prime})\log\frac{n(\text{disc}(X;m_{X})=k,\text{disc}(Y;m_{Y})=k^{\prime})}{n}$
$\displaystyle\phantom{=}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X}m_{Y},n)+L^{c\to
d}(m_{X})+L^{c\to d}(m_{Y})$
$\displaystyle=\mathcal{L}^{d}(\text{disc}(x^{n};m_{X}),\text{disc}(y^{n};m_{Y});\text{DISC}(M^{m_{X},m_{Y}}_{X\leftarrow
C\to Y}))+L^{c\to d}(m_{X},n)+L^{c\to d}(m_{Y},n).$

This completes the proof in case $M=M_{X\leftarrow C\to Y}$.

###  C.3 Direct Case

In this case, we also employ two-stage coding for $L(z^{n};M_{X\to
Y},m_{X},m_{Y})$ with respect to function $f$. That is, for a given
$(m_{X},m_{Y})$, we first estimate the optimal function $\hat{f}$ using
maximum likelihood estimation (Algorithm 4), and then encode the data $z^{n}$
based on the histogram model $\mathsf{HIS}^{m_{X},m_{Y}}_{X\to Y}$ with
$\hat{f}$ fixed:

$\displaystyle L(z^{n};M_{X\to Y},m_{X},m_{Y})=L(\hat{f};M_{X\to
Y},m_{X},m_{Y})+L(z^{n};M_{X\to Y},m_{X},m_{Y},\hat{f}),$

where the first term on the right-hand side is given by Eq. (21) and the
second one is calculated by encoding $x^{n}$ and $(y-\hat{f}(x))^{n}$ based on
$\mathsf{HIS}^{m_{X}}$ and $\mathsf{HIS}^{m_{Y}}$, respectively.

Therefore, the codelength is formulated as follows:

$\displaystyle\mathcal{L}^{c}(z^{n};m_{X},m_{Y};M_{X\to Y})$
$\displaystyle=L(m_{X},m_{Y};M_{X\to Y})+L(\hat{f};M_{X\to
Y},m_{X},m_{Y})+L(z^{n};M_{X\to Y},m_{X},m_{Y},\hat{f})$
$\displaystyle=\log^{*}m_{X}+\log^{*}m_{Y}+\log
m_{Y}^{m_{X}-1}+\mathcal{SC}(x^{n};\mathsf{HIS}^{m_{X}})+\mathcal{SC}(y^{n}-\hat{f}(x^{n});\mathsf{HIS}^{m_{Y}}),$

By the results of Example 1 and Example 2, the total codelength in the above
is expressed as follows:

$\displaystyle\mathcal{L}^{c}(z^{n};m_{X},m_{Y};M_{X\to Y})$
$\displaystyle=\log^{*}m_{X}+\log^{*}m_{Y}+\log m_{Y}^{m_{X}-1}$
$\displaystyle\phantom{=}-\sum_{k=0}^{m_{X}-1}n(X\in
I^{X}_{k})\log\frac{n(X\in
I^{X}_{k})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)-n\log m_{X}$
$\displaystyle\phantom{=}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(Y-f(X)\in
I^{Y}_{k^{\prime}})\log\frac{n(Y-f(X)\in
I^{Y}_{k^{\prime}})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n)-n\log m_{Y}$
$\displaystyle=-\sum_{k=0}^{m_{X}-1}n(\text{disc}(X;m_{X})=k)\log\frac{n(\text{disc}(X;m_{X})=k)}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{X},n)$
$\displaystyle\phantom{=}-\sum_{k^{\prime}=0}^{m_{Y}-1}n(\text{disc}(Y-\hat{f}(X);m_{Y})=k^{\prime})\log\frac{n(\text{disc}(Y-\hat{f}(X);m_{Y})=k^{\prime})}{n}+\log\mathcal{C}_{\mathsf{CAT}}(K=m_{Y},n)$
$\displaystyle\phantom{=}+\log m_{Y}^{m_{X}-1}+L^{c\to d}(m_{X},n)+L^{c\to
d}(m_{Y},n)$
$\displaystyle=\mathcal{L}^{d}(\text{disc}(x^{n};m_{X}),\text{disc}(y^{n};m_{Y});\text{DISC}(M_{X\to
Y};m_{X},m_{Y}))+L^{c\to d}(m_{X},n)+L^{c\to d}(m_{Y},n).$

This completes the proof in case $M=M_{X\to Y}$. The same argument holds for
$M=M_{X\leftarrow Y}$

##  Appendix D Description of Tübingen Benchmark Pairs

We provide the detailed descriptions of every dataset [19] we employed in
Section 6.3.1.

#### Discrete Case:

Traffic Dataset (No. 47): This dataset focused on the relationship between the
type of day and traffic volume. $X$ represents the number of cars counted per
24 hours at various stations in Oberschwaben, Germany. $Y$ is categorical,
distinguishing between Sundays plus holidays (labelled as ’1’) and working
days (labelled as ’2’). The ground truth is $X\leftarrow Y$, suggesting that
the type of day influences the traffic volume.

Internet Connections and Traffic Dataset (No. 68): This dataset comes from a
time series study focusing on internet connections and traffic at the MPI for
Intelligent Systems. It features $X$ which represents the bytes sent at minute
and $Y$ which denotes the number of open HTTP connections during that same
minute. Measurements were taken every 20 minutes. The established ground truth
is $Y\text{(open HTTP connections) }$ causes $X\text{(bytes sent)}$.

Direction of Gabor Patches Dataset (No. 107): This dataset originated from a
psychophysics experiment involving human subjects and their perception of
Gabor patches (stripe patterns used in psychological experiments) displayed on
a screen. The Gabor patches were tilted either to the left or right, with
varying contrast levels. $X$ represents the contrast values, ranging from
0.0150 to 0.0500, in increments of 0.0025, and $Y$ is binary, indicating
whether the direction of the tilt was correctly identified or not. $X$ is
regarded as the cause of $Y$.

#### Mixed Case:

Milk Protein Dataset (No. 85): The pair0085 dataset used in our experiments is
a subset of the milk protein trial dataset by Verbyla and Cullis in 1990. The
dataset contains weekly measurements of the assayed protein content of milk
samples taken from 71 cows over a 14-week period. The cows were randomly
allocated to one of three diets: barley, mixed barley-lupins, and lupins. The
variables in the dataset are $X$, representing the time at which the weekly
measurement was taken (ranging from 1 to 14), and $Y$, representing the
protein content of the milk produced by each cow at time $X$. The ground truth
for our experiments was set as $X\to Y$. Note that the dataset does not
consider the effect of the diets on the protein content.

Electricity Consumption Dataset (No. 95): This comprises 9,504 hourly
measurements of total electricity consumption in MWh, denoted as $Y$, in a
region of Turkey. The variable $X$ represents the hour of the day during which
these measurements were taken. The ground truth for this dataset is set as
$X\to Y$, suggesting that the hour of the day is the driving factor for
electricity consumption.

NLSchools Dataset (No. 99): This dataset contains the information of 2287
Dutch eighth graders (about 11 years old) with features $X$(language test
score) and $Y$(social-economic status of pupil’s family). $X\to Y$ is regarded
as the ground truth.

#### Continuous Case:

Cardiac Arrhythmia Database (No. 23): The data, contributed in January 1998,
includes 452 instances with two attributes: age and weight. Age was
hypothesized to influence weight.

Solar Radiation and Air Temperature Dataset (No. 77): This contains daily
measurements of solar radiation in $W/m^{2}$ and the daily average temperature
of the air in Furtwangen, Black Forest, Germany. The dataset covers a time
period from January 1, 1985, to December 31, 2008, with a sample size of
8,401. Solar radiation is denoted by the variable $Y$, while the air
temperature is denoted by the variable $X$. The ground truth for this dataset
was set as $X\leftarrow Y$, indicating that solar radiation is the cause of
air temperature.

Brightness of screen Dataset (No. 101): This is from an experiment that was
performed to generate samples that were clearly unconfounded. $X$ is grey
value of a pixel randomly chosen from a fixed image. The grey value was
displayed by the color of a square on a computer screen. $Y$ is light
intensity seen by a photo diode placed several centimeters away from the
screen. $X$ was the cause of $Y$.