LiNGAM (ICA-Version) mit mathematischen Formeln und Python zu verstehen

Ich bin ein Anfänger im kausalen Denken.

Vor kurzem hatte ich die Gelegenheit, die neueste Methode der statistischen Kausalsuche zu lernen, aber mir wurde klar, dass ich, wenn ich die grundlegende Methode ** LiNGAM ** nicht verstehe, überhaupt nicht darüber sprechen kann, also habe ich die Formel ** LiNGAM ** verwendet. Ich werde versuchen, es zu verstehen, indem ich es lese und den Prozess mit Python nacheinander verfolge. Es gibt zwei Ansätze zur Schätzung von ** LiNGAM **, aber dieses Mal werden wir den Ansatz der unabhängigen Komponentenanalyse (ICA) verwenden.

Zur Implementierung ein Buch über statistische Kausalforschung (Professional Series für maschinelles Lernen) [1], ein Artikel von LiNGAM [2] und Code des ursprünglichen Github ) Wurde verwiesen.

Was Sie mit LiNGAM machen können

Angenommen, Sie haben $ N $ Beobachtungsdaten mit vier Beobachtungsvariablen $ x_1 $, $ x_2 $, $ x_3 $, $ x_4 $. Aus diesen Daten leiten wir einen linearen, nicht kreisförmigen Kausaldiagramm ab, dessen Fehler einer nicht-Gaußschen Verteilung folgen (siehe Abbildung unten).

LiNGAM-Modell

Die Formel für das Modell lautet wie folgt. Für $ p $ beobachtete Variablen $ x_1, x_2, x_3,…, x_p $ das LiNGAM-Modell $x_i = \displaystyle \sum_{j≠i} b_{i,j} x_j + e_i \qquad i=1,2,...,p\tag{1}$ $ E_i $ ist jedoch eine Fehlervariable und setzt Unabhängigkeit mit einer nicht-Gaußschen kontinuierlichen Verteilung voraus. Die Matrixdarstellung des Modells sieht folgendermaßen aus: ${\bf x} = {\bf Bx} + {\bf e} \tag{2}$ $ {\ bf B} $ ist eine Koeffizientenmatrix, eine quadratische Matrix von $ p × p $. Da wir für den Kausaldiagramm eine Nichtzirkulation annehmen, ist die Koeffizientenmatrix $ {\ bf B} $ ** eine untere Dreiecksmatrix, in der alle diagonalen Komponenten 0 sind (streng niedriger), wenn die Reihenfolge der beobachteten Variablen in der richtigen Reihenfolge neu angeordnet wird. Dreiecksmatrix) **.

SCHRITT 0 Erzeugung künstlicher Daten

Dieses Mal werde ich die folgenden künstlichen Daten erstellen und einen Kausaldiagramm finden. Erstellen Sie Daten basierend auf dem in der obigen Abbildung gezeigten Kausaldiagramm.

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

#SCHRITT 0 Künstliche Datengenerierung
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"]

#Sortieren Sie die beobachteten Variablen in Text
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)

Nachdem Sie die exakte untere Dreiecksmatrix B als Lösung erstellt haben, ändern Sie die Reihenfolge der beobachteten Variablen in Texto. Hier ist die Reihenfolge $ (x_4, x_1, x_3, x_2) $. Dies liegt daran, dass die korrekte Reihenfolge der beobachteten Variablen in der realen Welt nicht bekannt ist.

SCHRITT 1 Schätzen Sie die Restaurationsmatrix W.

Wenn wir den Ausdruck $ (2) $ transformieren, erhalten wir:

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

Da hier der Fehlervariablenvektor $ \ bf e $ unabhängig und nicht Gaußsch ist, kann diese Gleichung als ICA-Modell betrachtet werden, und wenn die inverse Matrix der Matrix $ \ bf A $ die Wiederherstellungsmatrix $ \ bf $ ist, dann $ \ bf Es wird {W = IB} $. Was die ICA findet, kann jedoch in Zeilenreihenfolge und Skalierung vom ursprünglichen $ \ bf W $ abweichen. Daher wird die Rekonstruktionsmatrix $ {\ bf W_ {ICA}} $ von ICA unter Verwendung der Substitutionsmatrix $ {\ bf P} $ in der Zeilenreihenfolge und der Diagonalmatrix $ \ bf D $, die die Skala angibt, geschätzt

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

