## A Bayesian estimation approach to analyze non-Gaussian data-generating processes with latent classes

Naoki Tanaka ∗

Shohei Shimizu ∗

Takashi Washio ∗

## Abstract

A large amount of observational data has been accumulated in various fields in recent times, and there is a growing need to estimate the generating processes of these data. A linear non-Gaussian acyclic model (LiNGAM) based on the non-Gaussianity of external influences has been proposed to estimate the data-generating processes of variables. However, the results of the estimation can be biased if there are latent classes. In this paper, we first review LiNGAM, its extended model, as well as the estimation procedure for LiNGAM in a Bayesian framework. We then propose a new Bayesian estimation procedure that solves the problem.

## 1 Introduction

Several methods have recently been proposed to discover a complete causal structure, that is, all the causal directions, under the assumption that disturbance variables have non-Gaussian distributions. However, the estimation results can be biased if there are 'latent classes.' Latent classes are unobserved discrete variables that have more than one observed child variables. Data that has been generated from different processes are mixed in the presence of latent classes. Several methods have been proposed to estimate the causal structure in the presence of latent classes [12], but all of these are affected by local optima. Therefore, in this paper, we propose a new estimation approach that can solve this problem.

The structure of this paper is as follows. In Section 2, we briefly review the data generating process for estimating causal structure, (LiNGAM, short for Linear Non-Gaussian Acyclic Model) [11], the LiNGAM mixture model [12], existing estimation approaches to the LiNGAM mixture model [12], and the BayesLiNGAM model. In Section 3, we extend previous research dedicated to solving the problem, and test the performance of our new method through experiments with simulated data in Section 4. Section 5 concludes the paper.

∗ The Institute of Scientific and Industrial Research, Osaka University

## 2 Background

## 2.1 Linear non-Gaussian acyclic model (LiNGAM)

We begin by introducing the basic LiNGAM model [11]. LiNGAM is a causal model with the following four assumptions.

- 1. The relations between the observed variables x i ( i = 1 , . . . , n ) can be represented by a directed acyclic graph (DAG), as shown in Fig. 1.
- 2. x i is assigned a value by a linear function of the values already assigned to the variables constituting its parents in the DAG as well as a 'disturbance' (noise) term e i and an optimal constant term µ i :

$$x _ { i } = \sum _ { k ( j ) < k ( i ) } b _ { i j } ( x _ { j } - \mu _ { j } ) + e _ { i } + \mu _ { i } \quad \quad ( 1 )$$

where k i ( ) is a causal order of x i in the DAG (so that if there is a directed edge from x j to x i in the DAG, k j ( ) &lt; k i ( )), b ij represents the strength of the connection between x j and x i in the DAG.

- 3. The e i are all continuous random variables that follow non-Gaussian distributions with zero means and non-zero variances, and the e i are mutually independent, i.e., p e , . . . , e ( 1 n ) = ∏ i p i ( e i ).
- 4. The dataset D = { x 1 , . . . , x N } (each x contains components x i ) are observed, and each data vector x is generated according to the process described above with the same DAG, coefficients b ij , constants µ i , and disturbances e i , sampled independently from the same distributions. Note that this assumption implies that there are either no unobserved (latent) confounders [10] (hidden variables) or that the x i are causally sufficient [13].

Figure 1: An example of LiNGAM described by a DAG.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1408.0337v1-figure-9c3730bb-7321-4362-af0e-370495649c48.jpg)

## 2.2 LiNGAM mixture model

The LiNGAM mixture model [12] is an extension of the basic LiNGAM to cases where the final assumption of the basic LiNGAM model collapses because of the presence of latent classes (hidden groups). If they exist, each latent class has a different structure. 'Structure' here consists of the causal order k i ( ), the connection strengths b ij , the densities of disturbances p i , and the means of observed variables µ i . Therefore, structure in a class c can be described using k ( c ) ( ), i b ( c ) ij , p ( c ) i , and µ ( c ) i (Note that the disturbances e ( c ) i are generated from p ( c ) i ). The data within each class ( c = 1 , . . . , l ) are assumed to be generated by the basic LiNGAM. Thus, the data-generating model within class c can be described as follows:

$$x _ { i } = \sum _ { k ^ { ( c ) } ( j ) < k ^ { ( c ) } ( i ) } b _ { i j } ^ { ( c ) } ( x _ { j } - \mu _ { j } ^ { ( c ) } ) + \mu _ { i } ^ { ( c ) } + e _ { i } ^ { ( c ) }$$

