[Python] PCA scratch dans l'exemple de "Introduction à la méthode d'analyse multivariée"

introduction

J'ai essayé de suivre la preuve de l'analyse en composantes principales, mais je ne comprenais pas comment l'utiliser, j'ai donc décidé d'approfondir ma compréhension en grattant l'exemple 4 du chapitre 9 du livre "Introduction to Multivariate Analysis" publié par Science. Cet article est destiné à être ce rappel.

Quoi utiliser

・ VSCode (... eh bien, n'importe quel éditeur va bien) ・ Python version 3.7.0 ·Bibliothèque -Numpy (pour gérer la matrice) -Matplotlib (utilisé dans le diagramme de dispersion) ・ Introduction à la méthode d'analyse multivariée (Science) ← Dans ce qui suit, ce livre sera appelé un manuel.

couler

Vérifiez facilement la politique de scratch. ① Saisie des données ② Standardisation des données ③ Calculez la matrice des coefficients de corrélation ④ Calculer la valeur propre et le vecteur propre ⑤ Effectuer l'analyse des composants principaux et compresser en 2 dimensions ⑥ De plus, le facteur de charge est également indiqué. ⑦ Dessinez un diagramme de dispersion de chaque chargement factoriel et du score de la composante principale

Matrice de données

Les données de l'exemple, je les ai saisies à la main.

ex9_4.txt

86 79 67 68
71 75 78 84
42 43 39 44
62 58 98 95
96 97 61 63
39 33 45 50
50 53 64 72
78 66 52 47
51 44 76 72
89 92 93 91

code

Je vais mettre le code.

scratch.py

import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt

#Entrez la taille de l'échantillon(ligne)
n = int(input())

#Entrer des données
a = [input().split() for l in range(n)]
data = np.array(a, dtype=float)
line = int(data.shape[0])
row = int(data.shape[1])
print(line, row)

means = np.empty(row)
sds = np.empty(row)

#Calculer la moyenne
sum = 0
for i in range(row):
    for j in range(line):
        sum +=  data[j][i]
    else:
        means[i] = sum / line
        sum = 0
print("moyenne:")
print(means)

#Calculer l'écart type
for i in range(row):
    for j in range(line):
        sum += (data[j][i]-means[i]) ** 2
    else:
        sds[i] = (sum / (line - 1)) ** 0.5
        sum = 0
print("écart-type:")
print(sds)

#Standardisation
standardized_data = np.zeros((line,row))
for i in range(row):
    for j in range(line):
        standardized_data[j][i] = (data[j][i] - means[i]) / sds[i]  
print(standardized_data)

#Calculer la matrice des coefficients de corrélation
cc_matrix = np.zeros((row, row))
for i in range(row):
    for j in range(row):
        if(cc_matrix[i][j] == 0):
            for k in range(line):
                sum += standardized_data[k][i] * standardized_data[k][j]
            else:
                cc_matrix[i][j] = sum / (line - 1)
                cc_matrix[j][i] = cc_matrix[i][j]
                sum = 0

#Problème de valeur unique
w,v = LA.eig(cc_matrix)
#print(w)
#print(v)

#Compression en 2 dimensions par analyse en composantes principales
main_component = np.zeros((line, 2))
z_1, z_2 = 0, 0
for i in range(line):
    for j in range(row):
        z_1 += standardized_data[i][j] * v[j][0]
        z_2 += standardized_data[i][j] * v[j][1]
    else:
        main_component[i][0] = z_1
        main_component[i][1] = z_2
        z_1, z_2 = 0, 0
print("Composant principal:")
print(main_component)

#Facteur de chargement
factor_load = np.zeros((2, row))
for i in range(2):
    for j in range(row):
        factor_load[i][j] = (w[i] ** 0.5) * v[j][i]
print("Facteur de chargement:")
print(factor_load)

#Dessinez un diagramme de dispersion

#Diagramme de dispersion du chargement factoriel
for label in range(row):
    plt.scatter(factor_load[0][label], factor_load[1][label], marker="$"+str(label+1)+"$")
    plt.title('factor_load')
    plt.xlabel('pc1')
    plt.ylabel('pc2')
plt.show()

#Diagramme de dispersion des scores des composants principaux
for label in range(line):
    plt.scatter(main_component[label][0], main_component[label][1], marker="$"+str(label+1)+"$")
    plt.title('principal component')
    plt.xlabel('pc1')
    plt.ylabel('pc2')
plt.show()

Chaque processus du code est expliqué ci-dessous.

Saisie des données

Tout d'abord, utilisez input () pour saisir le nombre de lignes de la matrice à créer à partir du clavier. Il s'agit de la taille de l'échantillon dans les données d'origine. Ensuite, saisissez les données de la matrice avec le clavier et gérez-les avec numpy. Je pense qu'il existe une autre bonne méthode pour cette méthode de saisie (mais cette fois, cette partie est comparée à l'analyse des composants principaux cibles. Je pensais que c'était une question triviale, et je ne voulais pas passer trop de temps, alors j'ai décidé de le mettre sur l'étagère et de continuer.) Ensuite, obtenez le nombre de lignes et de colonnes de la matrice créée. Ces types l'utilisent beaucoup. De plus, avec ce code, le vecteur moyen et le vecteur d'écart type sont initialisés à ce stade.

#Entrez la taille de l'échantillon(ligne)
n = int(input())

#Entrer des données
a = [input().split() for l in range(n)]
data = np.array(a, dtype=float)
line = int(data.shape[0])
row = int(data.shape[1])
print(line, row)

means = np.empty(row)
sds = np.empty(row)

Standardisation des données

À partir de la matrice de données, trouvez la moyenne et l'écart type de chaque colonne et utilisez-les pour créer une matrice qui standardise les données d'origine. (Standardized_data) À partir de là, c'est essentiellement une tempête de phrases, j'ai beaucoup utilisé ma tête en indice. Notez également que lors du calcul de l'écart type, le nombre à diviser doit être n-1 au lieu de n, sinon la valeur ne sera pas celle du manuel.

#Calculer la moyenne
sum = 0
for i in range(row):
    for j in range(line):
        sum +=  data[j][i]
    else:
        means[i] = sum / line
        sum = 0
print("moyenne:")
print(means)

#Calculer l'écart type
for i in range(row):
    for j in range(line):
        sum += (data[j][i]-means[i]) ** 2
    else:
        sds[i] = (sum / (line - 1)) ** 0.5
        sum = 0
print("écart-type:")
print(sds)

#Standardisation
standardized_data = np.zeros((line,row))
for i in range(row):
    for j in range(line):
        standardized_data[j][i] = (data[j][i] - means[i]) / sds[i]  
print(standardized_data)

Calculer la matrice des coefficients de corrélation

Nous apprécions les équations du manuel (9.2) et (9.3) et calculons la matrice des coefficients de corrélation. Bien entendu, nous confirmons la preuve sur papier avant d'utiliser les équations. De plus, nous essayons de réduire la quantité de calcul en utilisant le fait que la matrice des coefficients de corrélation est une matrice symétrique.

#Calculer la matrice des coefficients de corrélation
cc_matrix = np.zeros((row, row))
for i in range(row):
    for j in range(row):
        if(cc_matrix[i][j] == 0):
            for k in range(line):
                sum += standardized_data[k][i] * standardized_data[k][j]
            else:
                cc_matrix[i][j] = sum / (line - 1)
                cc_matrix[j][i] = cc_matrix[i][j]
                sum = 0

Problème de valeur unique

Je vais calculer la valeur propre, le vecteur propre ... mais je ne pouvais pas penser à comment gratter ici, alors j'ai décidé de le jeter dans la fonction linalg.eig et de continuer.Je mets la valeur propre dans w et le vecteur propre dans v. Masu. Si vous avez un bon moyen de gratter ici, faites-le moi savoir.

#Problème de valeur unique
w,v = LA.eig(cc_matrix)
#print(w)
#print(v)

Enfin l'analyse des composants principaux

C'est le sujet principal Cette fois, je veux vérifier le flux jusqu'au diagramme de dispersion, j'ai donc décidé de supprimer deux dimensions, c'est-à-dire deux composants principaux, et j'ai écrit le code. car la déclaration est difficile ...

#Compression en 2 dimensions par analyse en composantes principales
main_component = np.zeros((line, 2))
z_1, z_2 = 0, 0
for i in range(line):
    for j in range(row):
        z_1 += standardized_data[i][j] * v[j][0]
        z_2 += standardized_data[i][j] * v[j][1]
    else:
        main_component[i][0] = z_1
        main_component[i][1] = z_2
        z_1, z_2 = 0, 0
print("Composant principal:")
print(main_component)

Facteur de chargement

En passant, c'est un facteur de chargement.

#Facteur de chargement
factor_load = np.zeros((2, row))
for i in range(2):
    for j in range(row):
        factor_load[i][j] = (w[i] ** 0.5) * v[j][i]
print("Facteur de chargement:")
print(factor_load)

Nuage de points

Je vais tracer le score de la composante principale et la quantité de chargement de facteur sur lesquels j'ai travaillé dur sur le diagramme de dispersion. L'axe des x est le premier composant principal et l'axe des y est le deuxième composant principal.

#Dessinez un diagramme de dispersion

#Diagramme de dispersion du chargement factoriel
for label in range(row):
    plt.scatter(factor_load[0][label], factor_load[1][label], marker="$"+str(label+1)+"$")
    plt.title('factor_load')
    plt.xlabel('pc1')
    plt.ylabel('pc2')
plt.show()

#Diagramme de dispersion des scores des composants principaux
for label in range(line):
    plt.scatter(main_component[label][0], main_component[label][1], marker="$"+str(label+1)+"$")
    plt.title('principal component')
    plt.xlabel('pc1')
    plt.ylabel('pc2')
plt.show()

Nuage de points

J'ai pu dessiner à peu près la même chose qu'un manuel. Est-ce bien

À la fin

J'ai écrit un article sur la qiita pour la première fois, je pense qu'il y a beaucoup de choses difficiles à lire. Aussi, parce que je suis immature, je pense qu'il y a beaucoup de tsukkomi dans le code ci-dessus. Je suis très bienvenu à tsukkomi pour eux. Plutôt, je l'ai posté sur qiita pour obtenir tsukkomi.

Les références

Yasushi Nagata, Masahiko Muchinaka Introduction to Multivariate Analysis Method Science, 10 avril 2001