Es wird sein. Lassen Sie uns nun tatsächlich nach $ {\ bf W_ {ICA}} $ in Python fragen. Da ICA diesmal jedoch nicht das Thema ist, werde ich mit Fast ICA in scikit-learn schnell danach fragen.

#STEP1 W_ICA(PDW)Schätzen
ICA = FastICA(random_state=0)
ICA.fit(X)
W_ICA = ICA.components_

SCHRITT 2 Schätzen Sie die Substitutionsmatrix P.

Wenn wir den Ausdruck $ (4) $ transformieren, erhalten wir:

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

Da $ \ bf B $ aufgrund von Nichtzirkulation eine diagonale Komponente von $ 0 $ hat, beträgt die diagonale Komponente von $ \ bf {IB} $ $ 1 $ und die diagonale Komponente der Skalenmatrix $ \ bf {D} $. Ist nicht $ 0 $, also ist die diagonale Komponente von $ \ bf {D (IB)} $ nicht $ 0 $. Daher muss $ \ bf W_ {ICA} $ auf der linken Seite von Gleichung $ (5) $ ersetzt werden, damit $ 0 $ nicht zur diagonalen Komponente kommt. Wenn daher eine solche Substitutionsmatrix $ \ bf \ tilde {P} $ verwendet wird, ist die Schätzmatrix

\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}

Geschätzt von. Da das Problem, eine solche Substitutionsmatrix zu finden, als lineares Zuordnungsproblem angesehen wird, kann es durch das ungarische Verfahren oder dergleichen gelöst werden. Nun, es ist eine Implementierung, aber wie üblich ist die ungarische Methode nicht das Thema, daher wird sie mit scipy schnell gelöst.

#STEP2 P_tilde_Hutschätzung
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)

Wenn beim Aufnehmen der inversen Zahl 0 erscheint, explodiert diese, sodass ein Minutenwert ε addiert wird.

SCHRITT 3 Skalenmatrix D schätzen

Wie üblich ist die diagonale Komponente von $ \ bf {I-B} $ $ 1 $, daher ist die diagonale Komponente von $ {\ bf D (I-B)} $ dieselbe wie die diagonale Matrix $ \ bf D $. Aus der Gleichung $ (5) $ kann daher gesagt werden, dass die Diagonalmatrix $ \ bf D $ gleich der Diagonalkomponente von $ {\ bf P ^ {-1} W_ {ICA}} $ ist. Daher ist die Schätzmatrix von $ \ bf {D} $

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

In Python geschrieben, sieht es so aus.

#STEP3 D_Hutschätzung
D_hat = np.diag(np.diag(P_tilde_hat@W_ICA))

SCHRITT 4 Schätzen Sie die Rekonstruktionsmatrix W und die Koeffizientenmatrix B.

Schätzen Sie die Wiederherstellungsmatrix $ \ bf {W} $ aus der Gleichung $ (5) $

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

Und die Schätzmatrix ist die Koeffizientenmatrix $ \ bf {B} $

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

Es wird sein.

#STEP4 W_hat,B_Hutschätzung
W_hat = np.linalg.inv(D_hat)@P_tilde_hat@W_ICA
I = np.eye(num_of_features)
B_hat = I - W_hat

SCHRITT 5 Erraten der kausalen Ordnung

Zu diesem Zeitpunkt ist die LiNGAM-Modellformel selbst, die durch die Formel $ (2) $ dargestellt wird, fast erhalten, aber die kausale Richtung selbst zum Zeichnen des Graphen ist unbekannt. Daher müssen wir eine Substitution finden, die $ \ bf {B} $ zu einer strengen unteren Dreiecksmatrix macht, indem wir die Reihenfolge der beobachteten Variablen in Gleichung $ (2) $ vertauschen. Das heißt, multiplizieren Sie die Substitutionsmatrix $ \ bf {\ ddot {P}} $ in der Gleichung $ (2) $.

$ \ Bf \ ddot {P} B \ ddot {P} ^ {\ mathrm {T}} $ liegt zu diesem Zeitpunkt so nahe wie möglich an der exakten unteren Dreiecksmatrix $ \ bf {\ ddot {P}} Suchen Sie nach $.

Sie können eine solche Substitutionsmatrix durch Erzwingen der Suche finden, aber die Anzahl der möglichen Substitutionsmatrizen beträgt $ p! $. Es gibt also eine Grenze. In Referenz [3] wird daher ein Verfahren vorgeschlagen, um dies bei hoher Geschwindigkeit zu erhalten.