Moreover, the data vectors x are assumed to be generated by the following mixture density:

$$p ( \mathbf x | \Theta ) = \sum _ { c = 1 } ^ { l } \left \{ \begin{matrix} p ( \mathbf x | \theta ^ { ( c ) } ) \\ \theta ^ { ( c ) } \\ \end{matrix} ) \right \} P ( c ) \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \ \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \	 & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad &\quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad & \quad \\ \quad \quad \quad$$

where l is the number of classes, and Θ = [ θ (1) , . . . , θ ( ) l ], and θ ( c ) contain all the parameters in (2), i.e., k ( c ) ( j ) , k ( c ) ( ) i , b ( c ) ij , µ ( c ) j , µ ( c ) i , and p ( c ) i ( e ( c ) i ). Note that if only one latent class exists, the LiNGAM mixture model is equivalent to basic LiNGAM.

## 2.3 Existing estimation approach to LiNGAM mixture model

In this subsection, we briefly introduce the existing estimation approach to the LiNGAM mixture model. In this approach, the LiNGAM mixture model is transformed into an Independent Component Analysis (ICA) [6]) model, and an ICA algorithm (for example, [7] [9]) is applied to it [12].

The LiNGAM mixture model (Equation (2)) can be written in matrix form as follows:

$$x = B ^ { ( c ) } x + ( \mathbf I - B ^ { ( c ) } ) \mu ^ { ( c ) } + \mathbf e ^ { ( c ) } \quad \quad \quad ( 4 ) \\ \mathbf r _ { \lambda } \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad. \quad.$$

where B ( c ) , µ ( c ) , and e ( c ) collect b ( c ) ij , µ ( c ) i , and e ( c ) i , respectively, and I is the identity matrix. We can then obtain the following equation by solving Equation (4) for x :

$$\mathbf x = \mu ^ { ( c ) } + \mathbf A ^ { ( c ) } \mathbf e ^ { ( c ) }$$

where A ( c ) = ( I -B ( c ) ) -1 . We can estimate B ( c ) by calculating A ( c ) using the ICA algorithm and computing B ( c ) . However, there are two important indeterminacies that ICA cannot solve: the order and the scaling of the independent components. In order to solve them, we need to permute and normalize A ( c ) appropriately before using it to compute B ( c ) (See [11] for more details).

## 2.4 The BayesLiNGAM method

Existing estimation methods for the LiNGAM mixture model [7] [9] are affected by local optima, because of which we propose a Bayesian approach. In this subsection, we summarize the BayesLiNGAM method [5], which is an estimation method of basic LiNGAM.

In Bayesian inference, the inference with the highest posterior probability for some hypotheses is selected. In BayesLiNGAM, the hypotheses are different possible DAGs. BayesLiNGAM outputs the DAG with the highest posterior probability of all possible DAGs given the data. The posterior probabilities can be calculated by Bayes' theorem:

$$P ( G _ { m } | \mathbf D ) = \frac { p ( \mathbf D | G _ { m } ) P ( G _ { m } ) } { p ( \mathbf D ) }$$

Here, G m are the different possible DAGs, m = { 1 , . . . , N g } , where N g denotes the number of different DAGs on n variables, and D is the observed dataset D = { x 1 , . . . , x N } , N : sample size). Due to our assumption of independent and identically distributed (IID) data, p ( D ) = ∏ N s =1 p ( x s ). We then need to specify the likelihoods p ( D | G m ), prior probabilities P G ( m ), and normalization constant p ( D ) to compute the posterior probabilities.

The prior probabilities P G ( m ) first incorporate any domain of knowledge and prior information. If we have no knowledge or information about the G m , all values of P G ( m ) are equal, i.e., P G ( m ) = 1 N g .

Furthermore, p ( D ) is a constant that simply normalizes posterior probabilities. Hence, we can obtain p ( D ) as follows:

$$P ( \mathbf D ) = \sum _ { m = 1 } ^ { N _ { g } } p ( \mathbf D | G _ { m } ) P ( G _ { m } )$$

The likelihoods p ( D | G m ) can be computed by marginalizing over p ( DΘ | , G m ) of Θ :

$$P ( \mathbf D | G _ { m } ) = \int p ( \mathbf D | \Theta, G _ { m } ) p ( \Theta | G _ { m } ) d \Theta$$

