Bonjour, c'est long pour écrire un blog @ kenmatsu4. J'ai écrit l'article le 23ème jour du Calendrier de l'Avent Stan.

Dans ce blog, je vais essayer d'utiliser Stan, une méthode appelée Lasso graphique, qui estime une matrice de précision (matrice inverse d'une matrice distribuée co-distribuée) avec régularisation L1. Le code complet a été téléchargé sur GitHub.

1. Générer des données de test

Tout d'abord, générez un nombre aléatoire qui suit une distribution normale multivariée. Cette fois, nous allons générer 300 données en 6 dimensions avec la moyenne et la variance suivantes. Puis forcez $ x_6 $ et $ x_4 $, puis $ x_6 $ et $ x_5 $ à corréler, de sorte que $ x_4 $ et $ x_5 $ aient une corrélation indirecte. À l'origine, il n'y avait pas de $ x_4 $ et $ x_5 $, mais les valeurs de $ x_4 $ et $ x_5 $ évoluent en conjonction avec la fluctuation de $ x_6 $ sous l'influence de $ x_6 $, donc des variables qui ne sont pas corrélées à l'origine entre elles Apparaîtra avoir une corrélation. ← (*)

`python`


    m = np.zeros(M)
    #Matrice de covariance pour la génération aléatoire d'une distribution normale multivariée
    cov = [[   1, .29, .49, .10, .30,  0],
           [ .29,   1,-.49, .49, .40,  0],
           [ .49,-.49,   1,   0,   0,  0],
           [ .10, .49,   0,   1,   0,  0],
           [ .30, .40,   0,   0,   1,  0],
           [   0,   0,   0,   0,   0,  1]]

    #Créer 6 variables
    X = st.multivariate_normal.rvs(mean=m, cov=cov, size=size, random_state=random_state)

    #Corréler x4 et x6 (*）
    X[:,3] += 0.6*X[:,5]
    
    #Corréler x5 et x6 (*）
    X[:,4] += 0.6*X[:,5]

** Figure: Image de corrélation indirecte ** スクリーンショット 2016-12-23 20.43.46.png

Ce qui suit est un diagramme de dispersion des données obtenues avec pairplot.

La matrice de corrélation et la matrice de corrélation partielle sont calculées comme suit. La corrélation partielle peut être utilisée comme valeur qui indique le degré de corrélation directe, à l'exclusion de l'effet de corrélation indirecte comme indiqué dans la figure précédente. J'aimerais voir immédiatement à quoi ressemble la corrélation partielle.

La matrice de corrélation partielle peut être exprimée à l'aide de l'élément $ \ lambda_ {ij} $ de la matrice de précision (matrice inverse de la matrice distribuée co-distribuée) $ \ Lambda $.

\hat{\rho}_{ij} = {-\lambda_{ij} \over \sqrt{\lambda_{ii}\lambda_{jj}}}

Ce sera.

La droite est la matrice de corrélation partielle, mais si vous calculez à l'aide de la matrice de co-dispersion de dispersion qui est calculée directement à partir des données, elle sera affectée par le bruit et la valeur n'est pas claire. Il ne semble pas que la structure des données soit très visible. Il semble que toutes les variables sont liées ... Cela ne devrait pas être le cas.

Par conséquent, je voudrais exclure l'influence du bruit en utilisant la régularisation L1, c'est-à-dire le Lasso, et utiliser la matrice de co-dispersion de dispersion estimée et la matrice de précision.

Graphical Lasso Un graphe de Markov qui suppose une distribution normale multivariée est appelé un modèle graphique gaussien, et la relation entre les variables peut être considérée comme un modèle graphique utilisant la matrice de précision $ \ Lambda $, qui est un paramètre de cette distribution. Autrement dit, si l'élément i, j $ \ lambda_ {ij} $ de la matrice de précision est non nul, il existe une corrélation directe entre $ x_i $ et $ x_j $. Celui dans un tel état est représenté par le graphique suivant et s'appelle un modèle graphique.

Si les données $ \ boldsymbol {x} $ suivent une distribution normale multivariée, la distribution est

