Journal of Machine Learning Research 14 (2013) 111-152

Submitted 10/11; Revised 8/12; Published 1/13

Pairwise Likelihood Ratios for Estimation of
Non-Gaussian Structural Equation Models
Aapo Hyvärinen

AAPO . HYVARINEN @ HELSINKI . FI

Dept of Computer Science and HIIT
Dept of Mathematics and Statistics
University of Helsinki
Helsinki, Finland

Stephen M. Smith

STEVE @ FMRIB . OX . AC . UK

FMRIB (Oxford University Centre for Functional MRI of the Brain)
Nuffield Dept of Clinical Neurosciences
University of Oxford
Oxford, UK

Editor: Peter Spirtes

Abstract
We present new measures of the causal direction, or direction of effect, between two non-Gaussian
random variables. They are based on the likelihood ratio under the linear non-Gaussian acyclic
model (LiNGAM). We also develop simple first-order approximations of the likelihood ratio andanalyze them based on related cumulant-based measures, which can be shown to find the correctcausal directions. We show how to apply these measures to estimate LiNGAM for more than two
variables, and even in the case of more variables than observations. We further extend the method
to cyclic and nonlinear models. The proposed framework is statistically at least as good as existingones in the cases of few data points or noisy data, and it is computationally and conceptually very
simple. Results on simulated fMRI data indicate that the method may be useful in neuroimaging
where the number of time points is typically quite small.
Keywords: structural equation model, Bayesian network, non-Gaussianity, causality, independent
component analysis

1. Introduction
Estimating structural equation models (SEMs), or linear Bayesian networks is a challenging problem with many applications in bioinformatics, neuroinformatics, and econometrics. If the data is
Gaussian, the problem is fundamentally ill-posed. Recently, it has been shown that using the nonGaussianity of the data, such models can be identifiable (Shimizu et al., 2006). This led to the
Linear Non-Gaussian Acyclic Model, or LiNGAM.
The original method for estimating LiNGAM was based on first applying independent component analysis (ICA) to the data and then deducing the network connections from the results of ICA.
However, we believe that it may be possible to develop better methods for estimating LiNGAM
directly, without resorting to ICA algorithms.
c 2013 Aapo Hyvärinen and Stephen M. Smith.

H YV ÄRINEN AND S MITH

A framework called DirectLiNGAM was, in fact, proposed by Shimizu et al. (2011) to provide
an alternative to the ICA-based estimation. DirectLiNGAM was shown to give promising results
especially in the case where the number of observed data points is small compared to the dimensionof the data. It can also have algorithmic advantages because it does not need gradient-based iterative
methods. An essential ingredient in DirectLiNGAM is a measure of the causal direction between
two variables.
An alternative approach to estimating SEMs is to first estimate which variables have connections and then estimate the direction of the connection. While a rigorous justification for such an
approach may be missing, this is intuitively appealing especially in the case where the amount ofdata is limited. Determining the directions of the connections can be performed by considering eachconnection separately, which requires, again, analysis of the causal direction between two variables.
Such an approach was found to work best by Smith et al. (2011) which considered causal analysis
of simulated functional magnetic resonance imaging (fMRI) data, where the number of time points
is typically small. A closely related approach was proposed by Hoyer et al. (2008), in which the PC
algorithm was used to estimate the existence of connections, followed by a scoring of directions by
an approximate likelihood of the LiNGAM model; see also Ramsey et al. (2011).
Thus, we see that measuring pairwise causal directions is a central problem in the theory of
LiNGAM and related models. In fact, analyzing the causal direction between two non-Gaussian
random variables (with no time structure) is an important problem in its own right, and was considered in the literature before the advent of LiNGAM (Dodge and Rousson, 2001).
In this paper, we develop new measures of causal direction between two non-Gaussian randomvariables, and apply them to the estimation of LiNGAM. The approach uses the ratio of the likelihoods of the models corresponding to the two directions of causal influence. A likelihood ratio
is likely to provide a statistically powerful method because of the general optimality properties oflikelihood. We further propose first-order approximations of the likelihood ratio which are easy
to compute and have simple intuitive interpretations. They are also closely related to higher-order
cumulants and may be more resistant to noise. The framework is also simple to extend to cyclic or
nonlinear models.
The paper is structured as follows. The measures of causal directions are derived in Section 2.
In Section 3 we show how to apply them to estimating the model with more than two variables.
The extension to cyclic models is proposed in Section 4 and an extension to a nonlinear model in
Section 5. We report simulations with comparisons to other methods in Section 6, experiments on
simulated brain imaging data in Section 7, and results on a publicly available benchmark data set in
Section 8. Section 9 concludes the paper. Preliminary results were published by Hyvärinen (2010).

2. Finding Causal Direction Between Two Variables
In this section, we present our main contribution: new measures of causal direction between two
random variables.
This section is structured as follows: We first define the problem in Section 2.1. We derive
the likelihood ratio in Section 2.2. We propose a general-purpose approximation for the likelihood
ratio in Section 2.3. The connection to mutual information is explained in Section 2.4. We derivea particularly simple approximation for the likelihood ratio in Section 2.5, and propose an instancefor the case of sparse, symmetric densities. A theoretical analysis of the approximation based on
cumulants is given in Section 2.6. We give intuitive interpretations of the approximations in Sec112

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