Here, Θ collects all the parameters in Equation (1) (i.e., k j , k ( ) ( ) i , b ij , µ j , µ i , and p i ), p ( DΘ | , G m ) denotes the likelihood of the model (here, basic LiNGAM), and p ( Θ | G m ) denotes the prior distributions of Θ . Note that if the G m are given, we can find the parents of x i in the G m and determine the causal order k j ( ) and k i ( ). That is, the first term of the right side in Equation (1) can be specified by the G m . We introduce the likelihood of the model and prior distribution for Θ in Section 3.

## 3 Proposed method

## 3.1 Model

We assume that the data within each class ( c = 1 , . . . , l ) are generated by the basic LiNGAM. Thus, the data-generating model in class c can be described by Equation (2), and the probability densities of the data x are the same as in Equation (3). We then use the Bayesian approach to estimate the LiNGAM mixture model in the same manner as BayesLiNGAM in Section 2.4.

Here, we need to parameterize the densities p i to specify p ( DΘ | , G m ) in Equation (8). Due to the assumption in LiNGAM that the e i are all continuous random variables with non-Gaussian distributions, we use a generalized Gaussian distribution [8] that incorporates a shape parameter. The generalized Gaussian distribution is symmetric and includes Gaussian, Laplacian, continuous uniform, and several non-Gaussian distributions. The probability density function is as follows:

$$p _ { i } ( e _ { i } ) = \frac { \lambda _ { i } \sqrt { \frac { \Gamma ( 3 / \lambda _ { i } ) } { \Gamma ( 1 / \lambda _ { i } ) } } \exp ( - ( \sqrt { \frac { \Gamma ( 3 / \lambda _ { i } ) } { \Gamma ( 1 / \lambda _ { i } ) } } \frac { | e _ { i } | } { \sigma _ { i } } ) ^ { \lambda _ { i } } ) } { 2 \sigma _ { i } \Gamma ( 1 / \lambda _ { i } ) }$$

Here, σ i represent the standard deviations, λ i are the shape parameters, and Γ() denotes the Gamma function.

Moreover, we need to transform the density of the e i into that of the x i . The LiNGAM mixture model (Equation (2)) can be written in matrix form as follows:

$$\mathbf x \ & = \ \mathbf f ( \mathbf e ^ { ( c ) } ) \\ & = \ \mu ^ { ( c ) } + ( \mathbf I - \mathbf B ^ { ( c ) } ) ^ { - 1 } \mathbf e ^ { ( c ) } \\ & \quad \dots \ \dots \ \dots \quad \mathbf w \ \nu$$

where f () is a mapping vector, and B ( c ) , e ( c ) , and µ ( c ) collect b ( c ) ij , e ( c ) i , and µ ( c ) i , respectively. I is the identity matrix.

$$p ( x | \theta ^ { ( c ) }, G _ { m } ) & = \frac { 1 } { | \det J f ( \mathbf f ^ { - 1 } ( x ) ) | } p _ { \mathbf e } ^ { ( c ) } ( \mathbf f ^ { - 1 } ( x ) ) \quad \quad ( 1 1 ) \\ \tau \cdot \quad \tau \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdots$$

The density of x within class c ( p ( x | θ ( c ) , G m ) ) is obtained from the density of e ( c ) ( p ( c ) e ( e ( c ) ) ) as follows [6]:

where J f is a Jacobian matrix. From equations (10) and (11) with the assumption of acyclicity, det J f ( f -1 ( x )) equals to one, and we can obtain the following equation:

$$p ( x | \theta ^ { ( c ) }, G _ { m } ) = p _ { \mathbf e } ^ { ( c ) } ( \mathbf e ^ { ( c ) } )$$

According to the assumption in LiNGAM that the e i are mutually independent,

$$p _ { \text{e} } ^ { ( c ) } ( \text{e} ) \ = \ \prod _ { i = 1 } ^ { n } p _ { i } ^ { ( c ) } ( e _ { i } ^ { ( c ) } ).$$

Then, the probability density of x within each class c , using equations (2), (12), and (13), is given by

$$. 3 ), & \text{ is given by} \\ & p ( x | \theta ^ { ( c ) }, G _ { m } ) \\ & = \prod _ { i = 1 } ^ { n } p _ { i } ^ { ( c ) } ( e _ { i } ^ { ( c ) } ) \\ & = \prod _ { i = 1 } ^ { n } p _ { i } ^ { ( c ) } ( x _ { i } - \mu _ { i } ^ { ( c ) } - \sum _ { k ^ { ( c ) } ( j ) < k ^ { ( c ) } ( i ) } b _ { i j } ^ { ( c ) } ( x _ { j } - \mu _ { j } ^ { ( c ) } ) ). \\ \ n \, \text{ now specify} \, p ( \text{D} | \Theta. G _ { m } ) \, \text{ with the assumption of IID data and equations}$$