\mathcal{N}(\boldsymbol{x} | {\bf 0}, \Lambda^{-1}) \sim {|\Lambda|^{1/2} \over (2\pi)^{M/2}} \exp\left( -{1 \over 2} \boldsymbol{x}^\intercal \Lambda \boldsymbol{x} \right)

Peut être exprimé comme. Cependant, comme nous avons décidé de considérer qu'il existe une corrélation directe si une valeur autre que 0 est incluse telle quelle, la structure du graphique ci-dessus a une corrélation directe entre la plupart des variables, étant donné que du bruit est ajouté à la plupart des résultats d'estimation. Cela finira par l'être. Il est nécessaire de concevoir de manière à obtenir une matrice de précision la plus éparse possible. Une solution éparse peut être obtenue en supposant une distribution a priori de la distribution de Laplace $ p (\ Lambda) $ dans la matrice de précision $ \ Lambda $. La distribution postérieure est

p(\Lambda|\boldsymbol{x}) \propto p(\Lambda)\prod_{n=1}^{N}\mathcal{N}(\boldsymbol{x}^{(n)}|{\bf 0}, \Lambda)

Par conséquent, en prenant un journal et en le différenciant avec $ \ Lambda $

\ln |\Lambda| - \mathrm{tr}(S\Lambda)-\alpha\|\Lambda\|_1

Par conséquent, définir la distribution antérieure de la distribution de Laplace signifie que la régularisation L1 est appliquée.

3. Trouvez la solution du lasso graphique en utilisant Scikit-Learn

Scikit-Learn implémente Graph Lasso qui implémente ce Lasso graphique. .. Une méthode d'optimisation appelée méthode de descente de coordonnées est utilisée pour cela. Essayons ceci en premier.

C'est très simple à mettre en œuvre. Ajustez simplement comme d'habitude: détendu: Vous pouvez alors obtenir la matrice distribuée co-distribuée et la matrice de précision.

`python`


alpha = 0.2 #Paramètre de régularisation L1
model = GraphLasso(alpha=alpha,
                     max_iter=100,                     
                     verbose=True,
                     assume_centered = True)

model.fit(X)
cov_ = model.covariance_ #Matrice co-distribuée distribuée
prec_ = model.precision_ #Matrice de précision

La matrice de co-distribution de dispersion et la matrice de précision obtenues sont les suivantes.

De plus, la matrice de corrélation calculée à partir de celle-ci et la matrice de corrélation partielle sont les suivantes. Parmi les corrélations $ x_4, x_5, x_6 $ créées par la force, vous pouvez voir que la corrélation entre $ x_4 $ et $ x_5 $ est de 0 sur la matrice de corrélation partielle.

Maintenant, je suis heureux de pouvoir trouver la matrice de corrélation partielle à partir de la matrice de précision clairsemée, mais ce blog est un article de Stan Advent Calendar 2016. Et j'ai écrit plus tôt que cette régularisation L1 peut être interprétée comme supposant une distribution de Laplace dans la distribution antérieure de $ \ Lambda $. Ensuite, vous pouvez essayer ceci avec Stan.

4. Trouvez la solution du lasso graphique en utilisant Stan

Alors, écrivons du code stan et faisons la même chose que ce lasso graphique. En stan, la distribution de Laplace est appelée distribution double exponentielle, nous allons donc l'utiliser. Le code stan est ci-dessous.

`glasso.stan`


data {
  int N;           // Sample size
  int P;           // Feature size
  matrix[N, P] X;  // Data
  real alpha;      // scale parameter of double exponential (L1 parameter)
}
parameters {
  corr_matrix[P] Lambda; // Covariance matrix
}
model {
  vector[P] zeros;
  for (i in 1:P) {
     zeros[i] = 0;
  }
  
  // Precision matrix follows laplace distribution
  to_vector(Lambda) ~ double_exponential(0, 1/alpha);
  
  for (j in 1:N){
    // X follows multi normal distribution
    X[j] ~ multi_normal(zeros, inverse(Lambda));
  }
}
generated quantities {
  matrix[P, P] Sigma;
  Sigma = inverse(Lambda);
}

