LiNGAM (version ICA) à comprendre avec des formules mathématiques et Python

Je suis un débutant dans le raisonnement causal.

Récemment, j'ai eu l'occasion d'apprendre la dernière méthode de recherche causale statistique, mais je me suis rendu compte que si je ne comprenais pas la méthode de base ** LiNGAM **, je ne pouvais pas en parler du tout, j'ai donc utilisé la formule ** LiNGAM **. Je vais essayer de le comprendre en le lisant et en suivant le processus avec Python un par un. Il existe deux approches pour estimer ** LiNGAM **, mais cette fois, nous utiliserons l'approche d'analyse en composantes indépendantes (ICA).

Pour la mise en œuvre, un livre sur l'exploration causale statistique (série professionnelle d'apprentissage automatique) [1], un article de LiNGAM [2] et code du Github original ) A été mentionnée.

Ce que vous pouvez faire avec LiNGAM

Par exemple, supposons que vous ayez des données d'observation $ N $ avec quatre variables d'observation $ x_1 $, $ x_2 $, $ x_3 $, $ x_4 $. À partir de ces données, nous dérivons un graphe causal linéaire, non circulaire et dont les erreurs suivent une distribution non gaussienne (voir figure ci-dessous).

Modèle LiNGAM

La formule du modèle est la suivante. Pour les variables observées $ p $ $ x_1, x_2, x_3,…, x_p , le modèle LiNGAM $x_i = \displaystyle \sum_{j≠i} b_{i,j} x_j + e_i \qquad i=1,2,...,p\tag{1}$$ Cependant, $ e_i $ est une variable d'erreur et suppose l'indépendance avec une distribution continue non gaussienne. La représentation matricielle du modèle ressemble à ceci: ${\bf x} = {\bf Bx} + {\bf e} \tag{2}$ $ {\ bf B} $ est une matrice de coefficients, une matrice carrée de $ p × p $. Puisque nous supposons la non-circulation pour le graphe causal, si l'ordre des variables observées est réarrangé dans le bon ordre, la matrice de coefficients $ {\ bf B} $ est ** une matrice triangulaire inférieure dans laquelle toutes les composantes diagonales sont 0 (strictes inférieures). Matrice triangulaire) **.

STEP0 Génération de données artificielles

Cette fois, je vais créer les données artificielles suivantes et trouver un graphe causal. Créez des données basées sur le graphique causal illustré dans la figure ci-dessus.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import FastICA
from sklearn.linear_model import LassoLarsIC, LinearRegression
from scipy.optimize import linear_sum_assignment
from graphviz import Digraph

#STEP0 Génération de données artificielles
num_of_data = 1000
num_of_features = 4

_B_true = np.asarray([[0,0,0,0], [2,0,0,0], [1,3,0,0], [4,5,3,0]])

e1 = np.random.uniform(size=num_of_data)
e2 = np.random.uniform(size=num_of_data)
e3 = np.random.uniform(size=num_of_data)
e4 = np.random.uniform(size=num_of_data)

x1 = e1
x2 = _B_true[1,0] * x1 + e2
x3 = _B_true[2,0] * x1 + _B_true[2,1] * x2 + e3
x4 = _B_true[3,0] * x1 + _B_true[3,1] * x2 + _B_true[3,2] * x3 + e4

_X = np.asarray([x1, x2, x3, x4]).T
_columns = ["x1", "x2", "x3", "x4"]

#Trier les variables observées en texte
order_shuffled = [3,0,2,1]
X = _X[:,order_shuffled]
columns = [_columns[i] for i in order_shuffled]
P_os = np.zeros_like(_B_true)
for i,j in enumerate(order_shuffled):
    P_os[i,j] = 1
B_true = P_os@_B_true@P_os.T

DF_X = pd.DataFrame(X,columns=columns)

Après avoir créé la matrice triangulaire inférieure exacte B pour être la solution, réorganisez l'ordre des variables observées en texto. Ici, la commande est $ (x_4, x_1, x_3, x_2) $. En effet, l'ordre correct des variables observées n'est pas connu dans le monde réel.