Ersetzen Sie bei dieser Methode zuerst $ p (p + 1) / 2 $ -Komponenten durch 0 aus $ \ bf \ hat {B} $ in aufsteigender Reihenfolge des absoluten Werts. Danach werden die beobachteten Variablen neu angeordnet, um zu bestimmen, ob es sich um eine untere Dreiecksmatrix handelt oder nicht. Wenn nicht, wird die nächst kleinere Komponente von $ \ bf \ hat {B} $ durch 0 ersetzt und die Bestimmung erneut durchgeführt.

Die Methode der Beurteilung der unteren Dreiecksmatrix (= Erraten der kausalen Ordnung), (1) Lassen Sie zuerst die Zeile von $ \ bf \ hat {B} $, deren Komponenten alle 0 sind, die Zeile $ i $ sein. (Wenn es nicht existiert, ist das Urteil Flase) ② Fügen Sie $ i $ zur Liste der kausalen Ordnungen hinzu ③ Löschen Sie die Zeile $ i $ und die Spalte $ i $ von $ \ bf \ hat {B} $ und kehren Sie als neue Matrix $ \ bf \ hat {B} $ zu ① zurück.

Die vollständige Liste der kausalen Ordnungen lautet Vorfahr → Enkelkind. Nun zur Python-Implementierung. Ehrlich gesagt war dieser Teil zu kompliziert, um ihn zu verarbeiten, also habe ich den Originalcode als Referenz verwendet ...

#SCHRITT 5 Kausale Reihenfolge(causal order)Vermuten

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)/Ersetzen Sie 2 durch 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):    #Strikte Beurteilung der unteren Dreiecksmatrix
    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("Wird keine strenge untere Dreiecksmatrix")
        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

In diesen Daten haben wir eine kausale Ordnung gefunden, die mit einem Schuss zu einer unteren Dreiecksmatrix wird, sodass der Code keine Schleife durchläuft. Eine ordnungsgemäße Implementierung finden Sie im Original.

SCHRITT 6 Fertigstellen und Beschneiden

Wenn die Anzahl der Daten ausreichend ist, ist die in SCHRITT 5 abgeleitete Koeffizientenmatrix in Ordnung. Wenn nicht, scheint das Modell durch Verwendung des Beschneidens als endgültiges Ende verbessert zu werden.

Das Buch schlägt vor, adaptives Lasso als eine der Schnittmethoden zu verwenden. Adaptives Lasso macht die Variablenauswahl konsistent, indem der Lasso-Regularisierungsterm mit $ \ bf w $ gewichtet wird. Es wurde vorgeschlagen, dass dieses $ \ bf w $ die Umkehrung des durch lineare Regression geschätzten Koeffizienten verwendet. Das adaptive Lasso wird für jede beobachtete Variable gemäß der in STEP5 erhaltenen kausalen Reihenfolge durchgeführt, und die endgültige Lösung wird erhalten.

Bei der Implementierung wird der Koeffizient zuerst durch Ausführen einer linearen Regression geschätzt, das Gewicht wird aus dem geschätzten Koeffizienten erhalten, das Gewicht wird auf die Eingabedaten von Lasso angewendet und der geschätzte Betrag, der die Ausgabe von Lasso ist, wird gewichtet.

#SCHRITT 6 Beschneiden
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]    #Predictor Variablen
    target = causal_order[i]    #Zielvariable
    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

Visualisierung

Schließlich wird das erhaltene Modell wie in der folgenden Abbildung dargestellt visualisiert. Von links das wahre Diagramm, das in STEP5 erhaltene Diagramm und das in STEP6 verarbeitete Diagramm. Dieses Mal sind die Daten einfach und wir haben genügend Daten vorbereitet, sodass sie ab STEP5 gut funktionieren.

Graphviz wurde zur Visualisierung verwendet. Ich habe den folgenden Diagrammerstellungscode verwendet.

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

abschließend

Wie oben erwähnt, verglich ein Anfänger des kausalen Denkens Bücher [1] und Papiere [2], untersuchte den Code der Hauptfamilie eingehend und arbeitete hart daran, LiNGAM zu lesen.

Verweise

** [1] ** Shohei Shimizu. Statistische Kausalsuche (professionelle Serie des maschinellen Lernens). 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).