We can now specify p ( DΘ | , G m ) with the assumption of IID data and equations (3) and (14):

$$\cos \liminf _ { s = 1 } ^ { N } p ( x ^ { s } | \Theta, G _ { m } ) & & ( 1 5 ) \\ p ( x | \Theta, G _ { m } ) & & ( 1 5 ) \\ = \sum _ { c = 1 } ^ { l } \left \{ p ( x | \theta ^ { ( c ) }, G _ { m } ) \right \} P ( c ) \\ = \sum _ { c = 1 } ^ { l } \left \{ \prod _ { i = 1 } ^ { n } p _ { i } ^ { ( c ) } ( x _ { i } - \mu _ { i } ^ { ( c ) } - \sum _ { k ^ { ( c ) } ( j ) < k ^ { ( c ) } ( i ) } b _ { i j } ^ { ( c ) } ( x _ { j } - \mu _ { j } ^ { ( c ) } ) ) \right \} \\ \times P ( c ) & & ( 1 6 ) \\ 2 \ \Pr i o r \ \text{ parameters and hyperparameters}$$

## 3.2 Priors for parameters and hyperparameters

We use multinomial distribution for P c ( ) in Equation (3). The distribution represents, for N independent trials, each of which leads to one of the possible events c = 1 , . . . , l , the number of times that each event occurs ( z = [ z (1) , . . . , z ( ) l ]). w ( c ) is the probability of the occurrence of each event. Generally, the probability mass function of a multinomial distribution is as follows:

$$P ( \mathbf z ) = \frac { N! } { \prod _ { c = 1 } ^ { l } z ^ { ( c ) }! } \prod _ { c = 1 } ^ { l } ( w ^ { ( c ) } ) ^ { z ^ { ( c ) } } & & ( 1 7 ) \\ > 0 \text{ and } \sum _ { \varepsilon = 1 } ^ { l } w ^ { ( c ) } = 1. \text{ When we use the above distribution in}$$

where w ( c ) &gt; 0 and ∑ l c =1 w ( c ) = 1. When we use the above distribution in Equation (3), the number of trials N is one. Then, z can be regarded as the indicator vector and corresponds to variable c . For example, if x belongs to class 2 (i.e., c = 2), z (2) = 1 and the other components are 0. From the above, we can use a multinomial distribution for P c ( ) in the following form:

$$P ( c ) & = w ^ { ( c ) }$$

We use a Dirichlet distribution for the parameters of the multinomial distribution ( w (1) , . . . , w ( ) l ) because the former is a conjugate prior for the latter

and is typically used in similar contexts. The density function of the Dirichlet distribution is as follows:

$$p ( \mathbf w ) = \frac { \Gamma ( \sum _ { c = 1 } ^ { l } a ^ { ( c ) } ) } { \prod _ { c = 1 } ^ { l } \Gamma ( a ^ { ( c ) } ) } \prod _ { c = 1 } ^ { l } ( w ^ { ( c ) } ) ^ { a ^ { ( c ) } - 1 } & & ( 1 9 ) \\ \intertext { r } \, = \, w ^ { ( 1 ) } \dots \, w ^ { ( l ) } \text{ and } a ^ { ( c ) } ( > \, 0 ) \text{ are concentration parameters. } \text{ We }$$

where w = w (1) , . . . , w ( ) l and a ( c ) ( &gt; 0) are concentration parameters. We can generate a random vector w by normalizing independent gamma random variables γ (1) , . . . , γ ( ) l with shape parameters a ( c ) (concentration parameters in the Dirichlet distribution) and scale parameter 1 [4]:

$$w ^ { ( c ) } \ & = \ \frac { \gamma ^ { ( c ) } } { \sum _ { c = 1 } ^ { l } \gamma ^ { ( c ) } } \\ p ( \gamma ^ { ( c ) } ) \ & = \ \frac { 1 } { \Gamma / \wedge ( c ) \wedge } ( \gamma ^ { ( c ) } ) ^ { a ^ { ( c ) } - 1 } \exp ( \gamma ^ { ( c ) } )$$