ÉTAPE1 Estimer la matrice de restauration W

En transformant l'expression $ (2) $, on obtient:

\begin{align}
{\bf x} &= {\bf (I-B)}^{-1} {\bf e}  \\
&= {\bf A} e \tag{3}
\end{align}

Ici, puisque le vecteur de variable d'erreur $ \ bf e $ est indépendant et non gaussien, cette équation peut être considérée comme un modèle ICA. Cela devient {W = IB} $. Cependant, ce que l'ICA trouve peut différer de l'original $ \ bf W $ dans l'ordre et l'échelle des lignes. Par conséquent, la matrice de reconstruction $ {\ bf W_ {ICA}} $ estimée par ICA en utilisant la matrice de substitution $ {\ bf P} $ dans l'ordre des lignes et la matrice diagonale $ \ bf D $ indiquant l'échelle est

\begin{align}
{\bf W_{ICA}} &={\bf PDW}  \\
&= {\bf PD(I-B)}  \tag{4}
\end{align}

Ce sera. Maintenant, demandons en fait $ {\ bf W_ {ICA}} $ en python. Cependant, comme ICA n'est pas le thème cette fois, je vais le demander rapidement avec Fast ICA dans scikit-learn.

#STEP1 W_ICA(PDW)Estimation
ICA = FastICA(random_state=0)
ICA.fit(X)
W_ICA = ICA.components_

ÉTAPE2 Estimer la matrice de substitution P

En transformant l'expression $ (4) $, on obtient:

\begin{align}
{\bf P^{-1}W_{ICA}} &={\bf D(I-B)} \tag{5}
\end{align}

Ici, puisque $ \ bf B $ a une composante diagonale de $ 0 $ due à la non-circulation, la composante diagonale de $ \ bf {IB} $ est $ 1 $, et la composante diagonale de la matrice d'échelle $ \ bf {D} $. N'est pas $ 0 $, donc la composante diagonale de $ \ bf {D (IB)} $ n'est pas $ 0 $. Par conséquent, $ \ bf W_ {ICA} $ doit être remplacé sur le côté gauche de l'équation $ (5) $ pour que $ 0 $ ne vienne pas à la composante diagonale. Par conséquent, si une telle matrice de substitution $ \ bf \ tilde {P} $ est utilisée, la matrice d'estimation est

\begin{align}
{\bf \hat{\tilde{P}} = \mathop{\rm arg~min}\limits_{\bf \tilde{P}} \sum_{i=1}^p \frac{1}{|(\bf {\tilde{P} \hat{W}_{ICA}})_{ii}| }}  \tag{6}
\end{align}

Estimé par. Puisque le problème de trouver une telle matrice de substitution est considéré comme un problème d'allocation linéaire, il peut être résolu par la méthode hongroise ou similaire. Eh bien, c'est une implémentation, mais comme d'habitude, la méthode hongroise n'est pas le thème, donc elle sera résolue rapidement avec scipy.

#STEP2 P_tilde_Estimation du chapeau
epsilon = 1e-16

cost = 1 / np.abs(W_ICA +  epsilon)
row_ind, col_ind = linear_sum_assignment(cost)

P_tilde_hat = np.zeros_like(W_ICA)
for i,j in  zip(row_ind,col_ind):
    P_tilde_hat[i,j] = 1
P_tilde_hat = np.linalg.inv(P_tilde_hat)

Si 0 apparaît lors de la prise du nombre inverse, il explosera, donc une valeur minute ε est ajoutée.

ÉTAPE 3 Matrice d'échelle d'estimation D

Comme d'habitude, la composante diagonale de $ \ bf {I-B} $ est $ 1 $, donc la composante diagonale de $ {\ bf D (I-B)} $ est la même que la matrice diagonale $ \ bf D $. Donc, à partir de l'équation $ (5) $, on peut dire que la matrice diagonale $ \ bf D $ est égale à la composante diagonale de $ {\ bf P ^ {-1} W_ {ICA}} $. Par conséquent, la matrice d'estimation de $ \ bf {D} $ est