tion 2.7, and discuss their noise-tolerance in Section 2.8. Finally, we show how to use the likelihood
ratio approximations in the case of skewed variables in Section 2.9.
For the benefit of the reader, we have further created Table 3 in the Conclusion on page 150 thatlists the main new measures proposed in this paper.
2.1 Problem Definition
Denote the two observed random variables by x and y. Assume they are non-Gaussian, as well as
standardized to zero mean and unit variance. Our goal is to distinguish between two causal models.
The first one we denote by x → y and define asy = ρx + d
where the disturbance d is independent of x, and the regression coefficient is denoted by ρ. The
second model is denoted by y → x and defined as
x = ρy + e
where the disturbance e is independent of y. The parameter ρ is the same in the two models because
it is equal to the correlation coefficient. Note that these models belong to the LiNGAM family(Shimizu et al., 2006) with two variables. In the following, we assume that x, y follow one of these
two models.
Note that in contrast to Dodge and Rousson (2001) or Dodge and Yadegari (2010), we do notassume that d or e are normal, or have zero cumulants. We make no assumptions on their distributions. It is not even necessary to assume that they are non-Gaussian; it is enough that x and y are
non-Gaussian. (This is related to the identifiability theorem in ICA which says that one of the latentvariables can be non-Gaussian, see Comon, 1994).
2.2 Likelihood Ratio
An attractive way of deciding between the two models is to compute their likelihoods and their ratio.
Consider a sample (x1 , y1 ), . . . , (xT , yT ) of data. The likelihood of the LiNGAM in which x → y was
given by Hyvärinen et al. (2010) as#"
yt − ρxt) − T log(1 − ρ2 )log L(x → y) = ∑ Gx (xt ) + Gd ( p21−ρtwhere Gx (u) = log px (u), and Gd is the standardized log-pdf of the residual when regressing y on x.
The last term here is a normalization term due to the use of standardized log-pdf Gd . From this we
can compute the likelihood ratio, which we normalize by T1 for convenience:R=

1
1
log L(x → y) − log L(y → x)TT1
xt − ρytyt − ρxt= ∑ Gx (xt ) + Gd ( p) − Gy (yt ) − Ge ( p). (1)T t
1 − ρ2
1 − ρ2

We can thus compute R and decide based on it what the causal direction is. If R is positive, we
conclude x → y, and if it is negative, we conclude y → x.
113

H YV ÄRINEN AND S MITH

To use (1) in practice, we need to choose the G’s and estimate ρ. The statistically optimal wayof estimating ρ would be to maximize the likelihood, but in practice it may be better to estimate
it simply by the conventional least-squares solution to the linear regression problem. Nevertheless,maximization of likelihood might be more robust against outliers, because log-likelihood functions
often grow more slowly than the squaring function when moving away from the origin.
Choosing the four log-pdf’s Gx , Gy , Gd , Ge could, in principle, be done by modelling the relevant
log-pdf’s by parametric (Karvanen and Koivunen, 2002) or non-parametric (Pham and Garrat, 1997)
methods, which will be discussed in more detail below. However, for small sample sizes suchmodelling can be very difficult. In the following, we provide various parametric approximations.
2.3 Maximum Entropy Approximations of Likelihood Ratio
The likelihood ratio has a simple information-theoretic interpretation which also means we can use
well-known entropy approximations for its practical computation in the case where we do not wantto postulate functional forms for the G’s.
If we take the asymptotic limit of the likelihood ratio, we obtain
ˆ d ) + H(y) + H(ê/σe )R −→ −H(x) − H(d/σ(2)

where we denote differential entropy by H, the estimated residuals by dˆ = y − ρx, ê = x − ρy, andthe variances of the estimated residuals by σ2d , σ2e .
Thus, we can approximate the likelihood ratio using any general, possibly non-parametric, approximations of differential entropy. For example, we can use the maximum entropy approximations by Hyvärinen (1998) which are computationally simple. In fact, we only need to approximate
one-dimensional differential entropies, which is much simpler than approximating two-dimensional
entropies.
One version of the approximations by Hyvärinen (1998) is given byĤ(u) = H(ν) − k1 [E{log cosh u} − γ]2 − k2 [E{u exp(−u2 /2)}]2(3)

where H(ν) = 12 (1 + log 2π) is the entropy of the standardized Gaussian distribution, and the otherconstants are numerically evaluated as
k1 ≈ 79.047,

k2 ≈ 7.4129,

γ ≈ 0.37457.

This approximation is valid for standardized variables; in fact, all the variables in (2) are standardized. The intuitive idea in this approximation is that since the Gaussian distribution has maximum
entropy among all distributions of unit variance, entropy can be approximated by a measure of nonGaussianity which is subtracted from H(ν). The sum of the second and third terms on the right-handside of (3) is a measure of non-Gaussianity (ignoring their negative signs) since the terms are thesquared differences of certain statistics from the corresponding values obtained for a Gaussian distribution. In fact, γ is defined as the expectation of log cosh for a standardized Gaussian distribution,
so the second term on the right-hand side is zero for a Gaussian distribution, just like the skewnesslike statistic measured by the last term.
114

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

The expression in (2) also readily gives a simple intuitive interpretation of the estimation of
causal direction. The (negative) entropies can all be interpreted as measures of non-Gaussianity,
since the variables are standardized. Thus, in (2) we essentially compute the sum of the nonGaussianities of the explaining variable and the resulting residual, and compare them for the two
directions. The directions which leads to maximum non-Gaussianity is chosen.1
2.4 Connection to Mutual Information
It is also possible to give an information-theoretic interpretation which connects the likelihood ratios
to independence measures.
A widely-used independence measure is mutual information, defined for two variables x, y as
I(x, y) = H(x) + H(y) − H(x, y)
where H denotes differential entropy. For a linear transformation
 
 
u
x
=A
,
v
y
we have the entropy transformation formulaH(u, v) = H(x, y) + log | det A|.
On the other hand, the transformation from x, y to x, d has unit determinant, since
  
 
x
1 0
x
=
.
y
a 1
d
Thus, we haveH(x, d) = H(x, y)and likewise for H(y, e). We can now consider the mutual information of the regressors and theresiduals in the two models, and in particular, compute the difference of the mutual informations to
see which one is smaller. In fact, the difference of the mutual informations is asymptotically equal
to the likelihood ratio R sinceI(x, d) − I(y, e) = H(x) + H(d) − H(x, d) − (H(y) + H(e) − H(y, e))ed= H(x) + H(d) − H(y) − H(e) = H(x) + H( ) − H(y) − H( ) − log σd + log σe
σd
σe
d
e
= H(x) + H( ) − H(y) − H( )σdσewhere the joint entropies H(x, e) and H(y, d) as well as the variances of the residuals (which areequal) cancel. Thus, our criterion is equivalent to evaluating the independence of x vs. d and y vs. e
using mutual information, and choosing the direction in which the regressor is more independent of
the residual.
Again, these developments show the important practical advantage that we only need to evaluate one-dimensional entropies, although the definition of mutual information contains a twodimensional entropy.
1. Note that this is not the same as the simple heuristic approach in which we only compute the non-Gaussianities of theactual variables x, y and assume that direction must be from the more non-Gaussian variable to the less non-Gaussian
one.

115

H YV ÄRINEN AND S MITH

2.5 First-Order Approximation of Likelihood Ratio
Next we develop some simple approximations of the likelihood ratio. Our goal is to find causality
measures which are simpler (conceptually and possibly also computationally) than the likelihood
ratio or its general approximation given above.
Let us make a first-order approximationy − ρx) = G(y) − ρ x g(y) + O(ρ2 )
G( p
2
1−ρ

where g is the derivative of G, and likewise for the regression in the other direction. Then, we get
the approximation R̃:R ≈ R̃ =1
ρG(xt ) + G(yt ) − ρxt g(yt ) − G(yt ) − G(xt ) + ρyt g(xt ) = ∑ −xt g(yt ) + g(xt )yt .
T∑
T
t
t

Pham and Garrat (1997) proposed a method for estimating the derivatives of log-pdf’s of randomvariables. Their method could be directly used for estimating g. However, since our main goal here
is to find methods which work for small sample sizes, we try to avoid such estimation of the g’s
which has potentially a very large number of parameters. Instead, here we assume that we havesome prior knowledge on the distributions of the variables in the model. In fact, a result wellknown in the theory of ICA is that it does not matter very much how we choose the log-pdf’s in themodel as long as they are roughly of the right kind (Hyvärinen et al., 2001). This claim is partly
justified by the cumulant-based analysis and simulations below.
In particular, very good empirical results are usually obtained by modelling any sparse (i.e.,
super-Gaussian, or positively kurtotic), symmetric densities by either the logistic densityπG(u) = −2 log cosh( √ u) + const.
2 3
or the Laplacian density(4)√
G(u) = − 2|u| + const.

where the additive constants are immaterial. The Laplacian density is not very often used in ICA
because its derivative is discontinuous at zero which leads to problems in maximization of the ICA
likelihood. However, here we do not have such a problem so we can use the Laplacian density as
well.
Thus, if we approximate all the log-pdf’s by (4), we get the “non-linear correlation”
R̃sparse = ρÊ{x tanh(y) − tanh(x) y}(5)πwhere we have omitted the constant 2√
which is close to one, as well as a multiplicative scaling
3

constant. Here, Ê means the sample average. This is the quantity we would use to determine thecausal direction. Under x → y, this is positive, and under y → x, it is negative.
2.6 Cumulant-Based Approach
To get further insight into the likelihood ratio approximation in (5), we consider a cumulant-based
approach which can be analyzed exactly. The theory of ICA has shown that cumulant-based approaches can shed light into the convergence properties of likelihood-based approaches. However,116

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

cumulant-based methods tend to be very sensitive to outliers, so their utility is mainly in the theoretical analysis; for analysing real data, the measure in (5) is preferred.
Here, an approach based on fourth-order cumulants is possible by definingR̃c4 (x, y) = ρÊ{x3 y − xy3 }(6)where the idea is that the third-order monomial analyzes the main nonlinearity in the nonlinear
correlation. In fact, we can approximate tanh by a Taylor expansion1tanh(u) = u − u3 + O(u5 ).
3(7)Then, first-order terms are immaterial because they produce terms like Ê{xy − xy} which cancel out,
and the third-order terms can be assumed to determine the qualitative behaviour of the nonlinearity.
Our main results of the cumulant-based approach is the following:
Theorem 1 If the causal direction is x → y, we haveR̃c4 = kurt(x)(ρ2 − ρ4 )(8)

where kurt(x) = E{x4 } − 3 is the kurtosis of x. If the causal direction is the opposite, we haveR̃c4 = kurt(y)(ρ4 − ρ2 ).

(9)

Proof Consider the fourth-order cumulant
C(x, y) = cum(x, x, x, y) = E{x3 y} − 3E{xy}where we assume the two variables are standardized. We have kurt(x) = C(x, x) = cum(x, x, x, x).
The nonlinear correlation can be expressed using this cumulant asR̃c4 = ρ[C(x, y) −C(y, x)]since the linear correlation terms cancel out. We use next two well-known properties of cumulants.
First, the linearity property says that for any two random variables v, w and constants a, b we havecum(v, v, v, av + bw) = a cum(v, v, v, v) + b cum(v, v, v, w)
and second, cum(v, w, x, y) = 0 if any of the variables v, w, x, y is statistically independent of theothers. Thus, assuming the causal direction is x → y, that is, y = ρx + d with x and d independent,
we haveR̃c4 = ρ[cum(x, x, x, ρx + d) − cum(x, ρx + d, ρx + d, ρx + d)]= ρ[ρcum(x, x, x, x) + cum(x, x, x, d)3− ρ cum(x, x, x, x) − 3ρ2 cum(x, x, x, d) − 3ρcum(x, x, d, d) − cum(x, d, d, d)= ρ[ρkurt(x) − ρ3 kurt(x)] = kurt(x)[ρ2 − ρ4 ]which proves (8). The proof of (9) is completely symmetric: exchanging the roles of x and y will
simply change the sign of the nonlinear correlation, and the kurtosis will be taken of y.
117

H YV ÄRINEN AND S MITH

The regression coefficient ρ is always smaller than one in absolute value, and thus ρ2 − ρ4 > 0.
Assuming that the relevant kurtosis is positive, which is very often the case for real data, the signof R̃c4 can be used to determine the causal direction in the same way as in the case of the likelihood
approximation R̃ in (5). Thus, the cumulant-based approach allowed us to prove rigorously that anonlinear correlation of the form (6) can be used to infer the causal direction, since it takes opposite
signs under the two models. Note that this nonlinear correlation has exactly the same algebraic
form as the likelihood ratio approximation (5); only the nonlinear scalar function is different. In
particular, this analysis shows that the exact form of the nonlinearity is not important: the cubic
nonlinearity is valid for all distributions of positive kurtosis.
If the relevant kurtosis is negative, a simple change of sign is needed. In general, we should thus
multiply R̃c4 by the sign of the kurtosis to obtain
R̃′c4 (x, y) = sign(kurt(x))ρÊ{x3 y − xy3 }.
Here, we get the complication that we have to choose whether we use the sign of the kurtosis of x
or y. Usually, however, the signs would be the same, and we might have prior information on theirsign, which is in most applications positive.2
Related cumulant-based measure were proposed by Dodge and Rousson (2001) and Dodge andYadegari (2010). Their fourth-order measures used the ratio of marginal kurtoses, as opposed tothe cross-cumulants we use here. They further assumed the disturbances to be Gaussian (or at least
to have zero cumulants), which makes their measures less general than ours. In fact, their method
relies on the fact that kurtosis is decreased by adding a Gaussian disturbance, but if the disturbanceis much more kurtotic than the regressor, the opposite can be the case.
2.7 Intuitive Interpretations
Next, we provide some intuitive interpretations of the obtained first-order approximations of thelikelihood ratio.
2.7.1 G RAPHICAL I NTERPRETATIONThe cumulants and nonlinear correlations have a simple intuitive intepretation. Let us consider the
cumulant first. The expectations E{x3 y} or E{xy3 } are basically measuring points where both x
and y have large values, but in contrast to ordinary correlation, they are strongly emphasizing large
values of the variable which is raised to the third power.
Assume the data follows x → y, and that both variables are sparse (super-Gaussian). Then, bothvariables simultaneously have large values mainly in the cases where x takes a large value, making
y large as well. Now, due to regression towards the mean, that is, |ρ| < 1, the value of x is typicallylarger than the value of y. Thus, E{x3 y} > E{xy3 }. This is why E{x3 y} − E{xy3 } > 0 under x → y.
The idea is illustrated in Figure 1.
2. In the general case where the (real) kurtoses of x and y are allowed to have different signs, we need to compute
two quantities: R̃′c4 (x, y) = sign(kurt(x))ρÊ{x3 y − xy3 } and R̃′c4 (y, x) = sign(kurt(y))ρÊ{y3 x − yx3 }. According to
the analysis above, the former quantity is positive if x → y, and the latter quantity is positive if y → x. However, inpractice, this does not lead to a simple decision rule since due to finite sample size, or violations of the model, it
could be that both of these quantities are positive, or none of them. In such cases, the decision rule should be definedso as to indicate that the causal direction could not be decided.

118

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM Syx

Figure 1: Intuitive illustration of the nonlinear correlations. Here, x → y and the variables are very
sparse. The nonlinear correlation E{x3 y} is larger than E{xy3 } because when both variables are simultaneously large (the “arm” of the distribution on the right and the left), x
attains larger values than y due to regression towards the mean.

This interpretation is valid for the tanh-based nonlinear correlation as well, because we can use
the function h(u) = u − tanh(u) instead of tanh to measure the same correlations but with opposite
sign. In fact, we haveR̃sparse = ρÊ{h(x)y − x h(y)}because the linear terms cancel each other. The function h is a soft thresholding function, and thus
has the same effect of emphasizing large values as the third power. Thus the same logic applies forh and the third power.
2.7.2 I NTERPRETATION AS I MPLICATION
Even if the data is not assumed to follow any particular model, the nonlinear correlation could beinterpreted as a logical implication. In general, if the existence of event A implies the existence of
event B, but there is no implication in the other direction, a causal influence from A to B might beinferred. Since A ⇒ B is equivalent to ¬B ⇒ ¬A, there has to be some clear distinction between theevents and their negations for this interpretation to be meaningful. We assume here that the events
are rare, that is, have small probabilities.
Now, let us consider the events A, defined as “x takes a very large value” and B, defined as “ytakes a relatively large value of the same sign as x”. Notice that because the regression coefficientis smaller than one, we cannot require y to take particularly large values. It is assumed here that thethresholds for deciding when a value is large are chosen so that both of these events are rare.
To investigate implication, we can consider how to refute it. To refute A ⇒ B, we should consider
cases where x takes a very large value but y takes a value of the opposite sign. This can be measured
by Ex3 (−y) where x3 looks at large values of x and the minus sign changes this into a measure of
how much large values of x coexist with y’s of opposite sign.
119

H YV ÄRINEN AND S MITH

Thus, Ex3 y − Exy3 can be seen as measuring of how much evidence we have to refute y ⇒ x(latter term) minus the evidence to refute x ⇒ y (negative of first term). If it is large, we accept the
implication x ⇒ y together with its causal interpretation.
It might be argued that the connection between causality and implication could also plausibly
be defined in the opposite direction: If A implies B as defined above, then B causes A. However,
we shall now argue that the interpretation we gave above follows naturally from the definition of a
SEM with two variables. Assume x → y and ρ > 0. If x is very large, y is likely to be large and of
the same sign, since it is not very probable that d would cancel out the effect of ax. Thus, we haveA ⇒ B when x causes y under the SEM framework.
2.8 Noise-Tolerance of the Nonlinear Correlations
An interesting point to note is that the cumulant in (6) is, in principle, immune to additive measurement noise. Assume that instead of the real x, y, we observe noisy versions x̃ = x + n1 and ỹ = y + n2
where the noise variables are independent of each other and x and y. By the basic properties ofcumulants (see proof of Theorem 1), the nonlinear correlations are not affected by the noise at all in
the limit of infinite sample size. Thus, our method in not biased by noise. This is in stark contrast to
ICA algorithms which are strongly affected by additive noise; thus ICA-based LiNGAM (Shimizu
et al., 2006) would not yield consistent estimators in the presence of noise.
To be more precise, we haveE{x̃3 ỹ} − E{x̃ỹ3 } = cum(x̃, x̃, x̃, ỹ) − cum(x̃, ỹ, ỹ, ỹ)= cum(x, x, x, y) − cum(x, y, y, y) = E{x3 y} − E{xy3 }due to the independence of n1 and n2 of the other variables and each other.
On the other hand, the estimation of ρ is strongly affected by the noise. This implies that R̃c4 is
not immune to noise. Nevertheless, measurement noise would only decrease the absolute value ofρ and not change its sign. Thus, the sign of R̃c4 is not affected by additive measurement noise in thelimit of infinite sample. This applies for both Gaussian and non-Gaussian noise.
The fact that the ρ is only a multiplicative scaling in the nonlinear correlations (6) or (5) must
be contrasted with its role in the likelihood ratio (1) where its effect is more complicated. Thus,
when ρ is underestimated due to measurement noise, it may have a stronger effect on the likelihood
ratio, while its effect on the nonlinear correlations is likely to be weaker. While this logic is quite
approximative, simulations below seem to support it.
On the other hand, the standardization of the variables is also affected by noise, in particular if
the noise variances are not equal. As long as the noise variances are equal, the error in standardization will affect the measures by a multiplicative constant only, effectively making the cumulants
smaller. Thus, the noise-tolerance of the cumulants may be useful in practice only if the variancesof the noise variables are equal.
2.9 Skewed Variables
Above, the likelihood ratio approximations and cumulants were developed for sparse, typicallysymmetrically-distributed variables. Here, we consider the extension to skewed variables. Again,
the underlying motivations is that if we know the distributions are skewed, we can use this priorknowledge to obtain particularly simple measures of causal direction. The cumulant-based analysis
120

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

is mainly for theoretical interest due to the sensitivity of cumulants to outliers; we provide a morerobust nonlinearity for analysing real data.
2.9.1 C UMULANT-BASED A PPROACH
The cumulant-based approach allows for a very simple extension of the framework to skewed variables. As a simple analogue to (6), we can define a third-order cumulant-based statistic as follows
R̃c3 (x, y) = ρÊ{x2 y − xy2 }.

(10)

The justification for this definition is in the following theorem, which is the analogue of Theorem 1:
Theorem 2 If the causal direction is x → y, we haveR̃c3 = skew(x)(ρ2 − ρ3 )(11)

and if the causal direction is the opposite, we haveR̃c3 = skew(y)(ρ3 − ρ2 ).

(12)

Proof Consider the third-order cumulant
C(x, y) = cum(x, x, y) = Ex2 ywhere we assume the two variables are standardized. We have skew(x) = C(x, x) = cum(x, x, x).
The nonlinear correlation can be expressed using this cumulant asR̃c3 = ρ[C(x, y) −C(y, x)].
Assuming the causal direction is x → y, we haveR̃c3 = ρ[cum(x, x, ρx + d) − cum(x, ρx + d, ρx + d)]= ρ[ρ cum(x, x, x) + cum(x, x, d) − ρ2 cum(x, x, x) − 2ρ cum(x, x, d) − cum(x, d, d)]= ρ[ρ skew(x) − ρ2 skew(x)] = skew(x)[ρ2 − ρ3 ]which proves (11). The proof of (12) is again completely symmetric.
To use the measure (10) in practice, we have to take into account the fact that we cannot assume,in general, the skewnesses of the variables to have some particular sign. In some applications this ispossible: For example, in resting-state fMRI data it might be safe to assume that the skewnesses are
all positive because it is much more common that the signals obtain large values due to activation
than due to inhibition (however, this point needs to be confirmed by empirical investigations offMRI data).
In the general case, we propose that before computing these nonlinear correlations, the signsof the variables are first chosen so that the skewnesses are all positive. This can be accomplishedsimply by multiplying the variables by the signs of their skewnesses to get a new variable x∗x∗ = sign(skew(x)) x
121

(13)

H YV ÄRINEN AND S MITH

and the same for y (this transformation has to be done before computing ρ). Now, we have a situation
similar to the previous measures: Under x → y, R̃′c3 (x, y) > 0. This is because again, |ρ| < 1, andtherefore ρ2 − ρ3 > 0 regardless of the sign of the coefficient. Likewise, for y → x, R̃′c3 (y, x) < 0.
Our measure is related to the directionality measure proposed by Dodge and Rousson (2001),
which in our notation would be:
R̃DR (x, y) = [Ê{x2 y}]2 − [E{xy2 }]2(14)

which has the advantage of of being particularly simple, and does not require the skewnesses to beof any particular sign. However, our measure is more closely related to likelihood ratios which may
give it some advantage in terms of statistical performance, as will be seen in the simulations below.
2.9.2 ROBUST, L IKELIHOOD -BASED A PPROACH
The skewed case might also be approached by defining a skewed log-pdf and using the methods
in previous sections. Unfortunately, in the theory of ICA, general-purpose skewed densities canhardly be found, and thus it is not clear how to define such densities and how generally they would
be applicable. Nevertheless, a likelihood-based approach is likely to be more robust against outliers
than the cumulant-based one (unless the model pdf has very light tails) which is why we develop
one here.
We propose the following nonlinearity:
gskew (x) = log cosh(max(x, 0))(15)which can be justified as follows. Consider the following family of pdf’s, defined using the derivative of the log-pdf(log p)′ (x) = gskew (x) − βx − α(16)where β and α are parameters. Let us take α and β so that we get a standardized pdf with zeromean and unit variance. Numerical calculations show that this is obtained by values which areapproximately α0 = 0.188 and β0 = 1.32. The ensuing pdf is illustrated in Figure 2.
Further numerical calculations show that the higher-order cumulants of the standardized pdf areboth positive: Skewness is approximately 0.37 and kurtosis 0.47.
Now, we can add any linear function and/or constant to (log p)′ without changing the value ofthe approximative likelihood ratio in (5). In particular, using the true derivative of log-pdf in (16) isequivalent to using the algebraically simpler gskew .
Thus, we obtain the following approximation for the likelihood ratio:R̃skrb (x, y) = ρÊ{gskew (x)y − xgskew (y)}(17)

with gskew defined in (15). Again, this applies for positively skewed variables only. If the skewnesses
are not known a priori, they can be made positive by (13).

3. Estimating a Network with More Than Two Variables
In this section, we consider the general case of more than two variables. We present two approaches:First, we use the pairwise analysis in a DirectLiNGAM framework, and second, we present a twostage method where the possible connections in a sparse graph are first pruned using covarianceinformation.
122

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

0.5

0

0.4
−5

0.3
−10

0.2
−15

0.1

0
−5

0

5

−20
−5

0

5

Figure 2: The pdf for robust modelling of skewed densities. Left: the pdf corresponding to thederivative of log-pdf in (16) is plotted (solid curve) with α and β chosen so that thedensity is standardized. For comparison, the Gaussian density of the same mean andvariance is plotted as well (dashed). Right: the logarithms of the same density functions.

3.1 Model Definition
Denote by x = (x1 , x2 , . . . , xn )T the vector of observed variables. The linear non-Gaussian acyclic
model (LiNGAM) proposed by Shimizu et al. (2006) can be expressed as
x = Bx + e
where e is the vector of disturbances, and B is the matrix that describes the influences of the xi oneach other; the diagonal of B is defined to be zero.
It was shown by Shimizu et al. (2006) that the model is identifiable under the following assumptions: a) the ei are non-Gaussian, b) the ei are mutually independent, and c) the matrix B
corresponds to a directed acyclic graph (DAG). It is well-known that the DAG property is equivalent to an existence of an ordering of the variables xi (not necessarily unique) in which there are
only connections “forward” in the ordering; if the variables are re-ordered according to the causal
ordering, the matrix B has all zeros above the diagonal.
3.2 Using Pairwise Measures in the DirectLiNGAM Framework
The first way to use the pairwise analysis developed above to estimate LiNGAM which has morethan two variables is to use the DirectLiNGAM framework (Shimizu et al., 2011).
3.2.1 F INDING ROOT OF G RAPH
In the DirectLiNGAM approach, we first compute the likelihood ratios of all different pairs ofvariables, and store the log-likelihood ratio for xi and x j as the (i, j)-th entry of a matrix M. Alternatively, we can use the likelihood ratio approximations which can be all subsumed under the
algebraic form
M = C ⊙ E{xg(x)T − g(x)xT }
(18)
123

