introduction

Cet article est un rappel de certaines tentatives faites dans le processus de prise en compte de la corrélation des données avec les variables d'erreur et les facteurs d'intrication.

Tout d'abord, je vais montrer la structure de cet article.

Structure causale des données ([Data Generation](# data-generation))
Corrélation / Corrélation partielle 1 ([Corrélation / Corrélation partielle [1]](# correlation - partial-correlation-1))
Estimation parcimonieuse ([Graphical Lasso](# graphical-lasso))
Corrélation / Corrélation partielle 2 ([Correlation / Partial Correlation [2]](# correlation - partial-correlation-2))
Analyse de régression ([Analyse de régression](# Analyse de régression))
Résumé (Summary)

Data Generation

En ce qui concerne l'échelle des données, «échelle nominale», «échelle d'ordre», «variable de catégorie», etc. ne sont pas incluses, et concernant la nature et la structure des données, «biais», «facteur d'enchevêtrement non observé», «patrouille dirigée», etc. pas ici. On suppose que toutes les données multivariées à traiter sont des valeurs continues et chacune suit une distribution normale. Nous normalisons la «granulométrie» des données.

import pandas as pd
import numpy as np
import math
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
import scipy as sp
from scipy import linalg
from scipy.stats import multivariate_normal
from sklearn.linear_model import LinearRegression, HuberRegressor, TheilSenRegressor, RANSACRegressor
from sklearn.svm import SVR
from sklearn.covariance import GraphicalLasso
from sklearn.preprocessing import StandardScaler 
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, mean_squared_log_error

<Génération des données à utiliser>

Génère des données avec la structure causale dirigée et non circulante suivante. Plutôt que l'effet de facteurs non observés, des nombres aléatoires qui suivent une distribution normale ou une distribution uniforme sont ajoutés à chaque nœud en tant que variables d'erreur.

np.random.seed(0)
data_ = np.array([[20, 1, 0.3242, 34, 1, 0.10, 52, 5], [52, 0, 0.6134, 63, 1, 0.35, 38, 8], 
                  [38, 0, 0.7136, 57, 0, 0.20, 64, 7], [63, 1, 0.4592, 49, 0, 0.80, 95, 2]])
mu = np.mean(data_, axis = 0)
cov = np.cov(data_, rowvar = 0, bias = 0)

w_9 = np.random.multivariate_normal(mu, cov, 1000)[:, 0][:, np.newaxis] + np.random.uniform(1.3, size = 1000)[:, np.newaxis]
w_8 = 0.32 * w_9 + np.random.multivariate_normal(mu, cov, 1000)[:, 3][:, np.newaxis]
w_7 = np.random.multivariate_normal(mu, cov, 1000)[:, 1][:, np.newaxis] + np.random.uniform(1.6, size = 1000)[:, np.newaxis]
w_6 = 0.50 * w_7 + np.random.multivariate_normal(mu, cov, 1000)[:, 4][:, np.newaxis]
w_5 = 0.53 * w_7 + np.random.multivariate_normal(mu, cov, 1000)[:, 5][:, np.newaxis]
w_4 = 0.48 * w_5 + 0.40 * w_6 + np.random.uniform(1.5, size = 1000)[:, np.newaxis]
w_3 = 0.28 * w_4 + np.random.uniform(1.7, size = 1000)[:, np.newaxis]
w_2 = 0.20 * w_5 + np.random.uniform(3.2, size = 1000)[:, np.newaxis]
w_1 = 0.31 * w_2 + np.random.uniform(1.5, size = 1000)[:, np.newaxis]
w_0 = 0.290 * w_1 + 0.327 * w_4 + 0.491 * w_8 + 0.136 * w_9 + 10.0 + np.random.multivariate_normal(mu, cov, 1000)[:, 2][:, np.newaxis]

mu = np.mean(w_0, axis = 0)
cov = np.cov(w_0, rowvar = 0, bias = 0)
list_ = np.where(w_0 <= 0)
for i in list_:
    normal = np.random.normal(mu, cov)
    if normal <= 0:
        w_0[i] = mu

    elif normal > 0:
        w_0[i] = normal

data_ = pd.DataFrame(np.concatenate([w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9], axis = 1))

Correlation / Partial Correlation [1]

Idéalement, la matrice de corrélation devrait donner le résultat que «w_7» a une corrélation entre «w_5» et «w_6», et la matrice de corrélation partielle devrait donner le résultat que la corrélation est une pseudo-corrélation. En même temps, je suis également intéressé par la manière dont les résultats de la corrélation entre «w_0» et «w_8» changent entre la matrice de corrélation et la matrice de corrélation partielle. Si la structure de données n'est pas compliquée comme la relation de «w_0» w_8 »« w_9 », même si le facteur d'intrication« w_9 »existe, s'il n'y a pas de relation causale entre« w_0 »et« w_8 », la matrice de corrélation partielle est utilisée. Le calcul doit donner le résultat que «w_0» et «w_8» n'ont pas de relation causale. D'autre part, dans les données à traiter, il y a une relation causale entre «w_0» et «w_8», et il y a «w_9», qui est un facteur d'intrication pour «w_0» et «w_8». Dans ce cas, voyons comment la corrélation de w_0`` w_8 est exprimée. Vérifiez la précision de la corrélation pour w_0 et w_1`` w_4`` w_8 w_9.

def partital_corr(data):
    temp_cov = data.cov()
    omega = np.linalg.inv(temp_cov)
    D = np.diag(np.power(np.diag(omega), -0.5))
    temp_pcorr = -np.dot(np.dot(D, omega), D) + 2 * np.eye(temp_cov.shape[0])
    mtx_pcorr = pd.DataFrame(temp_pcorr, columns = temp_cov.columns, index = temp_cov.index)
    return mtx_pcorr


partial_correlation = partital_corr(data_)

fig = plt.subplots(figsize=(9, 9))
heatmap = sns.heatmap(data_.corr(), annot = True, cmap = 'Blues', square = True)
heatmap_ = sns.heatmap(partial_correlation, annot = True, cmap = 'Blues', square = True)

Correlation(left) VS Partial Correlation(right)

	w_0``w_8	w_0``w_9	w_0``w_1	w_0``w_4	w_5``w_6	w_8``w_9
corrélation	0.96	0.66	-0.031	0.014	0.36	0.43
Corrélation partielle	1	1	0.28	0.24	-0.6	-1

Une interprétation est lors du calcul de la matrice de variance-covariance dans le processus de calcul de la matrice de corrélation partielle parce que les deux «w_5» et «w_6», qui sont affectés par le facteur d'intrication «w_7», affectent «w_4». Il semble que l'on puisse considérer que la relation entre les variables n'est pas bien mesurée. Pour une autre considération, Corrélation / Corrélation partielle [2] ajoute une structure de boucle et des valeurs aberrantes aux données.

Graphical Lasso Le lasso graphique est une estimation creuse qui calcule l'inverse d'une matrice de covariance distribuée creuse lorsque les données suivent une distribution normale multivariée. Il contient des éléments du modèle graphique gaussien et une estimation parcimonieuse. Considérez l'indépendance conditionnelle par rapport à d'autres variables dans la corrélation entre deux variables. Matrice inverse (matrice de précision) de matrice co-distribuée éparse en créant une fonction de vraisemblance logarithmique à partir de la matrice co-distribuée distribuée calculée à partir des données et en résolvant le problème de minimisation de la fonction de perte avec le terme de régularisation $ L1 $ ajouté. Estimer. Le lasso graphique ne convient pas aux données que vous traitez en devinant, mais je vais l'essayer.

ss = StandardScaler()
data_ss = ss.fit_transform(data_)
gl = GraphicalLasso(alpha = 0.1, max_iter = 100, )
gl.fit(data_ss)
gl.get_params()

# correlation
diag = np.zeros(shape = (10, 10))
for i in range(diag.shape[0]):
    for j in range(diag.shape[1]):
        if i == j:
            diag[i,j] = 1 / np.sqrt(gl.covariance_[i,j])
corr = np.dot(np.dot(diag, gl.covariance_), diag)

# partial correlation
omega = np.linalg.inv(gl.covariance_)
diag = np.diag(np.power(np.diag(omega), -0.5))
partial_corr = -np.dot(np.dot(diag, omega), diag)
partial_corr += 2 * np.eye(gl.covariance_.shape[0])

heatmap_gl =sns.heatmap(corr, annot = True, cmap = 'Blues', square = True)
figure_gl = heatmap_gl.get_figure()
figure_gl.savefig('graph_gl_corr')

heatmap_gl_ = sns.heatmap(partial_corr, annot = True, cmap = 'Blues', square = True)
figure_gl_ = heatmap_gl_.get_figure()
figure_gl_.savefig('graph_gl_partial_corr')

Correlation(left) VS Partial Correlation(right) Le lasso graphique n'est pas un algorithme qui prend en charge le bruit et l'intrication, il est donc considéré que ce n'est pas un choix approprié pour ces données.

Correlation / Partial Correlation [2] Créez une structure de boucle dans w_0`` w_4`` w_6 et créez des données w_1 w_4`` w_9 avec des valeurs aberrantes.

np.random.seed(0)
data_ = np.array([[20, 1, 0.3242, 34, 1, 0.10, 52, 5], [52, 0, 0.6134, 63, 1, 0.35, 38, 8], 
                  [38, 0, 0.7136, 57, 0, 0.20, 64, 7], [63, 1, 0.4592, 49, 0, 0.80, 95, 2]])
mu = np.mean(data_, axis = 0)
cov = np.cov(data_, rowvar = 0, bias = 0)

#w_11 = np.random.multivariate_normal(mu, cov, 1000)[:, 6][:, np.newaxis] 
#w_10 = 0.27 * w_11 + np.random.multivariate_normal(mu, cov, 1000)[:, 7][:, np.newaxis] 
w_9 = np.random.multivariate_normal(mu, cov, 1000)[:, 0][:, np.newaxis] + np.random.uniform(1.3, size = 1000)[:, np.newaxis]

noise_upper = np.random.normal(200, 30, 20)
noise_row = np.random.randint(0, 1000, 20)
w_9[noise_row] = noise_upper[:, np.newaxis]

w_8 = 0.32 * w_9 + np.random.multivariate_normal(mu, cov, 1000)[:, 3][:, np.newaxis]
w_7 = np.random.multivariate_normal(mu, cov, 1000)[:, 1][:, np.newaxis] + np.random.uniform(1.6, size = 1000)[:, np.newaxis]
w_6 = 0.50 * w_7 + np.random.multivariate_normal(mu, cov, 1000)[:, 4][:, np.newaxis] + 0.53 * w_0
w_5 = 0.53 * w_7 + np.random.multivariate_normal(mu, cov, 1000)[:, 5][:, np.newaxis]
w_4 = 0.48 * w_5 + 0.40 * w_6 + np.random.uniform(1.5, size = 1000)[:, np.newaxis]

noise_lower = np.random.normal(0.6, 0.2, 10)
noise_upper = np.random.normal(30, 2, 30)
noise_row = np.random.randint(0, 1000, 40)
w_4[noise_row] = np.concatenate([noise_lower, noise_upper])[:, np.newaxis]

w_3 = 0.28 * w_4 + np.random.uniform(1.7, size = 1000)[:, np.newaxis]
w_2 = 0.20 * w_5 + np.random.uniform(3.2, size = 1000)[:, np.newaxis]
w_1 = 0.31 * w_2 + np.random.uniform(1.5, size = 1000)[:, np.newaxis]

noise_upper = np.random.normal(30, 4, 20)
noise_upper_ = np.random.normal(50, 3, 20)
noise_row = np.random.randint(0, 1000, 40)
w_1[noise_row] = np.concatenate([noise_upper, noise_upper_])[:, np.newaxis]

w_0 = 0.290 * w_1 + 0.327 * w_4 + 0.491 * w_8 + 0.136 * w_9 + 10.0 + np.random.multivariate_normal(mu, cov, 1000)[:, 2][:, np.newaxis] #+ 0.426 * w_11

mu = np.mean(w_0, axis = 0)
cov = np.cov(w_0, rowvar = 0, bias = 0)
list_ = np.where(w_0 <= 0)
for i in list_:
    normal = np.random.normal(mu, cov)
    if normal <= 0:
        w_0[i] = mu

    elif normal > 0:
        w_0[i] = normal

data_ = pd.DataFrame(np.concatenate([w_0, w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8, w_9], axis = 1))

def partital_corr(data):
    temp_cov = data.cov()
    omega = np.linalg.inv(temp_cov)
    D = np.diag(np.power(np.diag(omega), -0.5))
    temp_pcorr = -np.dot(np.dot(D, omega), D) + 2 * np.eye(temp_cov.shape[0])
    mtx_pcorr = pd.DataFrame(temp_pcorr, columns = temp_cov.columns, index = temp_cov.index)
    return mtx_pcorr


partial_correlation = partital_corr(data_)

fig = plt.subplots(figsize=(9, 9))
heatmap = sns.heatmap(data_.corr(), annot = True, cmap = 'Blues', square = True)
heatmap_ = sns.heatmap(partial_correlation, annot = True, cmap = 'Blues', square = True)

Correlation(left) VS Partial Correlation(right) Puisque «w_0» et «w_4» sont connectés dans une structure en boucle, le résultat est que «w_4» est en corrélation avec «w_8» et «w_9» respectivement. La matrice de corrélation montre une corrélation positive et la matrice de corrélation partielle montre une corrélation négative. Dans le résultat de la matrice de corrélation, «w_9» montre indirectement une corrélation positive par rapport à «w_3» dans la corrélation à «w_3» w_4 »« w_9 ». Dans les résultats de la matrice de corrélation partielle, «w_3» et «w_9» ne sont pas corrélés, et «w_4» et «w_9» sont corrélés négativement. Vous pouvez voir certains endroits où la matrice de corrélation est compatible, mais pas la matrice de corrélation partielle. L'inverse est également vrai. Vous pouvez également voir que la matrice de corrélation montre une corrélation positive et la matrice de corrélation partielle montre une corrélation négative.

Il n'y a pas suffisamment de preuves pour donner une considération suffisante, mais une interprétation est que dans le cas de cette section, la matrice de corrélation et la matrice de corrélation partielle sont toutes deux des «variables d'erreur d'une certaine taille ou plus» et «un facteur (variable)». Deux ou plusieurs variables $ x_1, ..., x_n $ qui sont indirectement ou directement affectées par sont causalement ou corrélées, ou deux ou plus de $ x_1, ..., x_n $ sont autres Il semble que l'on puisse considérer que "en affectant la même variable" et "la structure en boucle" ne sont pas bien gérés.

analyse de régression

Puisque «w_0» est fait de sorte qu'il puisse être exprimé par une expression linéaire, une analyse de régression est effectuée. Si vous pouvez supposer «normalité», «dispersion égale» et «indépendance» pour la distribution d'erreur, et «linéarité» pour la relation entre la variable objective et la variable explicative, vous devriez avoir une forte motivation pour effectuer une analyse de régression. Je peux le faire. J'ai aussi essayé «HuberRegressor» et «SVR» avec «scikit-learn», mais «LinearRegression» est suffisant. En effet, les données traitées ne présentent pas de biais tels que des biais de sélection ou des facteurs d'intrication non observés. Lorsque nous faisons SVR, nous utilisons un noyau linéaire qui pondère le vecteur d'entrée. La précision peut être obtenue même si la régression est effectuée par la méthode d'ajustement du degré de liberté à l'aide d'un noyau polymorphe correspondant à plusieurs éléments ou la méthode utilisant un noyau gaussien. La régression Hoover est utilisée pour une régression robuste. Cette fois, la distribution d'erreur n'a pas d'asymétrie ou de bruit important, elle n'est donc pas considérée comme un choix approprié. Puisqu'il n'y a pas de motivation telle que "faire en sorte que la variable objective suive une distribution normale", "faire face aux valeurs aberrantes", et "la relation entre la variable objective et la variable explicative a tendance à être une fonction exponentielle", la conversion logarithmique n'est pas effectuée.

X_train, X_valid, y_train, y_valid = train_test_split(data_.drop(0, axis =1), data_[0], test_size = 0.4, random_state = 0)
data_ss = pd.DataFrame(data_ss)
X_train_ss, X_valid_ss, y_train_ss, y_valid_ss = train_test_split(data_ss.drop(0, axis =1), data_ss[0], test_size = 0.4, random_state = 0)
y_valid = y_valid.sort_index(ascending = True)
X_valid = X_valid.sort_index(ascending = True)
y_valid_ss = y_valid_ss.sort_index(ascending = True)
X_valid_ss = X_valid_ss.sort_index(ascending = True)

plt.hist(np.ravel(y_train), label = 'training data', color = 'c')
plt.hist(np.ravel(y_valid), label = 'valid data', color = 'k')
plt.savefig('graph_hist')

#Données non standardisées
linear = LinearRegression()
linear.fit(X_train, y_train)
# <SVR>
# svr_linear = SVR(kernel = 'linear', C = 1.0, epsilon = 0.1, gamma = 'scale')
# svr_linear.fit(X_train, np.ravel(y_train))
# scores = svr_linear.score(X_valid, y_valid)
# print('R^2 scores : %s' %scores)
# <HuberRegressor>
# hr = HuberRegressor(alpha = 0.0001, epsilon = 1.35, # fit_intercept = True, max_iter = 100, tol = 1e-05, warm_start = False)
# hr.fit(X_train, y_train)
pred = linear.predict(X_valid)
plt.plot(np.ravel(y_valid), label = 'true valid', color = 'k')
plt.plot(np.ravel(pred), label = 'pred valid', color = 'w')

diff = np.log1p(pred) - np.log1p(y_valid)
diff = np.where(np.isnan(diff) == True, 0, diff)
RMSLE = np.sqrt(np.mean(diff ** 2))

#Données standardisées
linear_ss = LinearRegression()
linear_ss.fit(X_train_ss, y_train_ss)
# <SVR>
# svr_linear_ss = SVR(kernel = 'linear', C = 1.0, epsilon = 0.1, gamma = 'scale')
# svr_linear_ss.fit(X_train_ss, np.ravel(y_train_ss))
# scores = svr_linear_ss.score(X_valid_ss, y_valid_ss)
# print('R^2 scores : %s' %scores)
# <HuberRegressor>
# hr_ss = HuberRegressor(alpha = 0.0001, epsilon = 1.35, # fit_intercept = True, max_iter = 100, tol = 1e-05, warm_start = False)
# hr_ss.fit(X_train_ss, y_train_ss)
pred_ss = linear_ss.predict(X_valid_ss)
plt.plot(np.ravel(y_valid_ss), label = 'true valid', color = 'k')
plt.plot(np.ravel(pred_ss), label = 'pred valid', color = 'w')

coef = pd.DataFrame(X_train.columns.values, columns = {'feature'})
coef['coef'] = linear.coef_.T
coef['coef_ss'] = linear_ss.coef_.T

<Résultats utilisant des données non standardisées>

print('MSE : %s' %mean_squared_error(y_valid, pred, squared = True))
print('RMSLE : %s' %RMSLE)
>> RMSE : 0.02983134938758614
>> RMSLE : 0.03423681936352853

coef.sort_values(by = 'coef', ascending = False)
>>
feature     coef
	8	0.480333
	4	0.322978
	1	0.284168
	9	0.133841
	6	0.039120
	5	0.016435
	3	-0.005714
	2	-0.009774
	7	-0.020549

<Résultats utilisant des données standardisées>

print('MSE : %s' %mean_squared_error(y_valid_ss, pred_ss, squared = True))
>> MSE : 0.000242675476861591

coef.sort_values(by = 'coef_ss', ascending = False)
>>
feature     coef
	8	0.652601
	9	0.332685
	1	0.194995
	4	0.109338
	6	0.020749
	5	0.000665
	3	-0.000551
	2	-0.000566
	7	-0.001129

Vérifiez les résidus.

#Données non standardisées
plt.scatter(pred, pred - y_valid, color = 'blue') 
plt.hlines(y = 0, xmin = -10, xmax = 150, color = 'c')
plt.xlabel('Prediction')
plt.ylabel('Residual error')
plt.grid()

#Données standardisées
plt.scatter(pred_ss, pred_ss - y_valid_ss, color = 'blue') 
plt.hlines(y = 0, xmin = -10, xmax = 150, color = 'c')
plt.xlabel('Prediction')
plt.ylabel('Residual error')
plt.grid()

Non Standardized(left) VS Standardized(right)

#Données non standardisées
plt.hist(pred - y_valid)

#Données standardisées
plt.hist(pred_ss - y_valid_ss)

Non Standardized(left) VS Standardized(right)

Les résidus des données standardisées et non standardisées sont dispersés dans la plage de $ ± 0,5 $ et $ ± 0,05 $ de $ y = 0 $, respectivement, et vous pouvez vérifier la normalité des résidus. La façon dont les résidus sont dispersés n'est pas irrégulière et est symétrique avec $ y = 0 $ comme limite. Les données normalisées ont une précision de prédiction plus élevée pour les données de vérification par MSE que les données non normalisées, mais lorsque la contribution du coefficient à la variable objective est confirmée par ordre décroissant, les données non normalisées 8 4 1 Vous pouvez voir que la séquence 9` est plus précise que les données standardisées, qui ont une précision de prédiction élevée pour les données de vérification. Même si des variables d'erreur et des valeurs aberrantes sont ajoutées aux données, si les données n'ont pas d'échelles multiples ou de facteurs d'intrication non observés, la précision est raisonnable. Même si des variables d'erreur et des valeurs aberrantes plus aléatoires et plus importantes sont attribuées, la contribution du coefficient à la variable objective est calculée de manière appropriée dans ces données, et le résultat est que l'exactitude du modèle n'est pas stable lorsqu'il y a des variables non observées. Il est sorti.

L'analyse régressive est effectuée à l'aide de «statsmodels.formula.api» pour les données qui prennent en compte les différences dans plusieurs échelles et la granularité lorsque l'accent est mis sur la prédiction plutôt que sur les relations causales, et «Cond pour la colinéarité multiple. Si le nombre de conditions est proche de 1, il est jugé que la colinéarité multiple est faible. Sinon, le coefficient d'expansion de la variance peut être calculé en utilisant variance_inflation_factor de statsmodels.stats.outliers_influence. .. Par exemple, lorsqu'il existe une colinéarité multiple entre deux variables, le plan de régression n'existe que sur la ligne droite car les variables sont disposées en ligne droite et la matrice inverse de la matrice co-distribuée distribuée n'existe pas. Je ne peux pas résoudre. Il en résulte que "l'écart type du coefficient devient grand", "la valeur $ t $ devient petite et affecte la valeur $ p $", "le coefficient de détermination prend une valeur grande", et ainsi de suite. En utilisant statsmodels.formula.api," déterminez l'effet de fond en utilisant la valeur $ p $, la valeur $ z $ et l'intervalle de confiance de la population "," Type de covariance "et" pour les données d'échantillon. Vous pouvez facilement vérifier le test de normalité (test Omnibus).

Résumé

J'ai examiné la corrélation et effectué une analyse de régression sur des données avec des variables d'erreur et des facteurs d'intrication. Bien qu'il s'agisse d'un mémorandum, j'estime que le contenu de l'article est incohérent par rapport au sujet. Si nous pouvons organiser correctement les données qui peuvent être utilisées comme article, ce sera plus intéressant que d'estimer la cause et l'effet par la modélisation d'équations structurelles.

Mémorandum sur la corrélation [Python]

introduction

analyse de régression

Résumé