$$p ( \gamma ^ { ( c ) } ) \ = \ \frac { \sum _ { 1 } ^ { - \cdot - \cdot } } { \Gamma ( a ^ { ( c ) } ) } ( \gamma ^ { ( c ) } ) ^ { a ^ { ( c ) } - 1 } \exp ( \gamma ^ { ( c ) } )$$

We use Gaussian distributions for b ( c ) ij and µ ( c ) i with zero and ϕ ( c ) as their means and v 2 and τ 2 as their variances respectively.

$$p ( b _ { i j } ^ { ( c ) } ) \ = \ \frac { 1 } { \sqrt { 2 \pi v ^ { 2 } } } \exp \left \{ - \frac { ( b _ { i j } ^ { ( c ) } ) ^ { 2 } } { 2 v ^ { 2 } } \right \} \quad \ \quad \ \ ( 2 2 )$$

$$p ( \mu _ { i } ^ { ( c ) } ) \ = \ \frac { 1 } { \sqrt { 2 \pi \tau ^ { 2 } } } \exp \left \{ - \frac { ( \mu _ { i } ^ { ( c ) } - \varphi ^ { ( c ) } ) ^ { 2 } } { 2 \tau ^ { 2 } } \right \}$$

We heuristically determine the value of ϕ . We use the result of the Gaussian mixture model estimation and an expectation-maximization (EM) algorithm [2] for ϕ because the greater the number of latent classes, the more difficult it is to estimate µ ( c ) i . The EM algorithm is a method to compute maximum likelihood solutions for models with latent variables. In the algorithm, parameters and responsibilities (conditional probabilities of latent variables given the data) are updated in turn until the change in the likelihood function falls below some threshold.

We use inverse gamma distributions [3] with shape parameters α , η , and χ and scale parameters β , ζ , and /epsilon1 for ( σ ( c ) i ) 2 , λ ( c ) i , and v 2 , respectively.

$$p ( ( \sigma _ { i } ^ { ( c ) } ) ^ { 2 } ) \ = \ \frac { \beta ^ { \alpha } } { \Gamma ( \alpha ) } ( \sigma _ { i } ^ { ( c ) } ) ^ { 2 ( - \alpha - 1 ) } \exp \left \{ \frac { - \beta } { ( \sigma _ { i } ^ { ( c ) } ) ^ { 2 } } \right \} \quad \quad ( 2 4 )$$

$$p ( \lambda _ { i } ^ { ( c ) } ) \ = \ \frac { \zeta ^ { \eta } } { \Gamma ( \eta ) } ( \lambda _ { i } ^ { ( c ) } ) ^ { - \eta - 1 } \exp \left ( \frac { - \zeta } { \lambda _ { i } ^ { ( c ) } } \right ) ^ { \mu \times \imath } \, ^ { \prime } \, \cdots \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

$$p ( v ^ { 2 } ) \ = \ \frac { \epsilon ^ { \chi } } { \Gamma ( \chi ) } v ^ { 2 ( - \chi - 1 ) } \exp \left ( \frac { - \epsilon } { v ^ { 2 } } \right )$$

We can generate inverse gamma random variables using gamma random variables. If a variable X has a gamma distribution with shape parameter α and

scale parameter β , Y = 1 /X has an inverse gamma distribution with shape parameter α and scale parameter 1 /β [3].

We determine the value of a ( c ) , α β , , η , ζ , χ /epsilon1 , , and τ arbitrarily. From the above, we can compute P ( D | G m ) in Equation (8) using ordinary Monte Carlo sampling [2] to compute the integral.

## 4 Simulation

In this simulation, we assume for simplicity that we know of the existence of a causal connection between two observed variables, and that its direction is the same for all classes. We can thus estimate which direction is true ( x ( c ) 1 → x ( c ) 2 or x ( c ) 1 ← x ( c ) 2 ). We generated 1,000 datasets under every combination of sample size ( N = 50 100 500) and the number of classes ( , , l = 2 4 6). , ,

The data within each class were generated by the basic LiNGAM (Equation(1)), following which we mixed them. In all datasets, the true model was x ( c ) 1 → x ( c ) 2 and the connection strength was b ( c ) ij = [ -1 5 . , -0 5] . ∩ [0 5 . , 1 5]. . The distributions of external influences e ( c ) i were randomly selected from the following three non-Gaussian distributions with one variance: Laplace, Uniform ([ -√ 3 , √ 3]) and Student-t (with five degrees of freedom).

With regard to the priors we did not use any prior information about G m , and thus P G ( m ) = 1 2 . We randomly selected the value of a ( c ) from among 3, 5, and 7. Our choice of values for the rest of the parameters was α = 3, β = 3, η = 3, ζ = 3, χ = 3, /epsilon1 = 3, and τ = 0 5. .