H YV ÄRINEN AND S MITH

where ⊙ is element-wise multiplication. The nonlinearity g is typically chosen so that it is g(u) =tanh(u) for symmetric sparse data and g(u) = −u2 or the function in (15) for skewed data. C isthe covariance matrix of the data; since the data is assumed standardized C equals the matrix ofcorrelation coefficients.
Now, for the variables xi which have no parents, all entries in the i-th row of M are non-negative,
neglecting random errors. (Note that there is no reason why there would be only one such “root”
variable.) This was shown to be exactly true for the cumulant-based approaches g(u) = −u3 andg(u) = −u2 (assuming that the kurtoses or skewnesses, respectively, are positive) and is true as afirst-order approximation based on (7) for g(u) = tanh(u). The reverse also holds if we assume
faithfulness.3
Thus, we first find the row, say with index i∗ , which is most likely to have all non-negative entries(the actual estimation procedure is considered below). Then, we regress (“deflate”) the variable xi∗
out of all the other variables (Shimizu et al., 2011). We iterate this procedure by computing M
again for the deflated x. By locating the row which is most likely to have only non-negative entries
in the newly computed M, we thus find a variable which has no parents except for possibly the firstvariable found in the previous step. Repeating this, we find variables which are next in the partial
order given by the DAG. Thus in the end we have the causal ordering of the variables.
After such estimation of the causal ordering, estimating the coefficients bi j is easy by just ordinary least-squares estimation (Shimizu et al., 2006).
Alternatively, we could use a simple approximation which is very simple and computationallyefficient. Instead of carrying out deflation by regression as described above, we simply remove theentries of the rows and columns corresponding to the already “found” variables in the matrix M, anditerate the procedure. Thus, we obtain the causal ordering directly from a single matrix of nonlinear
correlations, without any deflation. This is an approximation with no rigorous justification (because
when removing the root we should also remove its effect on all the entries of M) and it is likely
to be inconsistent. However, in simulations reported below it works quite well. It has the benefit
of being computationally extremely simple, and it gives a simple conceptual link between causal
ordering and the nonlinear correlations and cumulants.
3.2.2 AGGREGATING PAIRWISE M EASURES
To use the method just described we have to solve the problem of aggregating the pairwise measures.
We need to find the row which is most likely to be all non-negative up to random errors. Obviously,
we could just take the sums of the entries in each row and locate the maximum sum but this is not
likely to be optimal. So, we next develop a more principled way of aggregation.
Consider the mi j , j = 1, . . . , n for a fixed i, which are the estimates of pairwise likelihood ratios
or some approximations. Assume they are independent and have Gaussian distributions N(µi j , σ2 ),where the variances are assumed to be equal for simplicity. The variance σ2 is the estimation error
due to finite sample, and the µi j are the true values. The posterior of µi j given mi j is then Gaussian
3. For a variable x0 with no parents, any other variable is of the form x j = ax0 + d where a expresses the total effectcoming from x0 , and d is a sum of the inputs from other external influences, which are, by definition, independentof x0 . Thus, the pairwise model holds with a d independent of x0 and the pairwise measure is non-negative. On theother hand, consider xi which does have parents. Now, go backwards in the graph until you find a node x0 which has
no parents (in a DAG, such a variable is guaranteed to exist). According to the logic just given, we have xi = ax0 + d,again with an independent d. By faithfulness, a 6= 0. Since changing the direction simply changes the sign of ourmeasures, there will be a negative entry in the i-th row, and it has to be non-zero.

124

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