Voici le code Python qui l'appelle.

`python`


%time sm = pystan.StanModel(file='glasso.stan')
print('Compile finished.')

n_sample = 1000  #Nombre d'échantillons par chaîne
n_warm   = 1000  #Numéro utilisé pour l'échauffement
n_chain  = 4     #Nombre de chaînes
stan_data = {'N': X.shape[0], 'P': P, 'X': X, 'alpha': alpha}
%time fit = sm.sampling(data=stan_data, chains=n_chain, iter=n_sample+n_warm, warmup=n_warm)
print('Sampling finished.')

fit

Voici le résultat. Notez que le Rhat de l'élément diagonal de la matrice co-distribuée distribuée est devenu nan, mais la valeur ne semble pas étrange, et le Rhat de tous les autres éléments est 1,0.

`out`


Inference for Stan model: anon_model_31ac7e216f1b5eccff16f1394bd9827e.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

              mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
Lambda[0,0]    1.0     0.0    0.0    1.0    1.0    1.0    1.0    1.0   4000    nan
Lambda[1,0]  -0.31  8.7e-4   0.05   -0.4  -0.34  -0.31  -0.28  -0.21   2891    1.0
Lambda[2,0]  -0.42  7.1e-4   0.04   -0.5  -0.45  -0.42  -0.39  -0.33   3735    1.0
Lambda[3,0]   0.04  1.0e-3   0.05  -0.07 2.1e-3   0.04   0.07   0.14   2808    1.0
Lambda[4,0]  -0.12  9.1e-4   0.05  -0.22  -0.15  -0.11  -0.08-9.5e-3   3437    1.0
Lambda[5,0]   0.02  1.0e-3   0.06  -0.09  -0.01   0.02   0.06   0.13   3014    1.0
Lambda[0,1]  -0.31  8.7e-4   0.05   -0.4  -0.34  -0.31  -0.28  -0.21   2891    1.0
Lambda[1,1]    1.0 1.5e-189.0e-17    1.0    1.0    1.0    1.0    1.0   3633    nan
Lambda[2,1]   0.47  6.3e-4   0.04   0.39   0.44   0.47    0.5   0.55   4000    1.0
Lambda[3,1]  -0.31  7.6e-4   0.05   -0.4  -0.34  -0.31  -0.28  -0.22   3810    1.0
Lambda[4,1]  -0.19  9.4e-4   0.05  -0.29  -0.22  -0.19  -0.15  -0.08   3021    1.0
Lambda[5,1]    0.2  8.8e-4   0.05    0.1   0.16    0.2   0.23    0.3   3395    1.0
Lambda[0,2]  -0.42  7.1e-4   0.04   -0.5  -0.45  -0.42  -0.39  -0.33   3735    1.0
Lambda[1,2]   0.47  6.3e-4   0.04   0.39   0.44   0.47    0.5   0.55   4000    1.0
Lambda[2,2]    1.0 3.6e-188.7e-17    1.0    1.0    1.0    1.0    1.0    594    nan
Lambda[3,2]  -0.11  8.9e-4   0.05  -0.22  -0.15  -0.11  -0.08  -0.01   3623    1.0
Lambda[4,2]  -0.04  9.1e-4   0.05  -0.15  -0.08  -0.04-5.8e-3   0.07   3642    1.0
Lambda[5,2]   0.03  9.0e-4   0.05  -0.08-9.2e-3   0.03   0.06   0.13   3495    1.0
Lambda[0,3]   0.04  1.0e-3   0.05  -0.07 2.1e-3   0.04   0.07   0.14   2808    1.0
Lambda[1,3]  -0.31  7.6e-4   0.05   -0.4  -0.34  -0.31  -0.28  -0.22   3810    1.0
Lambda[2,3]  -0.11  8.9e-4   0.05  -0.22  -0.15  -0.11  -0.08  -0.01   3623    1.0
Lambda[3,3]    1.0 2.0e-181.2e-16    1.0    1.0    1.0    1.0    1.0   3553    nan
Lambda[4,3]  -0.02  9.3e-4   0.06  -0.13  -0.06  -0.02   0.02   0.09   3503    1.0
Lambda[5,3]  -0.38  7.5e-4   0.04  -0.47  -0.41  -0.38  -0.35  -0.29   3591    1.0
Lambda[0,4]  -0.12  9.1e-4   0.05  -0.22  -0.15  -0.11  -0.08-9.5e-3   3437    1.0
Lambda[1,4]  -0.19  9.4e-4   0.05  -0.29  -0.22  -0.19  -0.15  -0.08   3021    1.0
Lambda[2,4]  -0.04  9.1e-4   0.05  -0.15  -0.08  -0.04-5.8e-3   0.07   3642    1.0
Lambda[3,4]  -0.02  9.3e-4   0.06  -0.13  -0.06  -0.02   0.02   0.09   3503    1.0
Lambda[4,4]    1.0 2.0e-181.2e-16    1.0    1.0    1.0    1.0    1.0   3633    nan
Lambda[5,4]  -0.36  7.2e-4   0.05  -0.45  -0.39  -0.36  -0.33  -0.27   4000    1.0
Lambda[0,5]   0.02  1.0e-3   0.06  -0.09  -0.01   0.02   0.06   0.13   3014    1.0
Lambda[1,5]    0.2  8.8e-4   0.05    0.1   0.16    0.2   0.23    0.3   3395    1.0
Lambda[2,5]   0.03  9.0e-4   0.05  -0.08-9.2e-3   0.03   0.06   0.13   3495    1.0
Lambda[3,5]  -0.38  7.5e-4   0.04  -0.47  -0.41  -0.38  -0.35  -0.29   3591    1.0
Lambda[4,5]  -0.36  7.2e-4   0.05  -0.45  -0.39  -0.36  -0.33  -0.27   4000    1.0
Lambda[5,5]    1.0 2.2e-181.3e-16    1.0    1.0    1.0    1.0    1.0   3381    nan
Sigma[0,0]    1.31  1.1e-3   0.07   1.19   1.26    1.3   1.35   1.45   3507    1.0
Sigma[1,0]    0.26  1.3e-3   0.08   0.11   0.21   0.27   0.32   0.43   4000    1.0
Sigma[2,0]    0.45  1.3e-3   0.08   0.29   0.39   0.44   0.51   0.62   4000    1.0
Sigma[3,0]     0.1  1.2e-3   0.08  -0.05   0.05    0.1   0.15   0.25   4000    1.0
Sigma[4,0]    0.23  1.2e-3   0.08   0.09   0.18   0.23   0.28   0.38   4000    1.0
Sigma[5,0]    0.03  1.2e-3   0.08  -0.13  -0.02   0.03   0.08   0.18   4000    1.0
Sigma[0,1]    0.26  1.3e-3   0.08   0.11   0.21   0.27   0.32   0.43   4000    1.0
Sigma[1,1]    1.55  1.5e-3   0.09   1.38   1.48   1.54   1.61   1.74   4000    1.0
Sigma[2,1]   -0.56  1.4e-3   0.09  -0.74  -0.62  -0.56  -0.49  -0.39   4000    1.0
Sigma[3,1]    0.41  1.3e-3   0.08   0.24   0.35    0.4   0.46   0.57   4000    1.0
Sigma[4,1]    0.29  1.3e-3   0.08   0.14   0.24   0.29   0.34   0.46   4000    1.0
Sigma[5,1]   -0.04  1.3e-3   0.08   -0.2  -0.09  -0.04   0.02   0.13   4000    1.0
Sigma[0,2]    0.45  1.3e-3   0.08   0.29   0.39   0.44   0.51   0.62   4000    1.0
Sigma[1,2]   -0.56  1.4e-3   0.09  -0.74  -0.62  -0.56  -0.49  -0.39   4000    1.0
Sigma[2,2]    1.47  1.3e-3   0.08   1.32   1.41   1.46   1.52   1.65   4000    1.0
Sigma[3,2]  2.9e-3  1.3e-3   0.08  -0.15  -0.05 1.3e-3   0.06   0.16   4000    1.0
Sigma[4,2]    0.04  1.2e-3   0.08  -0.12  -0.02   0.03   0.09   0.19   4000    1.0
Sigma[5,2]    0.07  1.3e-3   0.08  -0.08   0.02   0.08   0.13   0.23   4000    1.0
Sigma[0,3]     0.1  1.2e-3   0.08  -0.05   0.05    0.1   0.15   0.25   4000    1.0
Sigma[1,3]    0.41  1.3e-3   0.08   0.24   0.35    0.4   0.46   0.57   4000    1.0
Sigma[2,3]  2.9e-3  1.3e-3   0.08  -0.15  -0.05 1.3e-3   0.06   0.16   4000    1.0
Sigma[3,3]    1.36  1.1e-3   0.07   1.23    1.3   1.35    1.4   1.51   4000    1.0
Sigma[4,3]    0.31  1.2e-3   0.08   0.17   0.26   0.31   0.36   0.47   4000    1.0
Sigma[5,3]    0.55  1.4e-3   0.09   0.39   0.49   0.55    0.6   0.73   4000    1.0
Sigma[0,4]    0.23  1.2e-3   0.08   0.09   0.18   0.23   0.28   0.38   4000    1.0
Sigma[1,4]    0.29  1.3e-3   0.08   0.14   0.24   0.29   0.34   0.46   4000    1.0
Sigma[2,4]    0.04  1.2e-3   0.08  -0.12  -0.02   0.03   0.09   0.19   4000    1.0
Sigma[3,4]    0.31  1.2e-3   0.08   0.17   0.26   0.31   0.36   0.47   4000    1.0
Sigma[4,4]    1.29  9.9e-4   0.06   1.19   1.25   1.29   1.33   1.43   4000    1.0
Sigma[5,4]    0.53  1.3e-3   0.08   0.38   0.47   0.52   0.58    0.7   4000    1.0
Sigma[0,5]    0.03  1.2e-3   0.08  -0.13  -0.02   0.03   0.08   0.18   4000    1.0
Sigma[1,5]   -0.04  1.3e-3   0.08   -0.2  -0.09  -0.04   0.02   0.13   4000    1.0
Sigma[2,5]    0.07  1.3e-3   0.08  -0.08   0.02   0.08   0.13   0.23   4000    1.0
Sigma[3,5]    0.55  1.4e-3   0.09   0.39   0.49   0.55    0.6   0.73   4000    1.0
Sigma[4,5]    0.53  1.3e-3   0.08   0.38   0.47   0.52   0.58    0.7   4000    1.0
Sigma[5,5]    1.42  1.3e-3   0.08   1.28   1.36   1.41   1.47   1.59   4000    1.0
lp__        -713.2    0.06   2.67 -719.3 -714.9 -712.9 -711.2 -709.0   1983    1.0