\begin{align}
{\bf{\hat D} = diag(\hat{\tilde{P}} \hat{W}_{ICA})} \tag{7}
\end{align}

Ecrit en python, ça ressemble à ça.

#STEP3 D_Estimation du chapeau
D_hat = np.diag(np.diag(P_tilde_hat@W_ICA))

ÉTAPE 4 Estimer la matrice de reconstruction W et la matrice de coefficients B

Estimer la matrice de restauration $ \ bf {W} $ à partir de l'équation $ (5) $

\begin{align}
{\bf \hat W = \hat D^{-1} \hat{\tilde{P}} \hat{W}_{ICA}}  \tag{8}
\end{align}

Et la matrice d'estimation est la matrice de coefficients $ \ bf {B} $

\begin{align}
{\bf \hat B = I - \hat W}  \tag{9}
\end{align}

Ce sera.

#STEP4 W_hat,B_Estimation du chapeau
W_hat = np.linalg.inv(D_hat)@P_tilde_hat@W_ICA
I = np.eye(num_of_features)
B_hat = I - W_hat

ÉTAPE 5 Deviner l'ordre causal

Au stade jusqu'à présent, la formule du modèle LiNGAM elle-même représentée par la formule $ (2) $ est presque obtenue, mais la direction causale elle-même pour dessiner le graphique est inconnue. Par conséquent, nous devons trouver une substitution qui fait de $ \ bf {B} $ une matrice triangulaire inférieure stricte en permutant l'ordre des variables observées dans l'équation $ (2) $. Autrement dit, multipliez la matrice de substitution $ \ bf {\ ddot {P}} $ dans l'équation $ (2) $.

$ \ Bf \ ddot {P} B \ ddot {P} ^ {\ mathrm {T}} $ à ce moment est aussi proche que possible de la matrice triangulaire inférieure exacte $ \ bf {\ ddot {P}} Recherchez $.

Vous pouvez trouver une telle matrice de substitution par recherche forcée, mais le nombre de matrices de substitution possibles est $ p! $, Il y a donc une limite. Par conséquent, la référence [3] propose une méthode pour obtenir cela à grande vitesse.

Dans cette méthode, remplacez d'abord $ p (p + 1) / 2 $ composants par 0 à partir de $ \ bf \ hat {B} $ dans l'ordre croissant de valeur absolue. Après cela, les variables observées sont réorganisées pour déterminer si elles deviennent ou non une matrice triangulaire inférieure, et sinon, la composante suivante plus petite de $ \ bf \ hat {B} $ est remplacée par 0 et la détermination est faite à nouveau.