with mean mi j and variance σ. Thus, the posterior log-probability that all of the µi j , j = 1, . . . , n are
positive can be calculated as
mi j
mi j
µi j − mi j>−
|mi j ) = ∑ log Φ()(19)
log ∏ P(µi j > 0|mi j ) = log ∏ P(
σ
σ
σ
j
j
j
where Φ is the cumulative distribution function of the standardized Gaussian distribution. Estimating σ is possible but we prefer to assume it is very small and make the following approximation:mi j
1) ≈ − 2 min(0, mi j )2σ
2σwhich can be seen to be quite accurate by a simple numerical comparison, and avoids numericalproblems in computing the logarithm of Φ for large negative values. Now, σ is simply a multiplicative scaling constant which can be ignored when comparing estimates of the log-probabilities in(19).
Thus, we propose the following way of aggregating the pairwise likelihood ratios. Compute foreach row of M
mi = − ∑ min(0, [M]i j )2
log Φ(

jwhich, intuitively speaking, punishes violations of the positivity. The index i∗ with maximum mi isthus taken as the estimate of a variable with no parents, that is, a first variable in the causal ordering.
3.3 Two-Stage Approach to Estimating a Sparse Model
If the matrix B is known to be sparse, we can use a two-stage method in which we first estimate theconnections in an undirected sense, and then find their directions using our pairwise method. This
two-stage method is interesting from the viewpoint of clearly dividing the estimating problem into
two parts.
We first find undirected connections by using any known method for estimating a Gaussian
undirected model (Spirtes et al., 1993). In the simplest case, this can be based on the inverse
covariance matrix, or the precision matrix. As is well-known in the theory of Gaussian graphicalmodels, there is an intimate connection between the non-zero entries in the precision matrix andthe existence of connections in the SEM—although the connection is not quite simple, especiallyfor directed graphs. In contrast, the direction of a connection cannot be easily determined from thecovariances, and is often unidentifiable, which was of course the original motivation for introducing
non-Gaussian models (Shimizu et al., 2006). Nevertheless, as a first approximation, we can prune
the set of candidate connections using the inverse covariance matrix, and apply our pairwise analysis
only on those connections which this covariance-based analysis indicates to be present.
In an estimated inverse covariance matrix, there are of course no exact zeros. Thus, we usebootstrapping to test if each entry is non-zero. That is, we draw bootstrap samples of the data,and compute the inverse covariance for each such sample. The ratio of the mean and the standard
deviation of the bootstrap estimates of any given entry is then compared with the relevant quantile
of a standardized Gaussian distribution.4 The test is made separately for each non-diagonal entry of
the inverse covariance matrix.
4. In the simulations below, we also tried methods for sparse estimation of the inverse covariance matrix. However, wefound that this simple testing procedure works by far the best. The sparse inverse estimation methods are, in fact,
developed for the case of a very large number of variables, and thus may not be useful in our simulations where wetypically have 5-10 variables only.

125

H YV ÄRINEN AND S MITH

Depending on the goal of the analysis, it may or may not be necessary to do corrections for
multiple testing. If we do such corrections, we can actually claim that the connections found arestatistically significant. However, this is obtained at the cost of a large number of false negatives. On
the other hand, if we simply consider the existence of the connections as another set of parameters
to estimate, it may be more advantageous not to make such corrections to reduce the overall error
rate. In fact, a false negative (setting an existing connection to zero) could be considered quite aserious error in this context, so we prefer to use a rather large α. In the simulations below, we thus
set the false positive rate α = 0.01 with no correction for multiple testing. Such corrections will ofcourse be needed if our aim were to claim that a particular connection exists, but if our goal hereis mainly the inference of the causal ordering, some false positives should not matter since they arelikely to correspond to small values of the coefficients anyway.
Then, for each of those significantly non-zero connections, we determine the direction of causality using our pairwise tests. There is no need to do any kind of deflation anymore. If we want toconvert the obtained estimates into a total ordering of the variables, we input those connections
which were not pruned to the ordering method presented by Shimizu et al. (2006).

4. Estimating Cyclic Models
An important generalization of the DAG framework would be to estimate cyclic models. Here, we
assume the following well-known generative model for the data. First, the external influences arrivein the system at time t = 0x(0) = ewhere x(t) is value a hypothetical dynamic system at time point t. Then, at subsequent time steps,
the external influence is completemented by feedback asx(t + 1) = e + Bx(t)where the matrix B has zero diagonal, which means we do not allow self-loops. Assuming that B
is stable in the sense that its largest eigenvalue is smaller than one in absolute value, we have in thelimit
x = ∑ Bk e = (I − B)−1 ek≥0

and thusx = Bx + e(20)

where B is now allowed to be cyclic. This gives a simple interpretation of a model of the form(20) in the case where B is allowed to be cyclic. As above, the ei are assumed independent andnon-Gaussian.
In fact, estimation of such a model by ICA is possible if B is small enough, namely if all its
entries are smaller than one in absolute value. Then, it is possible to estimate the model even by
ICA, since after estimating ICA, we can find the right permutation of the components based on
putting the largest entries of each row in the diagonal. Thus, the model is identifiable under theseassumptions. This is shown in detail in the following Theorem:
Theorem 3 Assume that the data follows the cyclic LiNGAM model in (21) with no self-loops.
Assume further that all the entries in the matrix B have absolute values smaller than one, and that
126

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

the dominant eigenvalue of B is smaller than one in absolute value. Then, the model is uniquely
identifiable, that is, the matrix B can be estimated from the data without any ambiguity.
Proof: The data actually follows the ICA model as
x = (I − B)−1 e.

(21)

The ICA model is known to be identifiable up to a) the ordering of the components and b) a scalarmultiplier for each of the components (Comon, 1994). The unidentifiability of the scalar multiplier
disappears here because by definition, the diagonal of the inverse of the mixing matrix has all ones
due to the diagonal matrix in (21). Thus, it was shown by Shimizu et al. (2006) that this implies
the idenitifiability of the LiNGAM model if we can solve the indeterminacy of the permutation.
Acyclicity was used for this purpose by Shimizu et al. (2006). Here, we use the assumption ofabsolute values smaller than one. In fact, consider the estimate of the inverse of the mixing matrix.
Normalize it by dividing each row by its maximum element. Then, it equals I − B up to a random
permutation of the rows. Due to our assumption of B, all non-diagonal entries in this matrix aresmaller than one in absolute value. Thus, the original (correct) permutation of the rows can befound by locating on each row i the unique entry which is equal to one. Denoting its column indexby j(i), the original matrix is given by permuting the rows of the matrix to the ordering given byj(i), that is, the ordering which puts the ones in the diagonal. This also suggests that we can estimate the model using the sparse graphs idea above. We
prune the inverse covariance matrix to find where there are (probably) connections, and then findthe directions of the connections using our pairwise measures. Using pairwise connections makes
sense if we further assume that there are no pairwise loops, that is, connections xi → x j and x j → xi
are not both non-zero. The main justification for this approach is that since the connections are
weak, one can assume that the cyclicity has little effect on local pairwise measures. However, an
exact convergence of such a method to the right parameter values does not seem possible to show
in general.

5. Estimation in Case of Nonlinear Relations
In this Section, we generalize our method to a nonlinear model.
5.1 Definition of Nonlinear Model
Another interesting extension of the linear causal models is obtained by considering nonlinearities
instead of non-Gaussianities (Hoyer et al., 2009). We define the two models as follows. The firstone, x → y, is given by
y = f (x) + d
where f is a nonlinear function, not necessarily invertible or even differentiable. The disturbanced is again independent of x. Both x and y are standardized to unit variance. The second model isdenoted by y → x and defined as
x = g(y) + e
where g is another nonlinear function, and e is a disturbance variable. Other approaches to inferring the causal direction with nonlinear relations were introduced by Zhang and Hyvärinen (2009),
Daniušis et al. (2010) and Mooij et al. (2010).
127

H YV ÄRINEN AND S MITH

5.2 Likelihood Ratio for Nonlinear Model
The likelihood of the model x → y can be obtained as the sum of the log-prior of the variable x andthe log-likelihood of the residual:
log p(x, y) = log px (x) + log pd (y − f (x)) = Gx (x) + Gd (dˆ) − log σdσdwhere we denote, like above, the variance of the standardized residual by σ2d , the log-pdf of thestandardized residual by Gd , and the log-pdf of x by Gx . Thus, like in the linear case, we obtainR=

xt − g(yt )yt − f (xt )1) − Gy (yt ) − Ge ()
Gx (xt ) + Gd (
T∑
σ
σe
d
t



− log σd + log σe . (22)An important difference to the linear case is that the variances of the residuals need not be equal,
σd 6= σe , so they do not cancel. In an information-theoretic formulation, we obtain asymptoticallyˆ d ) + H(y) + H(ê/σe ) − log σd + log σe .
R −→ −H(x) − H(d/σ(23)

We can approximate R using the same maximum entropy approximations (Hyvärinen, 1998) asin the linear case in Section 2.3. The only difference is that we need to add the log-variances of theresiduals to the expression. Thus, an important advantage of our approach is that we do not need any
measures of independence per se; estimation of one-dimensional differential entropies is sufficient.
On the other hand, it may be advantageous to adapt the approximation to the nonlinear case.
First, it does not seem useful to consider the prior non-Gaussianities of the variables, since a nonlinear mixing can change non-Gaussianities in completely unpredictable ways. This is unlike in thecase of ICA, where a linear mixing decreases non-Gaussianity. Second, we can assume that theresiduals tend to be sparse, and model them as Laplacian. This has the further advantage of making
the measure more robust to outliers.
Now, for a Laplacian variable, the scale parameter σ is most naturally estimated as the mean
absolute deviation (MAD), which is the maximum likelihood estimate. If we plug this estimate in
the likelihood ratio, and omit the priors on x and y, we have#√√1
2|yt − f (xt )|2|xt − g(yt )|R=− log σ̂d + log σ̂e
)+
−
σ̂d
σ̂e
T∑
t
√ σ̂d √ σ̂e= − 2 + 2 − log σ̂d + log σ̂e
σ̂d
σ̂e
"

which gives finally the following objective
ˆ + log Ê{|ê|}R̃mad = − log Ê{|d|}(24)where Ê denotes the sample average, and thus Ê{|.|} denotes the MAD. In other words, we have an
objective which simply compares the mean absolute deviations in the two cases.
The likelihood ratio depends on the estimated nonlinearities f , g. The estimation of f and g canbe done with classic least-squares estimation methods independently of any developments in thispaper. A large number of non-parametric methods have been developed in the literature, see, forexample, Hoyer et al. (2009) for an example.
128

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

5.3 Connection to Independence-Based Nonlinear Methods
In fact, our method has a close connection to the independence-based method by Hoyer et al. (2009),
generalizing the connection shown in Section 2.4. Using basic information-theoretic properties, wehave under x → yH(x, y) = H(x) + H(y|x) = H(x) + H(y − f (x)|x) = H(x) + H(d|x) = H(x, d)and likewise, this is equal to for H(y, e). Now, just like in the linear case, we can consider thedifference between the mutual informations of the regressors and residuals in the two directions,and obtain
I(x, d) − I(y, e) = H(x) + H(d) − H(y) − H(e) = H(x) + H(de
) − H(y) − H( ) + log σd − log σeσdσewhere two terms equal to h(x, y) cancel. Here, we see that asymptotically, our objective derived fromthe likelihood ratio is equal to the difference of the two mutual informations (with sign reversed).
Its sign tells which mutual information is larger, and in particular, in which direction the residual of
the regression is more independent. Thus, using the likelihood ratio is equivalent to using mutual
information as independence measure in the methods by Hoyer et al. (2009).
The developments given above thus show that when comparing independencies of the residuals
like Hoyer et al. (2009), it is not necessary to explicitly estimate mutual information; estimation of
one-dimensional entropies leads to an equivalent result.

6. Simulations
We conducted simulations comparing the different methods proposed in this paper, as well as previously proposed LiNGAM estimation methods. In all the simulations, we emphasize difficult conditions. In most of the simulations, this means the case where the number of observations is small; theexception being the simulations with added measurement noise. We also take weakly non-Gaussian
disturbances according to the logistic distribution in Equation (4), with the same aim of simulating
difficult conditions.
The methods were compared with three previously published methods:• LiNGAM estimated using ICA, as proposed in the original paper introducing LiNGAM by
Shimizu et al. (2006).5• DirectLiNGAM, specifically the kernel-based version proposed by Shimizu et al. (2011).6• In case of skewed data, we used the measure proposed by Dodge and Rousson (2001), given
in Equation (14).
The LiNGAM methods were implemented using the software found on the authors’ web sites.
We computed different performance indices for the methods. For acyclic models, we computed
5. Since basic FastICA, which is an integral part of the method, has convergence problems with the basic tanh nonlinearity in the case of a small sample size, we used the stabilized version by Hyvärinen (1999) obtained in the standard
FastICA package by the option “stabilize”. For skewed data, we used the skewness as a measure of non-Gaussianity.
6. We did not include the earliest version of DirectLiNGAM proposed by Shimizu et al. (2009) in the comparison
because in later simulations by Sogawa et al. (2010); Hyvärinen (2010), its performance was found clearly inferior
to that of the kernel-based version of DirectLiNGAM.