We selected the number of classes l as follows. We calculated the logmarginal likelihoods of the two models for all numbers of classes and selected the class with the largest log-marginal likelihood. The maximum number of classes tested was 2log N . This was motivated by the fact that when the sample size N approaches infinity in a Dirichlet process mixture model, the number of classes l approaches log N [1].

We compared our method with the existing method [12] to determine the one that can correctly estimate causal directions more times given the above datasets. However, we use [9] in the existing method [12] instead of [7]. We show the simulation results in Table 1. The numbers in the tables denote the number of times each method correctly estimated causal directions. As shown in the table, our method has more correct answers than the existing method.

## 5 Conclusion

In this paper we proposed a new estimation approach for discovering causal structure in the presence of latent classes. In a simulation run using artificial data, our method correctly estimated more causal directions than the existing method. Our plan for future research is to evaluate our method on a wide variety of real datasets.

|                 | Table   |   1: Simulation results Sample size |   1: Simulation results Sample size |     |
|-----------------|---------|-------------------------------------|-------------------------------------|-----|
|                 |         |                                  50 |                                 100 | 500 |
| Our method      | l = 2   |                                 913 |                                 947 | 981 |
|                 | l = 4   |                                 908 |                                 937 | 973 |
|                 | l = 6   |                                 922 |                                 957 | 967 |
| Existing method | l = 2   |                                 649 |                                 657 | 684 |
|                 | l = 4   |                                 663 |                                 655 | 729 |
|                 | l = 6   |                                 646 |                                 700 | 762 |

## References

- [1] C. E. Antoniak. Mixtures of dirichlet processes with applications to bayesian nonparametric problems. The annals of statistics , pages 11521174, 1974.
- [2] C. M. Bishop. Pattern Recognition and Machine Learning . Springer, 2006.
- [3] John D Cook. Inverse gamma distribution. online: h ttp://www. johndcook. com/inverse gamma. pdf, Tech. Rep , 2008.
- [4] L Devroye. Non-uniform Random Variate Generation . Springer-Verlag, 1986.
- [5] P. O. Hoyer and A. Hyttinen. Bayesian discovery of linear acyclic causal models. In Proc. 25th Conf. on Uncertainty in Artificial Intelligence (UAI2009) , pages 240-248, 2009.
- [6] A. Hyv¨ arinen, J. Karhunen, and E. Oja. Independent component analysis . Wiley, New York, 2001.
- [7] M. N. H. Mollah, M. Minami, and S. Eguchi. Exploring latent structure of mixture ICA models by the minimum β -divergence method. Neural Computation , 18:166-190, 2006.
- [8] S. Nadarajah. A generalized normal distribution. Journal of Applied Statistics , 32(7):685-694, 2005.
- [9] J. A. Palmer, S. Makeig, K. Kreutz-Delgado, and B. D. Rao. Newton method for the ICA mixture model. In Proc. IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP2008) , pages 1805-1808, 2008.
- [10] J. Pearl. Causality: Models, Reasoning, and Inference . Cambridge University Press, 2000.
- [11] S. Shimizu, P. O. Hoyer, A. Hyv¨rinen, and A. Kerminen. a A linear nongaussian acyclic model for causal discovery. J. Machine Learning Research , 7:2003-2030, 2006.

- [12] S. Shimizu and A. Hyv¨ arinen. Discovery of linear non-gaussian acyclic models in the presence of latent classes. In Proc. 14th Int. Conf. on Neural Information Processing (ICONIP2007) , pages 752-761, 2008.
- [13] P. Spirtes, C. Glymour, and R. Scheines. Causation, Prediction, and Search . Springer Verlag, 1993. (2nd ed. MIT Press 2000).
Naoki Tanaka Shohei Shimizu Takashi Washio
A Bayesian estimation approach to analyze non-Gaussian data-generating processes with latent classes
A large amount of observational data has been accumulated in various fields
in recent times, and there is a growing need to estimate the generating
processes of these data. A linear non-Gaussian acyclic model (LiNGAM) based on
the non-Gaussianity of external influences has been proposed to estimate the
data-generating processes of variables. However, the results of the estimation
can be biased if there are latent classes. In this paper, we first review
LiNGAM, its extended model, as well as the estimation procedure for LiNGAM in a
Bayesian framework. We then propose a new Bayesian estimation procedure that
solves the problem.