Samples were drawn using NUTS at Sat Dec 24 00:05:39 2016.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

C'est un peu difficile à voir avec juste cela, alors dessinons un graphique.

`python`


#Extraction des paramètres estimés
Lambda = fit.extract()["Lambda"]
Sigma  = fit.extract()["Sigma"]

#Calcul de l'estimation du PAE
EAP_Sigma  = np.mean(Sigma, axis=0)
EAP_Lambda = np.mean(Lambda, axis=0)

#Visualisation des estimations du PAE
plt.figure(figsize=(10,4))
ax = plt.subplot(121)
sns.heatmap(pd.DataFrame(EAP_Sigma), annot=EAP_Sigma, fmt='0.2f', ax=ax, xticklabels=label, yticklabels=label)
plt.title("Graphical Lasso with Stan: Covariance matrix")

ax = plt.subplot(122)
sns.heatmap(pd.DataFrame(EAP_Lambda), annot=EAP_Lambda, fmt='0.2f', ax=ax, xticklabels=label, yticklabels=label)
plt.title("Graphical Lasso with Stan: Precision matrix")
plt.savefig(img_path+"glasso_stan_cov_prec.png ", dpi=128)
plt.show()

`python`


#Calcul de la matrice de corrélation
EAP_cor = np.empty_like(EAP_Sigma)
for i in range(P):
    for j in range(P):
        EAP_cor[i, j] = EAP_Sigma[i, j]/np.sqrt(EAP_Sigma[i, i]*EAP_Sigma[j, j])
        