129

H YV ÄRINEN AND S MITH

1. The Spearman rank-correlation coefficient between the causal ordering given by the method
and the true ordering.
2. The percentage of connections for which a method correctly estimated the direction, considering only connections existing in the data-generating process. Here, the point was to look at
the abilities of the methods to find the directions locally, and thus the global ordering givenby the method was not used (except for DirectLiNGAM which essentially only computes a
global ordering and derives local ordering from that). For the ICA-based LiNGAM, we computed the measure sign(|bi j | − |b ji |) and used it in the same way as the signs of the pairwise
measures.
3. The percentage of data sets for which a method correctly estimated the first variable in thecausal ordering, that is, the variable with no parents.
For cyclic models, the comparison was based on the second measure only, since the other two are
not well-defined. Furthermore, we computed the CPU time needed for the computations.
Unless otherwise mentioned, the connection matrices were generated completely randomly, giving a fully connected DAG. The non-zero coefficients in the acyclic B had a uniform distribution in
the union of the intervals [−0.6, −0.2] and [0.2, 0.6].
6.1 Simulation 1: Sparse Influences
In the first, basic simulation, sample size and data dimension were varied so that there were in total
four different scenarios:
1. n = 5, T = 100, fully connected DAG
2. n = 2, T = 100, fully connected DAG
3. n = 5, T = 200, fully connected DAG
4. n = 5, T = 400, fully connected DAG
The disturbances had logistic distributions, with standard deviations equal to one. 2,000 data sets
were generated for each scenario; however, for DirectLiNGAM and ICA-LiNGAM only 1,000 wereused due to excessive computational demands.
To estimate the model, we used the following methods proposed above. First, the maximum
entropy approximation to the likelihood ratio in (3) was used in DirectLiNGAM with deflation.
Second, the LR first-order approximation matrix (18) was used in DirectLiNGAM with the nonlinearity g(u) = tanh(u) and with deflation. Third, the nonlinear correlations in (18) were used toestimate the causal ordering without any deflation, simply by locating the minimum of the row sumsof that matrix, removing the corresponding rows and columns, and so on, as described at the end ofsection 3.
See Figure 3 for the results. Typically, the tanh-based likelihood ratio approximation (“tanh”)
with deflation was the best. The method without deflation (“nodf”) gives, by definition, the same result for the total causal directions correct and first variable found, but looking at the rank-correlations,
we see that it is typically the second best. The maximum entropy approximation is usually the third
best. ICA-based LiNGAM is usually fourth but when there is more data, it can have very good130

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

performance. The (kernel-based) DirectLiNGAM (“kdir”) is typically last, although not necessarily
worse than ICA-based LiNGAM.
Regarding computational load, the methods proposed here are one to two orders of magnitude
faster than the others.
6.2 Simulation 2: Sparse Influences with NoiseIn the second simulation, we tested the noise-tolerance of the algorithms. The data dimension wasvaried from n = 2 to n = 8 and fully connected DAGs were used as above. The sample size was set
to T = 10, 000, which means we are now analyzing the statistical consistency7 of the method only
and neglecting random effects by taking a very large sample size. The noise was Gaussian and had
unit variance. The performance indices and algorithms are as in the first simulation. The results areshown in Figure 4. We can see that the tanh-based approximation is clearly the best, as predictedby our cumulant-based analysis. ICA-based LiNGAM, the maximum entropy approximations, andespecially kernel-based DirectLiNGAM seem to be more sensitive to noise.
6.3 Simulation 3: Skewed Influences
In the third simulation, we tested the performance of the methods with skewed data. We used thenonlinear correlation based on the third order cumulant (“skew”), introduced in Section 2.9, as wellas the robust measure in Equation (15), denoted by “skw2”.
We used two different skewed distributions for the disturbances. In both cases, the data was
obtained from a Gaussian mixture. One of the Gaussian distributions in the mixture had zero mean
and unit variance, while the other had mean equal to three and unit variance. The two distributions
we generated were distinguished by the amount of data points drawn from the two Gaussians. In
the first case (“pdf 1”), the “outlying” distribution with mean three generated 20% of the data, whilein the second case (“pdf 2”), it generated only 5%. Thus, pdf 2 was quite sparse whereas pdf 1 wasnot. We would then expect sparsity-based methods to work well with pdf 2 but not very well with
pdf 1. The data dimension were to n = 2, n = 5 and sample sizes T = 100, 200, respectively. DAGs
were generated to be fully connected.
The results are shown in Figure 5. We see that all the methods have very similar performance,
except the Dodge-Rousson measure which was somewhat worse. However, the computational loads
are very different, our two likelihood ratio approximations being faster than the earlier LiNGAM
methods by at least an order of magnitude.
6.4 Simulation 4: Skewed Influences with Noise
We further conducted a simulation with observational noise added to the skewed data. Again, we
fixed the sample size to T = 10, 000 and the noise variance to two (larger than above since thesemethods seem to be more tolerant to Gaussian noise), while the dimension and the skew data distribution were varied. We used only the skewed and sparse pdf 2. The results are in Figure 6. Here,
we start seeing clear differences in the statistical performances of the methods. In line with ourtheoretical analysis, the skewness cumulant-based method is the most resistant to noise. The robust
skewed LR approximation in Section 2.9.2 is second.
7. That is, convergence in the limit of infinite sample.

131

H YV ÄRINEN AND S MITH

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

tanh

nodf mxnt ICA
First variable found

kdir

0.5

tanh

nodf mxnt ICA

kdir

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=2,T=100
n=5,T=100
n=5,T=200
n=5,T=400

tanh

10

nodf mxnt ICA

1

kdir

tanh

nodf mxnt ICA

kdir

Figure 3: Simulation 1. Results of basic simulation with sparse, non-skewed data without noise.
Top left: Mean of rank-correlation coefficients between the estimated causal orderingand the true ordering. The error bars are standard errors of the mean. Top right: The
proportion of (really existing) connections for which the method estimated the direction
correctly (chance level is 50%). Bottom left: The proportion of data sets for which themethod estimated the first variable in the causal ordering correctly, that is, the variablewith no parents. Bottom right: Computation times of one run of the different algorithms
in milliseconds; note the logarithmic scale. Different colours are different data-generating
scenarios. The algorithms used are as follows:“tanh”: LR approximations in (18) based on tanh nonlinearity, combined with deflationin DirectLiNGAM;“nodf”: no deflation in likelihood ratio approximations, that is, ordering based on the LR
approximation matrix in (18) without any recomputation of the matrix;“mxnt”: maximum entropy approximation in (3) for likelihood ratios;“ICA”: LiNGAM estimated by ICA;
“kdir”: kernel-based DirectLiNGAM.

132

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

tanh

nodf mxnt ICA
First variable found

kdir

0.5

tanh

nodf mxnt ICA

kdir

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=2
n=3
n=5
n=8

tanh

10

nodf mxnt ICA

1

kdir

tanh

nodf mxnt ICA

kdir

Figure 4: Simulation 2, with noise. Legend as in Figure 3, and with T = 10, 000. The noise standard
deviations were all equal to one.

6.5 Simulations 5 and 6: Two-Stage Approach and Sparse Graphs
Next we investigated the utility of the two-stage approach of Section 3.3. We generated sparsegraphs only. The graphs were based on a simple “serial” structure x1 → x2 → . . . → xn with arandom connection strength in the same range as above. We further added 0, 1, or 2 connections
in random locations in the graph (preserving the DAG structure), the number of connections having
equal probabilities for the three values. The data sizes were 500, 900, 900, 900 and the numberof variables 5, 9, 15, 20, respectively. We used higher dimensions than above because otherwise
the networks could not be very sparse. In the testing for the existence of connections, we set thefalse-positive rate to α = 0.01 without correction for multiple testing, as motivated above.
In Simulation 5, we used sparse, non-skewed (logistic) influences, and in Simulation 6, skewed
influences as in Simulation 3. To add more realism to the simulations, we also added noise to thedata. The noise standard deviations were 0.2 in Simulation 5 and 0.6 in Simulation 6.
The results for Simulation 5 are shown in Figure 7. We can see that the two-stage method has
a performance which compares quite favourably with the other methods: ICA-LiNGAM and DirectLiNGAM perform quite badly with these combinations of sample size and dimension. Note that133

H YV ÄRINEN AND S MITH

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

skew skw2 mxnt ICA kdir D−RFirst variable found

0.5

skew skw2 mxnt ICA

kdir D−R

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=2,T=100,pdf 1
n=5,T=200,pdf 1
n=2,T=100,pdf 2
n=5,T=200,pdf 2

skew skw2 mxnt ICA

101

kdir D−R

skew skw2 mxnt ICA kdir D−R

Figure 5: Simulation 3, with skewed data. Legend as in Figure 3, with the following new algorithms:“skew”: cumulant-based LR approximation in (10), combined with deflation in DirectLiNGAM;“skw2”: the robust LR approximation proposed in Section 2.9.2; and“D-R”: the measure by Dodge and Rousson (2001).

for “total causal directions correct”, the two-stage method has, by definition, the same performanceas “tanh” and “nodf”. In fact, if our interest in only to discover the directions without bothering toestimate which variables are connected, or we are given perfect prior knowledge on which variables
are connected, there is in fact no need to do the pruning in the first stage of the method.
So, the utility of the new method (“icth”) is mainly seen in the mean rank correlations plot:
There is a modest improvement. The point here is that knowledge of which variables are connected
improves the estimation of the causal ordering (DAG structure) by showing which directionalities
should be used when pooling their information together, and which directionalities should be discarded (because the variables are not connected at all).
134

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

skew skw2 mxnt ICA kdir D−RFirst variable found

0.5

skew skw2 mxnt ICA

kdir D−R

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=2
n=3
n=5
n=8

skew skw2 mxnt ICA

101

kdir D−R

skew skw2 mxnt ICA kdir D−R

Figure 6: Simulation 4, with skewed data with noise. Legend as in Figure 5.

Interestingly, all the methods proposed in this paper are clearly superior to the methods proposed earlier (ICA-based LiNGAM and kernel-based DirectLiNGAM). Thus, the main utility ofthe present framework may indeed be in estimating directionality in sparse networks.
We carried out the same simulation for skewed influences using the skewed pdf 1. Results are
in Figure 8. When looking at methods using the same causality measure (“skew” vs. “icsk”, and“skw2” vs. “ics2”), we see that the pruning methods are better in terms of the mean rank correlations.
However, the maximum entropy method without pruning is actually the best.
6.6 Overview of Simulations 1–6
To provide a succinct overview of the simulations reported above, we averaged the performanceindices over the different scenarios (taking into account only scenarios in which the algorithm tookpart). Furthermore, we divided the simulations into three groups: basic data (simulations 1 and 2),
skewed data (simulations 3 and 4) and sparse connections either with sparse data (simulation 5) or
skewed data (simulation 6). We further averaged the performance indices inside these groups.
The results are shown in Figure 9.
135

H YV ÄRINEN AND S MITH

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

tanh nodf mxnt icth ICA
First variable found

kdir

0.5

tanh nodf mxnt icth ICA

kdir

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=9,T=500
n=9,T=900
n=15,T=900
n=20,T=900

tanh nodf mxnt icth ICA

10
1

kdir

tanh nodf mxnt icth ICA kdir

Figure 7: Simulation 5, with the two-stage pruning method and only sparse graphs. Legend as in
Figure 3, but now including the new algorithm “icth” which prunes the graph based on
inverse covariance and then estimates the directions using the same method as “tanh”.
(Note that only “icth” uses information on the pruned inverse covariance, other methods
are as in Simulation 1.)