La méthode du jugement de la matrice triangulaire inférieure (= deviner l'ordre causal), (1) Tout d'abord, laissez la ligne de $ \ bf \ hat {B} $ dont les composantes sont toutes 0 être la ligne $ i $. (S'il n'existe pas, le jugement est Flase) ② Ajoutez $ i $ à la liste d'ordre causal ③ Supprimez la ligne $ i $ et la colonne $ i $ de $ \ bf \ hat {B} $ et revenez à ① comme une nouvelle matrice $ \ bf \ hat {B} $.

La liste d'ordre causal complétée est ancêtre → petit-enfant. Maintenant pour l'implémentation python. Honnêtement, cette partie était trop compliquée à traiter, j'ai donc utilisé le code d'origine comme référence ...

#ÉTAPE 5 Ordre causal(causal order)Devine

pos_list = np.vstack(np.unravel_index(np.argsort(np.abs(B_hat), axis=None), B_hat.shape)).T
initial_zero_num = num_of_features * (num_of_features + 1) // 2

for i, j in pos_list[:initial_zero_num]:    #p(p+1)/Remplacer par 2 0
    B_hat[i, j] = 0

Mtx_lowtri = B_hat
causal_order = []
original_index = np.arange(num_of_features)

for i in range(num_of_features):    #Jugement strict de la matrice triangulaire inférieure
    idx_row_all_zero_list = np.where(np.sum(np.abs(Mtx_lowtri), axis=1) == 0)[0]
    if len(idx_row_all_zero_list) == 0:
        print("Ne devient pas une matrice triangulaire inférieure stricte")
        break
    idx_row_all_zero = idx_row_all_zero_list[0]
    causal_order.append(original_index[idx_row_all_zero])
    original_index = np.delete(original_index, idx_row_all_zero, axis=0)

    mask = np.delete(np.arange(len(Mtx_lowtri)), idx_row_all_zero, axis=0)
    Mtx_lowtri = Mtx_lowtri[mask][:, mask]

P_doubledot = np.zeros_like(B_hat)
for i,j in enumerate(causal_order):
    P_doubledot[j,i] = 1

Dans ces données, nous avons trouvé un ordre causal qui devient une matrice triangulaire inférieure avec un seul tir, de sorte que le code ne boucle pas. Voir l'original pour une mise en œuvre correcte.

ÉTAPE 6 Finition et élagage

Si le nombre de données est suffisant, la matrice de coefficients dérivée dans STEP 5 est correcte, mais sinon, il semble que le modèle puisse être amélioré en utilisant la taille comme finition finale.

Le livre propose d'utiliser le lasso adaptatif comme l'une des méthodes d'élagage. Adaptive Lasso rend la sélection de variable cohérente en pondérant le terme de régularisation Lasso avec $ \ bf w $. Il a été proposé que ce $ \ bf w $ utilise l'inverse des coefficients estimés par régression linéaire. Le lasso adaptatif est effectué pour chaque variable observée selon l'ordre causal obtenu dans STEP5, et la solution finale est obtenue.

Dans la mise en œuvre, le coefficient est d'abord estimé en effectuant une régression linéaire, le poids est obtenu à partir du coefficient estimé, le poids est appliqué aux données d'entrée de Lasso et la quantité estimée qui est la sortie de Lasso est pondérée.

#ÉTAPE 6 Élagage
B_pred = np.zeros_like(B_hat, dtype='float64')
lr = LinearRegression()
reg = LassoLarsIC(criterion='bic')

for i in range(1,num_of_features):
    #adaptive lasso
    predictors = causal_order[:i]    #Variables prédictives
    target = causal_order[i]    #Variable cible
    lr.fit(X[:, predictors], X[:, target])

    weight = 1/(np.abs(lr.coef_) + epsilon)
    reg = LassoLarsIC(criterion='bic')
    reg.fit(X[:, predictors] * weight, X[:, target])

    B_pred[causal_order[i], causal_order[:i]] = reg.coef_ * weight

Visualisation

Enfin, le modèle obtenu est visualisé comme le montre la figure ci-dessous. De gauche, le vrai graphique, le graphique obtenu dans STEP5 et le graphique traité dans STEP6. Cette fois, les données sont simples et nous avons préparé suffisamment de données, donc cela fonctionne bien comme de STEP5.

Graphviz a été utilisé pour la visualisation. J'ai utilisé le code de création de graphique suivant.

def Make_Graph(columns, B):
    Graph = Digraph(format="png")
    for i, j in enumerate(columns):
        for k, l in enumerate(columns):
            if B[i, k] != 0:
                Graph.edge(l, j, color="black", label= "    " + str(B[i, k].round(3)))
    return Graph

en conclusion

Comme mentionné ci-dessus, un débutant dans le raisonnement causal a comparé des livres [1] et des articles [2], a examiné minutieusement le code de la famille principale et a travaillé dur pour lire LiNGAM.

Les références

** [1] ** Shohei Shimizu. Recherche causale statistique (série professionnelle d'apprentissage automatique). Kodansha. [2] Shimizu et al. A linear non-Gaussian acyclic model for causal discovery. Journal of Machine Learning Research 7.Oct (2006): 2003-2030. [3] Hoyer et al. New permutation algorithms for causal discovery using ICA. In Proceedings of International Conference on Independent Component Analysis and Blind Signal Separation(2006).