#Calcul de la matrice de corrélation partielle
EAP_rho = np.empty_like(EAP_Lambda)
for i in range(P):
    for j in range(P):
        EAP_rho[i, j] = -EAP_Lambda[i, j]/np.sqrt(EAP_Lambda[i, i]*EAP_Lambda[j, j])
        
plt.figure(figsize=(11,4))
ax = plt.subplot(122)
sns.heatmap(pd.DataFrame(EAP_rho), annot=EAP_rho, fmt='0.2f', ax=ax, xticklabels=label, yticklabels=label)
plt.title("Partial correlation Coefficiant with stan")
#plt.savefig(img_path+"partial_corr_sklearn.png ", dpi=128)

ax = plt.subplot(121)
sns.heatmap(pd.DataFrame(EAP_cor), annot=EAP_cor, fmt='0.2f', ax=ax, xticklabels=label, yticklabels=label)
plt.title("Correlation Coefficiant with stan")
plt.savefig(img_path+"corr_pcorr_stan.png ", dpi=128)
plt.show()

De toute évidence, à cause de la simulation de nombres aléatoires, il n'y a aucun élément qui devient complètement 0. En ce sens, l'effet de la régularisation L1 ne peut pas être obtenu par simulation de Stan.

C'est de Scikit-Learn. En comparaison.