6.7 Simulation 7: More Variables than Observations
Next, we considered the case where there are more variables than observations, or at least the
number of variables is equal to the number of observations. We considered four scenarios, with n
ranging from 100 to 200 and T ranging from 100 to 400. In preliminary simulations, it turned outthat the problem was too difficult for logistic disturbances, so we used Laplacian disturbances here.
We only attempted to estimate the first two variables and not the whole causal ordering. Thevery first variables in the causal ordering can be considered to be the exogenous ones and thus
finding them is of special interest (Sogawa et al., 2011). We only used three of the new proposed
methods because none of implementations of the existing LiNGAM methods was such that it could
readily be used for this case.
The results are shown in Figure 10. While the performance of the methods is not very good, it
is very much above chance level (which would be 0.01 or less for finding the first variable). It is
136

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

skewskw2mxnt icsk ics2 ICA kdir
First variable found

0.5

skewskw2mxnt icsk ics2 ICA kdir
Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=9,T=500
n=9,T=900
n=15,T=900
n=20,T=900

101skewskw2mxnt icsk ics2 ICA kdirskewskw2mxnt icsk ics2 ICA kdir

Figure 8: Simulation 6, with skewed data, the two-stage pruning method and only sparse graphs.
Legend as in Figure 5, but now including the new algorithm “icsk” which prunes thegraph based on inverse covariance and estimates the directions based on the skewness
cumulant, and “ics2” which uses the robust skewness measure.

interesting to note that here the first-order approximation of likelihood is more than 100 times fasterthan the maximum entropy approximation.
6.8 Simulation 8: Cyclic Graphs
To test the new framework in the case of cyclic graphs, we created cyclic graphs by a simple ring
structure: x1 → x2 , . . . , xn → x1 . Further connections (0, 1, or 2) were added in random locations as
in Simulation 5 above. Such data were created according to the generating model in Section 4. We
further added noise with standard deviation 0.2. The dimensions of the data and the sample sizes
were as in Simulations 5 and 6. The influences had logistic distributions.
The only methods we compared were ICA-based LiNGAM and our two-stage pruning methods,
since the DirectLiNGAM methods cannot be used in the cyclic case. The results are shown in
Figure 11. The first observation is that both methods performed relatively well, obtaining 70%90% percent of the directions right. Our new method is slightly better than ICA-based LiNGAM.
137

H YV ÄRINEN AND S MITH

Basic sparse data

Skewed data

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

tanh nodf mxnt ICA kdir
Sparsely connected data

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

tanh nodf mxnt icth ICA

0

kdir

skew skw2 mxnt ICA kdir D−R
Sparsely connected, skewed data

skewskw2mxnt icsk ics2 ICA kdir

Figure 9: Overview of Simulations 1–6. Median correlations (blue, solid) and average directions
correct (green, dashed) are plotted averaged over different scenarios and similar simulations.

It should be emphasized here that our method assumes that there are no self-loops, so there is noindeterminacy in the results, as shown in Section 4.
6.9 Simulation 9: Nonlinear Relations
Finally, we performed simulations on the nonlinear model. We generated data from a modelx2 = αsign(x1 )|x1 |γ + d(25)where both x1 and d were standardized Gaussian. The exponents γ were given values 0.5 and 2, andthe parameter α was randomly drawn between 0.5 and 1.5. The sample sizes were either T = 200
or T = 500.
We then fitted the nonlinearity of the same functional form (25), that is, using the parameters α
and γ, to the data with a least-squares fit, and estimated the causal direction using the criterion in(22), or the criterion in (24). (Thus, we did not use a nonparametric model of the nonlinearity. See
Section 8 for estimation with non-parametric nonlinearities.) For comparison, we used the methods
138

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM SFirst variable foundTwo first variables found

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2
0

0.2

n=100,T=100
n=200,T=100
n=200,T=200
n=400,T=200

tanh

nodf

0

mxnt

tanh

nodf

mxnt

Computation time
100000
10000
1000
100
10
1

tanh

nodf

mxnt

Figure 10: Simulation 7, with more variables than observations. Legend as in Figure 3. Rankcorrelations and causal directions correct are omitted because we only computed thefirst two variables for lack of computation time.

tanh and maxent introduced above in a purely linear way (i.e., not fitting the nonlinear function
above, but just a linear function exactly as in previous simulations), to see if linear methods are able
to cope with this data.
Furthermore, we used the criterion of the original method by Hoyer et al. (2009), based on the
HSIC independence test by Gretton et al. (2008) of x (resp. y) and the residual in the regression ofy on x (resp. of x on y). This was implemented by code provided by A. Gretton,8 using the default
setting for the kernel width.
The results are shown in Figure 12. Our likelihood ratio methods both performed relatively
well, although the independence-based method by Hoyer et al. (2009) was arguably better than our
maximum entropy method. However, the HSIC-based method was 10-100 times slower due to theuse of kernel methods. The linear methods did not perform well at all.

8. Downloaded from http://www.gatsby.ucl.ac.uk/˜gretton/indepTestFiles/indep.htm.

139

H YV ÄRINEN AND S MITH

Total causal directions correct
Computation time

1
100000

0.9

10000

0.8

1000

0.7

100

0.6

10

0.5

icth

1

ICA

icth

ICA

Figure 11: Simulation 8, with cyclic sparse graphs. Legend (sample sizes and dimensions) as in
Figure 7.

Mean rank correlations
Computation time

1
100000

0.8

10000

0.6

1000

0.4

100

0.2

10

0

tanh

1

mxnt nlme hsic madtanhmxnt nlme hsic mad

Figure 12: Simulation 9, with nonlinear model. The new algorithms are “nlme”, the proposed
likelihood ratio method extended to the nonlinear case using maximum entropy approximation in (22); “mad”, a simplified and robustified approximation of the likelihood
ratio in (24); “hsic”, the original nonlinear method using independence (Hoyer et al.,2009). Blue: γ = 0.5, T = 200, Green: γ = 2, T = 200, Red: γ = 0.5, T = 500, Cyan:γ = 2, T = 500.

6.10 Simulation 10: Misspecified Disturbances
We further performed simulations in which the model is misspecified. First, we considered Simulation 1 with the following change: the disturbances generating the data had Laplacian distributions.
Everything else was identical to Simulation 1, including the assumed log-pdf’s and nonlinearities.
Thus, the distribution of the disturbances was not exactly known, and was misspecified in the estimation. We also added the basic skewness method in the set of algorithms.
140

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

tanh skew nodf mxnt ICA
First variable found

kdir

0.5

tanh skew nodf mxnt ICA

kdir

Computation time

1
100000

0.8

10000

0.6

1000

0.4
0.2
0

100
n=2,T=100
n=5,T=100
n=5,T=200
n=5,T=400

tanh skew nodf mxnt ICA

10
1

kdir

tanh skew nodf mxnt ICA kdir

Figure 13: Simulation 10. Like Simulation 1 but with Laplacian disturbances used in generating
the data, and the “skew” method added.

The results are in Figure 13. We see that the performance of most methods is actually better.
This was expected in light of the theory of ICA, where it is well-known that if the actual data is more
non-Gaussian than assumed in the estimation method, this is not a problem for most methods, andonly increases the performance of the method compared to the case of less non-Gaussian data. In
fact, the reason why we used the logistic distribution in generating data in many of the simulations
above was in order to make the problem more difficult. On the other hand, the reason for using thelogistic distribution in the algorithms is that it is widely used in ICA and has been empirically foundto work well, partly due to the fact that its log-pdf is smooth, unlike many other super-Gaussian
log-pdf’s including the Laplacian.
Of course, if the non-Gaussianity is completely misspecified in the estimation method, estimation with fixed nonlinearities will inevitably fail. This is why the skewness method was hardly above
chance level.
141

H YV ÄRINEN AND S MITH

6.11 Simulation 11: Latent Variables
We conducted a further simulation to gain some insight into the robustness of the different methods
to the existence of latent variables. We first created data x0 as in Simulation 1, with n = 4, T = 500.
Then, we added a latent variable to the data as
x = x0 + αbs̃
where s̃ is a latent variable with a standardized logistic distribution, b is a weight vector with elements drawn from a standardized Gaussian distribution, and α is the general strength of the latentvariable, which took the values [0, 0.25, 0.5, 1] in the different scenarios. (The value of α = 0 effectively means no latent variables and is provided for comparison.) The latent variable s̃ violates theassumption of LiNGAM of having only one (independent) external input for each variable xi .
The results are in Figure 14. Basically, we see that the latent variable deteriorates the performance of all the algorithms quite uniformly. It does not seem that any of the algorithms would be
more resistant, or more sensitive, to latent variables than the others.
Recently, the framework presented here was generalized to a model including Gaussian latentvariables by Chen and Chan (2012).

7. Experiments on Simulated fMRI Data
Since causal discovery experiments on real data are very difficult to validate, we use here brain imaging data which has been simulated using state-of-the-art biophysical models (Smith et al., 2011).
7.1 Simulation of fMRI Data
The simulations are described in detail by Smith et al. (2011); here we give a short summary.
Networks of varied complexity were used to simulate fMRI timeseries. The simulations were based
upon the dynamic causal modelling (DCM) forward model (Friston et al., 2003). DCM uses thenonlinear balloon model (Buxton et al., 1998) for the vascular dynamics, that is, the connection
between the neural activities and the measured signal, sitting above a simple neural network modelof the neural dynamics. Estimating causality from fMRI data is particularly challenging as thesignal-to-noise ratio is relatively poor, fMRI timeseries are fairly Gaussian, and the number of
timepoints is generally in the low hundreds.
We defined a number of nodes, which corresponded to brain regions. First, we generated theexternal inputs to the nodes, ui , which are not quite the same as the external influences in the SEM,
although related. They were binary (activity is “up” or “down”) and generated using a Poissonprocess that controls the likelihood of switching the state. Neural noise of standard deviation 1/20of the difference in height between the two states was added. The mean durations of the states
were 2.5s (up) and 10s (down), with the asymmetry representing longer average “rest” than “firing”
durations.
The neural activities zi were then simulated using the DCM neural network model, as definedbyż = σAz + Mu
where A is the matrix defining network dynamics and M contains the weights controlling how
the external inputs feed into the network (often just the identity matrix). The off-diagonal terms
142

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Mean rank correlations

Total causal directions correct

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0

tanh

nodf mxnt ICA
First variable found

kdir

0.5

tanh

nodf mxnt ICA

kdir

Computation time

1
100000

0.8

10000

0.6

1000

0.4

100

0.2

10

0

tanh

nodf mxnt ICA

1

kdir

tanh

nodf mxnt ICA

kdir

Figure 14: Simulation 11. Like Simulation 1, with n = 4, T = 500, but with a latent variable added.
The four scenarios (curves) correspond to different strengths of the latent variable, starting with zero strength in blue curve.

in A determine the network connections between nodes, and the diagonal elements are all set to-1, to model within-node temporal decay; thus σ controls both the within-node (neural) temporal
inertia/smoothing and the neural lag between nodes.
A central problem in fMRI is that the measured signal does not directly correspond to z. To simulate this, each node’s neural timeseries zi was fed through the nonlinear balloon model for vascular
dynamics responding to changing neural demand. The balloon model parameters were in general
set according to the prior means in DCM. However, it is known that the haemodynamic processes
vary across brain areas and subjects, resulting in different lags between the neural processes and the
BOLD data, with variations of up to at least 1s (Handwerker et al., 2004; Chang et al., 2008). We
therefore added randomness into the balloon model parameters at each node, resulting in variations
in HRF (haemodynamic response function) delay of standard deviation 0.5s. Finally, thermal white(measurement) noise of standard deviation 0.1–1% (of mean signal level) was added.
143

H YV ÄRINEN AND S MITH

Thus, we obtained the measured fMRI signals. They were sampled with a sampling intervalof 3s (in most simulations), corresponding to a typical time of repetition (TR) in brain imaging
literature.
The simulations comprised 50 separate realisations (or “subjects”), all using the same simulation
parameters, except for having independently generated external inputs and different HRF parametersat each node (as described above); furthermore, the connection strengths were slightly perturbatedfor each subject. Each “subject’s” data was a 10-minute fMRI session (200 timepoints) in many of
the simulations.
For a summary of the specifications for the 28 simulations see Table 1.
Sim

n

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

5
10
15
50
5
10
5
5
5
5
10
10
5
5
5
5
10
5
5
5
5
5
5
5
5
5
5
5

length(mins)
10
10
10
10
60
60
250
10
250
10
10
10
10
10
10
10
10
10
10
10
10
10
10
10
5
2.5
2.5
5

TR
(s)
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00
0.25
0.25
3.00
3.00
3.00
3.00
3.00
3.00
3.00
3.00

noise
(%)
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
1.0
0.1
1.0
0.1
1.0
0.1
0.1
1.0
0.1
0.1
0.1
1.0
1.0
0.1
0.1

HRF std
(s)
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.0
0.5
0.0
0.5
0.5
0.5
0.5
0.5
0.5
0.5
0.5

other factors

shared inputs
shared inputs
global mean confoundtimeseries mixed with each othernew random timeseries mixed inbackwards connections
cyclic connections
stronger connections
more connections

neural lag=100ms
neural lag=100ms
2-group test
nonstationary connection strengths
stationary connection strengths
only one strong external input

Table 1: Summary of the 28 fMRI simulations’ specifications (from Smith et al., 2011)

144

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

7.2 Estimation Methods for Simulated fMRI Data
We used pairwise measures with three different nonlinearities: tanh, skewness, and the robust measure of skewness. We did not estimate the existence of connections at all since that is not the main
topic of the paper: We only looked at the estimated directionalities for those connections whichreally existed in the simulated data.
Since the skewness of the data was mainly positive, we used this prior information of positive
skewness skewness-based measures. In other words, we skipped the skewness correction in (13).
For comparison, we used two methods by Patel et al. (2006) which were the most successfulof the many methods tested by Smith et al. (2011), as well as basic ICA-based LiNGAM (Shimizu
et al., 2006) which was applied on the whole data (not pairwise).
7.3 Results on Simulated fMRI Data
The goal is thus to recover the directionalities defined by the non-zero entries of the neural dynamics
matrix A by estimating the directionalities given by B in our SEM. We evaluated the results using
the same measures as Smith et al. (2011) to allow for direct comparison.
Our evaluations looked at the distribution of correctly estimated connections over the 50 simulated subjects. We concentrate here on evaluating methods for single subject (single session) data
sets, and only utilise multiple subjects’ data sets in order to characterise variability of results acrossmultiple random instantiations of the same underlying network simulation. This is in contrast to theapproach by Ramsey et al. (2011) who estimated the network over random subsets of 10 subjects,
which is an easier task, at least if the subjects are not very different.
The raw connection strengths bi j were converted into z-scores in order to make the plots morequalitatively interpretable, as the connection strengths are then more comparable across the different
methods. The conversion from raw connection strengths to z-scores was achieved by using a null
distribution of connection strengths, obtained by feeding in truly null timeseries data into each of theestimation methods. The null data was created by testing for connections between timeseries fromdifferent subjects’ data sets, which have no causal connections between them (i.e., we randomly
shuffled the subject labels for each node in the network). See Smith et al. (2011) for details. To
specifically look at estimated directionalities, we use the higher of the two directions’ measures tobe the estimated connection strength.
The results are shown in Figs 15-16. The distributions are over all 50 simulated “subjects” andover all correct network edges; higher is better. Note, however, that this plot does not take into
account the false positives, that is, the values estimated in the network matrix that should be empty,
and concentrates exclusively on the estimation of causal directions. The plots, known as “violin
plots”, are simply (vertically-oriented) smoothed histograms, reflected in the vertical axis for better
visualisation.
We see that the pairwise methods perform much better than Patel’s measures or ICA-based
LiNGAM on all the simulations. (The comparison to ICA-based LiNGAM may not be entirely fairsince it estimates more than just directionalities.) In fact, our methods perform extremely well inmost simulations. In all the simulations, the pairwise measures are the best, although in two cases
the performances of all methods are so close to chance level that any comparison is difficult. The
results are not very good in the following cases:• Simulation 13 which has backwards connections (i.e., both x → y and y → x) which is not surprising since it is against the basic philosophy of our modelling. However, the performance is145

H YV ÄRINEN AND S MITH

Caption for Figs. 15 and 16 on pages 147 and 148: The z-scores of the different measures used todetermine the directionality, computed over subjects and connections, are shown as violin plots (i.e.,
histograms rotated to be horizontal and made symmetric). If the directions are found completely
correctly, the violin plots are concentrated at the top. The blue dots show the the percentage ofcorrectly estimated directions. First, we have three pairwise methods, and for comparison, two
methods by Patel, as well as ICA-based LiNGAM. Each panel is one simulation.

clearly better than chance, which shows that our method is able to find the dominant directionsome of the time.
• Simulation 22 which has nonstationary connection strengths, which violates another basic
assumption of the model.
• Simulation 24 in which one of the inputs is strongly dominant. This is presumably because
the effective signal-to-noise ratio is too poor for many of the connections.
• Simulations 25-27 in which the number of data points is smaller (recording length is shorter),
performance being close to chance level for all methods.
Among the pairwise measures, there is no clear winner. However, the robust skewness measure
is most often the best (in those cases where a clear difference can be seen), and never much worsethen the other two.

8. Nonlinear Causal Discovery on Real Data
Finally, we applied our nonlinear methods on the Tübingen-UCI cause-effect data set9 which consists of real measurements in which the true direction of causation in known. We used a total of 81
data sets, consisting of the subset of those data that had exactly two variables. In addition to our
two new nonlinear methods, we applied the original HSIC-based methods by Hoyer et al. (2009),
as well as the linear likelihood ratio with maximum negentropy approximation of Section 2.3. The
relations in this data set are often quite nonlinear, and the linear methods are hardly above chancelevel (results not shown except for one method below), so we concentrate on the nonlinear methods
here.
The nonlinear regression was performed by first fitting a least-squares regression curve using a
Gaussian process as implemented in the fit gp package by J. Mooij, based on code by C. E. Rasmussen and H. Nickisch. We used only the first 1,000 data points due to the excessive computational
complexity of HSIC.
The results are shown in Table 2. The linear method, as well as our basic nonlinear method
using maximum entropy were hardly better than chance. The method by Hoyer et al. (2009) was
close to 62%. On the other hand, our simplest approximation using mean absolute deviation was
69% correct.
Presumably, one reason for the weak performance of our nonlinear method using maximum
entropy approximations was that many of the data sets have strong outliers. The MAD-based objective in Equation (24) is quite robust against them (although the nonlinear regression method was not
9. Data set can be found at http://webdav.tuebingen.mpg.de/cause-effect/.

146

10

100

5

75

0

50

−5

25

−10

0

ICA−LiNGAM

Patel’sPatel’s bin.75

PW−LR tanh

PW−LR r skew

PW−LR skew

% directions correctcausality (Zright − Zwrong)

ICA−LiNGAM

Patel’sPatel’s bin.75% directions correct0

PW−LR tanh

−10

PW−LR r skew

25

PW−LR skew

75

0

50

−5

25

−10

0

ICA−LiNGAM

5Patel’s100

Patel’s bin.75

10% directions correct

Simulation 6 (10 nodes, 60 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )PW−LR tanh

% directions correct

Simulation 14 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , cyclic connections)causality (Zright − Zwrong)% directions correct

Simulation 11 (10 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , bad ROIs (timeseries mixed with each other))−5

Simulation 9 (5 nodes, 250 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , shared inputs)
10

100

5

75

0

50

−5

25

−10

0

% directions correct

050

Simulation 12 (10 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , bad ROIs (new random timeseries mixed in))
10

100

5

75

0

50

−5

25

−10

0

% directions correct−100

ICA−LiNGAM

2575

ICA−LiNGAM−5

5Patel’s

50

100

Patel’s bin.75

010Patel’s

75

0Patel’s bin.75

5

ICA−LiNGAM

100

Patel’s bin.75

10% directions correct

Simulation 8 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , shared inputs)−10

PW−LR r skew

025

PW−LR tanh

−10

−5

PW−LR r skew

25

50

PW−LR tanh

−50PW−LR skew

50

75

PW−LR skew

0

ICA−LiNGAM

75

Patel’s bin.75

5Patel’s

100

PW−LR r skew

10

% directions correct

Simulation 5 (5 nodes, 60 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )causality (Zright − Zwrong)

ICA−LiNGAM

Patel’s bin.75

0Patel’s−10

PW−LR r skew

25

5

Simulation 3 (15 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )causality (Zright − Zwrong)

ICA−LiNGAM

Patel’s bin.75

Patel’s

PW−LR r skew

PW−LR tanh

PW−LR skew

−5

PW−LR tanh50

PW−LR skew

% directions correct

Simulation 13 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , backwards connections)

0100

PW−LR r skew

0

75

10

PW−LR skew

−105

% directions correct

25

100

causality (Zright − Zwrong)−5% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 10 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , global mean confound)

10

ICA−LiNGAM

0

Simulation 2 (10 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )Patel’s bin.75−100Patel’s

25−10

PW−LR r skew

−5

% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 7 (5 nodes, 250 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )25

Patel’s

0−5PW−LR r skew

−10

50

PW−LR tanh250PW−LR skew

−5

% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 4 (50 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )75

PW−LR tanh

05PW−LR skew

−10

100

PW−LR tanh25

10

PW−LR skew

−5

% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 1 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )causality (Zright − Zwrong)0causality (Zright − Zwrong)−10

causality (Zright − Zwrong)25

causality (Zright − Zwrong)−5

causality (Zright − Zwrong)

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Simulation 15 (5 nodes, 10 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s , stronger connections)Figure 15: Results on simulated fMRI data, first half. See page 146 for caption.

147

50

−5

25

−10

0

10

100

5

75

0

50

−5

25

−10

0

ICA−LiNGAM

Patel’sPatel’s bin.75

PW−LR tanh

PW−LR r skew

PW−LR skew

% directions correctcausality (Zright − Zwrong)

25

−10

0

% directions correct−5

ICA−LiNGAM

50

Patel’s0

Patel’s bin.75

75

PW−LR tanh5

PW−LR r skew

100

PW−LR skew

75

0

50

−5

25

−10

0

ICA−LiNGAM

5Patel’s100

Patel’s bin.75

10% directions correct

Simulation 21 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , 2−group test)

PW−LR tanh

% directions correct% directions correctcausality (Zright − Zwrong)

10

Simulation 24 (5 nodes, 10 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s , only one strong external input)
10

100

5

75

0

50

−5

25

−10

0

Simulation 27 (5 nodes, 2.5 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s )% directions correct

Simulation 26 (5 nodes, 2.5 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )0% directions correct

0−10

ICA−LiNGAM

75

25

Patel’s5−5Patel’s bin.75

100

ICA−LiNGAM

10% directions correct

Simulation 23 (5 nodes, 10 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s , stationary connection strengths)

50

PW−LR r skew

00

PW−LR tanh

−10

75

PW−LR r skew

25

5

PW−LR skew

−5

100

PW−LR skew

50

% directions correct

0

ICA−LiNGAM

75

Patel’s bin.75

5Patel’s

100

PW−LR r skew

10

causality (Zright − Zwrong)

0

ICA−LiNGAM−10Patel’s bin.75

25

Patel’s−5

PW−LR r skew

50

10

Simulation 18 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.0s )causality (Zright − Zwrong)

ICA−LiNGAM

Patel’s bin.75

Patel’s

PW−LR r skew

PW−LR tanh

PW−LR skew

0

PW−LR tanh75

PW−LR skew

5

Simulation 20 (5 nodes, 10 minute sessions, TR=0.25s, noise=0.1%,
HRFstd=0.0s , neural lag=100ms)

causality (Zright − Zwrong)

Simulation 25 (5 nodes, 5 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s )100