Les valeurs sont légèrement différentes, mais elles ont une structure similaire. La corrélation indirecte entre $ x_4 $ et $ x_5 $ disparaît également dans la matrice de corrélation partielle. L'histogramme échantillonné et le résultat de Scikit-Learn sont superposés et dessinés ci-dessous.

Visualisation des résultats d'échantillonnage

** Histogramme des résultats d'échantillonnage **

Version agrandie. La ligne rouge est le résultat de Scikit-Learn. Les lignes pointillées représentent les points à 2,5% et 97,5% de la distribution postérieure. Certains d'entre eux sont désactivés, mais vous pouvez voir qu'ils sont dans la section à un prix raisonnable. Donc, je pense qu'on peut dire que le résultat est presque le même (hors éléments diagonaux). Tout n'est pas dans la section de conviction, donc un peu plus de réglage peut être nécessaire.

Trace Plot

5. Conclusion

C'est un article que j'ai écrit parce que je savais que la régularisation L1 de Graphical Lasso était une distribution de Laplace de la distribution antérieure des paramètres, et je voulais l'essayer avec Stan. J'ai eu une matrice de corrélation partielle avec une structure similaire, mais il y a une légère différence avec le résultat de Scikit-Learn, donc j'aimerais enquêter un peu plus.

référence

"Détection d'anomalies et détection de changement (série professionnelle d'apprentissage automatique)" Takeshi Ide, Masaru Sugiyama Stan Modeling Language User’s Guide and Reference Manual 　⇒ http://www.uvm.edu/~bbeckage/Teaching/DataAnalysis/Manuals/stan-reference-2.8.0.pdf Coefficient de corrélation partielle 　⇒ http://www.ae.keio.ac.jp/lab/soc/takeuchi/lectures/5_Parcor.pdf

[PYTHON] [PyStan] Essayez le lasso graphique avec Stan.

1. Générer des données de test

python

3. Trouvez la solution du lasso graphique en utilisant Scikit-Learn

python

4. Trouvez la solution du lasso graphique en utilisant Stan

glasso.stan

python

out

python

python

Visualisation des résultats d'échantillonnage

5. Conclusion

référence

`python`

`python`

`glasso.stan`

`python`

`out`

`python`

`python`