Patel’s bin.75

010

ICA−LiNGAM−10Simulation 17 (10 nodes, 10 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s )Patel’s bin.75

25

0Patel’s−5% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 22 (5 nodes, 10 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s , nonstationary connection strengths)−10

PW−LR r skew

025

Patel’s−10−5PW−LR r skew

25

50

PW−LR tanh−5

% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 19 (5 nodes, 10 minute sessions, TR=0.25s, noise=0.1%,
HRFstd=0.5s , neural lag=100ms)

0

PW−LR skew

075

PW−LR tanh

−105PW−LR skew

25

100

PW−LR tanh−5

% directions correct

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

Simulation 16 (5 nodes, 10 minute sessions, TR=3.00s, noise=1.0%,
HRFstd=0.5s , more connections)

10

PW−LR skew

% directions correct0

causality (Zright − Zwrong)−10

causality (Zright − Zwrong)25

causality (Zright − Zwrong)−5

causality (Zright − Zwrong)

50

ICA−LiNGAM

0Patel’s bin.75

75

Patel’s

5

PW−LR r skew

100

PW−LR tanh

10

PW−LR skew

causality (Zright − Zwrong)

H YV ÄRINEN AND S MITH

Simulation 28 (5 nodes, 5 minute sessions, TR=3.00s, noise=0.1%,
HRFstd=0.5s )Figure 16: Results on simulated fMRI data, second half. See page 146 for caption.

148

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Method
mxnt
hsic
nlme
mad

% correct
50.6
61.7
53.1
69.1

Table 2: Nonlinear models applied on the Tübingen-UCI data set. The algorithms are “nlme”,
the nonlinear likelihood ratio with maximum entropy approximation; “mad”, the approximation using mean absolute deviations. For comparison: “hsic”, the original nonlinear
method using the HSIC measure of independence (Hoyer et al., 2009); and “mxnt”, thelinear method using the maximum entropy approximation.

made robust). The maximum entropy approximation might be greatly improved if the estimation of
the variances used in the maximum entropy objective were made robust against outliers. Furthermore, both of our methods might be improved if the nonlinear fitting used a robust criterion insteadof least squares.

9. Conclusion
We proposed very simple measures of the pairwise causal direction based on likelihood ratio tests
and their approximations. We started with general measures based on entropy approximation which
can accommodate different kinds of distributions, in Equations (2) and (3). Assuming that priorknowledge is available, we can develop more specific methods; for sparse variables we propose (5)and for skewed variables (17). We further showed how the measures can be extended to cyclic andnonlinear models. The different measures are recapitulated in Table 3.
We also showed how the pairwise measures can be used to estimate the whole Bayesian network
in two ways. This is possible either in the DirectLiNGAM framework, or by a two-stage method
based on first estimating the existence of the connections and then orienting them using the pairwise
measures.
We also proposed a cumulant-based version of the nonlinear correlations. It was shown thatthe cumulant gives the correct pairwise direction. This shows the utility of using cumulants intheoretical analysis, and gives an intuitive interpretation of a new kind of cumulant. The cumulantbased analysis also indicated the noise-robustness of the nonlinear correlation methods, which was
confirmed in the simulations. However, in practice the cumulant-based methods may suffer fromsensitivity to outliers and thus their utility may be mainly in theoretical analysis.
The proposed measures seem to be particularly useful in the case where the number of data
points is small compared to the dimension of the data, or the data is noisy. In such a case, the statistical performance of our methods is clearly superior to ICA-based LiNGAM and, to a lesser extent,
DirectLiNGAM. The new methods are also computationally much faster than DirectLiNGAM. The
importance of estimating causal networks with few data points has been recently highlighted bySmith et al. (2011) in the context of brain imaging. In fact, applied to the simulations by Smith et al.
(2011), the new pairwise measures were clearly better then the methods originally tested.
149

H YV ÄRINEN AND S MITH

Proposed measures for linear acyclic model (LiNGAM)
Assumptions on non-Gaussianity
None
SparseEquation for main new pairwise measure(2) with (3) (5)Equation for previous measure by Dodge and Rousson ——Skewed(17)(14)

Proposed measures for extensions of LiNGAM
Extension type
Cyclic case
Nonlinear case
Equation for new pairwise measure Any LiNGAM measure (23) with (3); or (24)
Table 3: The pairwise measures proposed in this paper recapitulated. The new cumulant-based
approximations (Equations 6 and 10) have been omitted since they are mainly for theoretical analysis and not for practical use. Some of the measures by Dodge and Rousson(2001); Dodge and Yadegari (2010) would otherwise fit the “sparse” category but they assume the disturbances to be Gaussian and are thus less general. In the cyclic case, no new
pairwise measures were introduced and it was merely proposed that the LiNGAM measures can be directly used even in the cyclic case. Of the measures highlighted above, thesparse-LiNGAM measure in Equation (5) was proposed in an earlier report on this work
(Hyvärinen, 2010) while all others are new.

Thus, when estimating the LiNGAM model, it may be important to choose a suitable algorithm
depending on data dimension, sample size, noise level, the distributions of the external influences,and other relevant factors.
Basic code for the pairwise measures is distributed on the Internet.10

Acknowledgments
We are grateful to Shohei Shimizu and Patrik Hoyer for deep and insightful comments on the
manuscript, as well as to Christian Beckmann and Mark Woolrich for interesting discussions. We
would also like to thank an anonymous referee for improving the derivation in Section 5.3. This
work was supported by Academy of Finland, Computational Science Program and the Finnish
Centre-of-Excellence in Algorithmic Data Analysis.

References
R.B. Buxton, E.C. Wong, and L.R. Frank. Dynamics of blood flow and oxygenation changes during
brain activati on: the balloon model. Magnetic Resonance in Medicine, 39:855–864, 1998.
C. Chang, M.E. Thomason, and G.H. Glover. Mapping and correction of vascular hemodynamic
latency in the BOLD signal. NeuroImage, 43:90–102, 2008.
10. Code can be found at http://www.cs.helsinki.fi/u/ahyvarin/code/pwcausal/.

150

PAIRWISE L IKELIHOOD R ATIOS FOR N ON -G AUSSIAN SEM S

Z. Chen and L. Chan. Causal discovery for linear non-gaussian acyclic models in the presence of latent gaussian confounders. In Proc. Int. Conf. on Latent Variable Analysis and Signal Separation,
pages 17–24, 2012.
P. Comon. Independent component analysis—a new concept? Signal Processing, 36:287–314,
1994.
P. Daniušis, D. Janzing, J. Mooij, J. Zscheischler, B. Steudel, K. Zhang, and B. Schölkopf. Inferring
deterministic causal relations. In Proc. 26th Conference on Uncertainty in Artificial Intelligence(UAI2010), 2010.
Y. Dodge and V. Rousson. On asymmetric properties of the correlation coefficient in the regression
setting. The American Statistician, 55:51–54, 2001.
Y. Dodge and I. Yadegari. On direction of dependence. Metrika, 72:139–150, 2010.
K. J. Friston, L. Harrison, and W. Penny. Dynamic causal modelling. NeuroImage, 19(4):1273–
1302, 2003.
A. Gretton, K. Fukumizu, C.-H. Teo, L. Song, B. Schölkopf, and A. Smola. A kernel statistical testof independence. In Advances in Neural Information Processing Systems, volume 20. MIT Press,
2008.
D.A. Handwerker, J.M. Ollinger, and M. D’Esposito. Variation of BOLD hemodynamic responses
across subjects and brain regions and their effects on statistical analyses. NeuroImage, 21:1639–
1651, 2004.
P. O. Hoyer, A. Hyvärinen, R. Scheines, P. Spirtes, J. Ramsey, G. Lacerda, and S. Shimizu. Causal
discovery of linear acyclic models with arbitrary distributions. In Proc. 24th Conf. on Uncertainty
in Artificial Intelligence (UAI2008), pages 282–289, Helsinki, Finland, 2008.
P. O. Hoyer, D. Janzing, J. Mooij, J. Peters, and B. Schölkopf. Nonlinear causal discovery withadditive noise models. In Advances in Neural Information Processing Systems, volume 21, pages
689–696. MIT Press, 2009.
A. Hyvärinen. New approximations of differential entropy for independent component analysis andprojection pursuit. In Advances in Neural Information Processing Systems, volume 10, pages
273–279. MIT Press, 1998.
A. Hyvärinen. Fast and robust fixed-point algorithms for independent component analysis. IEEE
Transactions on Neural Networks, 10(3):626–634, 1999.
A. Hyvärinen. Pairwise measures of causal direction in linear non-gaussian acyclic models. In Proc.
Asian Conf. on Machine Learning, JMLR W&CP, volume 13, pages 1–16, Tokyo, Japan, 2010.
A. Hyvärinen, J. Karhunen, and E. Oja. Independent Component Analysis. Wiley Interscience,
2001.
A. Hyvärinen, K. Zhang, S. Shimizu, and P. O. Hoyer. Estimation of a structural vector autoregression model using non-gaussianity. J. of Machine Learning Research, 11:1709–1731, 2010.
151

H YV ÄRINEN AND S MITH

J. Karvanen and V. Koivunen. Blind separation methods based on pearson system and its extensions.
Signal Processing, 82(4):663–573, 2002.
J. M. Mooij, O. Stegle, D. Janzing, K. Zhang, and B. Schölkopf. Probabilistic latent variable models
for distinguishing between cause and effect. In J. Lafferty, C. K. I. Williams, J. Shawe-Taylor,
R.S. Zemel, and A. Culotta, editors, Advances in Neural Information Processing Systems 23(NIPS*2010), pages 1687–1695, 2010.
R.S. Patel, F.D. Bowman, and J.K. Rilling. A Bayesian approach to determining connectivity of thehuman brain˙ Human Brain Mapping, 27:267–276, 2006.
D.-T. Pham and P. Garrat. Blind separation of mixture of independent sources through a quasimaximum likelihood approach. IEEE Trans. on Signal Processing, 45(7):1712–1725, 1997.
J. D. Ramsey, S. J. Hanson, and C. Glymour. Multi-subject search correctly identifies causal connections and most causal directions in the DCM models of the Smith et al. simulation study.
NeuroImage, 58(3):838–848, 2011.
S. Shimizu, P. O. Hoyer, A. Hyvärinen, and A. Kerminen. A linear non-Gaussian acyclic model forcausal discovery. J. of Machine Learning Research, 7:2003–2030, 2006.
S. Shimizu, A. Hyvärinen, Y. Kawahara, and T. Washio. A direct method for estimating a causal
ordering in a linear non-gaussian acyclic model. In Proc. 25th Conference on Uncertainty in
Artificial Intelligence (UAI2009), pages 506–513, Montréal, Canada, 2009.
S. Shimizu, T. Inazumi, Y. Sogawa, A. Hyvärinen, Y. Kawahara, T. Washio, P. O. Hoyer, andK. Bollen. DirectLiNGAM: A direct method for learning a linear non-gaussian structural equation model. J. of Machine Learning Research, 12:1225–1248, 2011.
S. M. Smith, K. L. Miller, G. Salimi-Khorshidi, M. Webster, C. F. Beckmann, T. E. Nichols, J. D.
Ramsey, and M. W. Woolrich. Network modelling methods for fMRI. NeuroImage, 54:875–891,
2011.
Y. Sogawa, S. Shimizu, Y. Kawahara, and T. Washio. An experimental comparison of linear nongaussian causal discovery methods and their variants. In Proc. Int. Joint Conf. on Neural Networks(IJCNN2010), Barcelona, Spain, 2010.
Y. Sogawa, S. Shimizu, A. Hyvärinen, T. Washio, T. Shimamura, and S. Imoto. Estimating exogenous variables in data with more variables than observations. Neural Networks, 24(8):875–880,
2011.
P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search. Springer-Verlag, 1993.
K. Zhang and A. Hyvärinen. On the identifiability of the post-nonlinear causal model. In Proc.
25th Conference on Uncertainty in Artificial Intelligence (UAI2009), pages 647–655, Montréal,Canada, 2009.

152