introduction

Après avoir passé le test statistique de niveau 2 en novembre 2019, j'ai commencé à apprendre régulièrement avec l'intention de viser également le niveau 1. D'après ce que j'entends, il semble nécessaire de maintenir l'analyse multivariée afin de passer en quasi-première année. Donc, je vais lire le livre de référence et apprendre l'analyse multivariée, mais je pense qu'il est préférable de bouger mes mains pour approfondir ma compréhension, donc je vais essayer diverses méthodes d'analyse utilisant Python. Je déciderai.

Cependant, le titre a été réglé sur [Test statistique niveau 2 / niveau 1] car il peut être un peu impliqué dans le problème de lecture des résultats d'analyse du logiciel de traitement statistique (R), que j'avais en quelque sorte traversé dans l'étude du niveau 2. Il y a des questions telles que "le nombre de variables indépendantes + 1" réduit le degré de liberté du coefficient de régression de l'analyse de régression multiple, et "le coefficient n'est pas 0 au niveau de signification de xx%", mais l'histoire autour de cela est ~~cette fois. </ del> Il est lié au contenu à introduire.~~

De plus, nous n'utiliserons pas de bibliothèques statistiques ici, mais uniquement des bibliothèques telles que numpy et pandas (car les utiliser dans une boîte noire n'aide pas à comprendre ...).

Livre de référence

Katsuya Hasegawa, "Multivariate Analysis from Zero" (2004) Chapitre 2

Matière

Cette fois, nous traiterons et utiliserons les données de «Statut d'utilisation par route» de JR East. Les données suivantes, séparées par des virgules, ont été créées en formatant le nom de l'itinéraire et les kilomètres d'exploitation, le nombre moyen de passagers qui transitent (exercice 2018) et les revenus du transport de passagers (exercice 2018) dans le fichier PDF. ("Ofunato Line", "Kesennuma Line" et "Yamada Line" avec des données incomplètes sont exclues)

`jreast2018.csv`


Nom de la ligne,Km affaires,Nombre moyen de passants,Revenus du transport de passagers
Ligne Yamate,20.6,1134963,111014
Ligne Saikyo,5.5,752645,15032
Ligne principale de Tokaido,169.2,367765,238493
Ligne Yokohama,42.6,231706,38214
Ligne principale Sobu,145.4,207099,117700
Ligne Keiyo,54.3,181483,38201
Ligne Negishi,22.1,177099,16618
Ligne Nanbu,45,163148,29946
Ligne principale centrale,247.8,160328,170577
Ligne Takasaki,74.7,116646,34095
Ligne Musashino,105.5,114847,43177
Ligne principale de Tohoku,572,84567,183688
Ligne Joban,351,71035,100757
Ligne Oume,37.2,63151,8989
Ligne Yokosuka,23.9,62379,6107
Ligne Kawagoe,30.6,56191,6629
Ligne Sotobo,93.3,34851,12038
Ligne Sagami,33.3,29643,4372
Ligne Gokaichi,11.1,24712,1025
Ligne Sengoku,49,20497,4466
Ligne Uchibo,119.4,20483,8544
Ligne Ito,16.9,16907,1849
Nouvelle ligne blanche,27.3,15978,1686
Ligne Narita,119.1,15183,8901
Ligne Tsurumi,9.7,14133,569
Ligne Shinonoi,66.7,12465,4717
Les deux lignes de cheveux,84.4,11346,3176
Ligne Hachiko,92,9338,3612
Ligne principale Shinetsu,175.3,9233,6260
Ligne Senzan,58,9046,2176
Ligne Tokin,13.8,8075,373
Ligne Tazawako,75.6,7023,4421
Ligne Mito,50.2,7011,1255
Ligne Echigo,83.8,6021,1862
Lumière du soleil,40.5,5806,863
Ligne Joetsu,164.4,5389,3464
Ou Main Line,484.5,4983,14111
Ligne Sazawa,24.3,3342,244
Ligne Oito,70.1,3140,937
Ligne Yahiko,17.4,2338,150
Ligne Azuma,55.3,2327,615
Ligne principale Hagoshi,271.7,2194,2731
Ligne Oshika,26.4,1877,146
Ligne ouest d'Iwakoshi,175.6,1745,1198
Ligne Mizugun,147,1666,760
Ligne Karasuyama,20.4,1457,83
Ligne Iwakoshi Est,85.6,1385,371
Enroulement de pierre,44.7,1255,190
Ligne Kashima,17.4,1221,87
Ligne Koumi,78.9,1194,420
Ligne Kururi,32.2,1094,96
Ligne Rikuha Est,94.1,906,296
Ligne Hachinohe,64.9,883,269
Ligne Kamaishi,90.2,764,308
Ligne Gono,147.2,631,352
Ligne Iiyama,96.7,598,199
Ligne Ominato,58.4,578,168
Ligne Tsugaru,55.8,464,112
Ligne de guirlande,106.9,382,124
Ligne Yonesaka,90.7,379,105
Ligne Rikuha Nishi,43,345,67
Ligne Kitakami,61.1,311,70
Ligne Tadami,135.2,280,140

Enregistrez ce fichier sous ** jreast2018.csv avec ** le code de caractère UTF-8.

Entraine toi

Ci-dessous, nous avons confirmé l'opération avec Python 3.6.8.

Saisie des données

Tout d'abord, vous devez importer les données. Importez des données au format CSV à l'aide de pandas.

import numpy as np
import pandas as pd
df = pd.read_csv("jreast2018.csv")
print(df) #Vérifiez le contenu des données importées
Nom de la ligne Business km Nombre moyen de passagers de passage Revenus du transport de passagers
0 Ligne Yamate 20.6  1134963  111014
1 Saikyo Ligne 5.5   752645   15032
2 Ligne principale Tokaido 169.2   367765  238493
3 Ligne 42 de Yokohama.6   231706   38214
4 Ligne principale Sobu 145.4   207099  117700
..    ...    ...      ...     ...
58 ligne de guirlande 106.9      382     124
59 Ligne 90 Yonesaka.7      379     105
60 Rikuha Nishi Ligne 43.0      345      67
61 Kitakami Ligne 61.1      311      70
62 Tadami Ligne 135.2      280     140

Trouvez le coefficient de la droite de régression (section 2-1)

Tout d'abord, la variable indépendante (variable explicative) $ x $ est le "nombre moyen de passants" et la variable dépendante (variable objective) $ y $ le "revenu du transport de passagers", et tracez une ligne droite $ y = ax + b $ qui exprime bien la relation entre les deux. Je pense à tirer. Lorsque vous tracez une ligne droite, considérez une ligne droite qui minimise la somme des carrés des résidus de chaque point de données (la différence entre la valeur observée (valeur réalisée) de la variable dépendante et la valeur théorique (valeur prédite)). C'est la soi-disant «méthode du carré minimum». Pour chaque point de données $ (x_1, y_1), ..., (x_n, y_n) $, si la valeur observée de la variable dépendante est $ y_i $ et la valeur théorique est $ \ hat {y} _i = ax_i + b $ La somme des carrés des résidus $ L $ est

\begin{align}
L &= \sum_{i=1}^n (y_i - \hat{y}_i)^2 \\
&= \sum_{i=1}^n (y_i - ax_i - b)^2
\end{align}

Ce sera. Lorsque $ L $ prend la valeur minimale, le coefficient différentiel partiel sera de 0 pour $ a et b $ respectivement.

\begin{align}
\frac{\partial L}{\partial a} &= \sum_{i=1}^n 2(y_i - ax_i - b)(-x_i) = 0 \\
\frac{\partial L}{\partial b} &= \sum_{i=1}^n 2(y_i - ax_i - b)(-1) = 0 \\
\end{align}

Si vous organisez cela pour $ a et b $

\begin{align}
\left(\sum_{i=1}^n x_i^2\right) a + \left( \sum_{i=1}^n x_i \right) b &= \sum_{i=1}^n x_i y_i \\
\left( \sum_{i=1}^n x_i \right) a + nb &= \sum_{i=1}^n y_i \\
\end{align}

Lorsqu'il est résolu comme une équation simultanée pour $ a, b $

\begin{align}
a  &= \frac{n\sum_{i=1}^n x_i y_i - \sum_{i=1}^n x_i \sum_{i=1}^n y_i}{n \sum_{i=1}^n x_i^2 - \left( \sum_{i=1}^n x_i \right)^2} \\
 b &=  \frac{ \sum_{i=1}^n x_i^2 \sum_{i=1}^n y_i - \sum_{i=1}^n x_i \sum_{i=1}^n x_i y_i}{n \sum_{i=1}^n x_i^2 - \left( \sum_{i=1}^n x_i \right)^2} \\
\end{align}

Il peut être obtenu. Il est difficile de calculer ces $ a et b $ à la main (même si vous utilisez une calculatrice), mais vous pouvez facilement les calculer à l'aide d'un programme.

n = len(df)
X = df["Nombre moyen de passants"]
y = df["Revenus du transport de passagers"]
a = (n * (x*y).sum() - x.sum() * y.sum()) / (n * (x**2).sum() - x.sum() ** 2)
b = ((x**2).sum() * y.sum() - x.sum() * (x*y).sum()) / (n * (x**2).sum() - x.sum() ** 2)
print(a) # Output: 0.12665116273608049
print(b) # Output: 11411.585376160469

De là, on peut voir que les revenus annuels du transport de passagers ont tendance à augmenter d'environ 0,127 million de yens (= 127 000 yens) pour chaque augmentation du nombre moyen de passagers de passage. Dessinons également un graphique.

import matplotlib as mpl
mpl.use('Agg')
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.scatter(x, y, color="tab:blue")
ax.plot(x, a*x + b, color="tab:orange")
fig.savefig("output.png ")

Eh bien, je n'ai pas l'impression de très bien l'analyser. D'après le graphique, les données ne semblent pas très linéaires. Le point en bas à droite où la valeur de l'axe des x est d'environ 750 000 est la «ligne Saikyo». Si vous regardez attentivement les données originales, le nombre moyen de personnes qui transitent est important, mais comme les kilomètres d'exploitation sont aussi courts que 5,5 km [^ saikyo], les revenus du transport de passagers sont faibles, n'est-ce pas?

[^ saikyo]: La "Ligne Saikyo" ici ne fait pas référence à la "Ligne Saikyo" (Osaki-Omiya) dans les informations passagers, mais seulement à la section (Ikebukuro-Akaha) dont le nom officiel de l'itinéraire est la "Ligne Akaha". Est affiché dans la note de bas de page (* C) sur la page source.

Au fait, quel était le «nombre moyen de passants»? (À partir de "État d'utilisation par route")

○ Le «nombre moyen de passagers de passage» représente le nombre de passagers par 1km par jour, et est calculé par le calcul suivant basé sur le «Rapport sur la performance des entreprises ferroviaires» que nous rapportons chaque année au ministère du Territoire, des Infrastructures, des Transports et du Tourisme. Je suis. [Nombre moyen de passagers de passage] = [Kilométrage des transporteurs de passagers au cours de l'exercice pour chaque itinéraire * A] ÷ [Kilomètres d'affaires au cours de l'exercice pour l'itinéraire concerné] ÷ [Jours ouvrables de l'exercice]

Étant donné que ce sont les «kilomètres de transport de passagers» (nombre d'utilisateurs x distance) qui génèrent réellement des revenus, le «nombre moyen de passants x kilomètres d'affaires x nombre de jours ouvrables» (distance d'utilisation totale tout au long de l'année) est utilisé comme variable indépendante. Est considéré comme plus approprié. Il y a 365 jours en 2018, calculons donc le nombre de jours ouvrables à 365.

x = df["Nombre moyen de passants"] * df["Km affaires"] * 365
y = df["Revenus du transport de passagers"]
n = len(df)
a = (n * (x*y).sum() - x.sum() * y.sum()) / (n * (x**2).sum() - x.sum() ** 2)
b = ((x**2).sum() * y.sum() - x.sum() * (x*y).sum()) / (n * (x**2).sum() - x.sum() ** 2)
print(a) # Output: 1.082862145743171e-05
print(b) # Output: 331.32038547946183
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.scatter(x, y, color="tab:blue")
ax.plot(x, a*x + b, color="tab:orange")
fig.savefig("output.png ")

C'est un sentiment agréable. En supposant que le kilomètre passager annuel est de $ x $ et que les revenus annuels du transport de passagers (millions de yens) sont de y $, la valeur théorique de $ y $ $ \ hat {y} $ peut être calculée comme suit.

\hat{y} = (1.083 \times 10^{-5}) x + 331.3

Il est calculé que chaque fois que le chemin de fer est utilisé sur 1 km, le revenu augmentera d'environ 10,83 yens. Le tarif régulier de JR (3 compagnies à Honshu) au kilomètre (tarif par 1 km) est calculé à 16,20 yens (départ et arrivée mutuels sur la ligne principale, si elle est de 300 km ou moins, hors taxes), mais il est régulier [^ tarif] Si vous achetez un billet, le tarif sera réduit, de sorte que le revenu que vous gagnerez en utilisant 1 km sera inférieur à 16,20 yens. En outre, il existe en fait divers modèles de calcul des tarifs tels que de courtes distances de 10 km ou moins, de longues distances de plus de 300 km et des lignes de transport locales avec peu d'utilisateurs, et les données ne vont pas directement.

[^ fare]: [JR East: Passenger Business Rules> Volume 2 Passenger Business-Chapter 3 Passenger Fares / Prices-Section 2 Ordinary Passenger Fares](https://www.jreast.co.jp/ryokaku/02_hen / 03_syo / 02_setsu / index.html)

Trouvez le facteur de décision (section 2-2)

Maintenant, trouvons le facteur de décision $ R ^ 2 $ pour voir dans quelle mesure cette équation de régression s'adapte. Comme la fluctuation que la valeur observée $ y_i $ a

Variation totale: Fluctuation à partir de la moyenne de tous les points de données (somme des carrés de l'écart) --Variation de retour: Fluctuation à partir de la valeur théorique obtenue par l'équation de régression, c'est-à-dire variation expliquée par l'équation de régression (somme des carrés de l'écart)
Fluctuation résiduelle: Fluctuation non expliquée par l'équation de régression (somme des carrés de l'écart)

Compte tenu de la somme de tous les points de données, (variation totale) = (variation de retour) + (variation résiduelle) est valable. Ici, le rapport de la variation de régression à la variation totale est le coefficient de détermination $ R ^ 2 $. Plus le facteur de décision est grand (plus proche de 1), meilleure est l'équation de régression qui correspond aux données.

TSS = ((y - y.mean()) ** 2).sum()     #Fluctuation totale
ESS = ((a*x+b - y.mean()) ** 2).sum() #Fluctuation régressive
RSS = ((y - a*x+b) ** 2).sum()        #Fluctuation résiduelle
R_2 = ESS / TSS                       #Coefficient de décision
print(TSS) # 139066190955.65082
print(ESS) # 138362318781.54144
print(RSS) #    731535019.9636117
print(R_2) # 0.9949385816259694

Le facteur de décision est de 0,995, ce qui montre qu'il s'intègre bien. Je ne pense pas que ce soit si étrange parce que le calcul du tarif d'origine est basé sur une méthode proportionnelle à la distance. Cela devrait être (fluctuation totale) = (fluctuation de retour) + (fluctuation résiduelle), mais il semble qu'il y ait une erreur due à la perte d'informations parce que les données avec des nombres de chiffres différents ont été traitées dans un désordre.

Le coefficient de décision obtenu ici est en fait égal au carré du coefficient de corrélation de la variable indépendante et de la variable dépendante. Essayons de trouver et de confirmer le coefficient de corrélation (de l'échantillon).

cov_xy = ((x - x.mean()) * (y - y.mean())).sum() / (n - 1) #Covariance de l'échantillon
var_x  = ((x - x.mean()) ** 2).sum() / (n - 1)             #Exemple de distribution de variables indépendantes
var_y  = ((y - y.mean()) ** 2).sum() / (n - 1)             #Exemple de distribution des variables dépendantes
corr   = cov_xy / np.sqrt(var_x * var_y)                   #Coefficient de corrélation
print(corr ** 2) # 0.994938581625969

Certes, il semble correspondre au coefficient de décision obtenu précédemment (il n'est pas exactement égal en raison d'une erreur de calcul).

Parce que c'est devenu long

La prochaine fois. Considérons maintenant le cas où il y a deux variables indépendantes ou plus (multivariées), puis voyons comment gérer l'équation de régression de manière statistique. Suivant → [Test statistique 2e année / quasi 1e année] Formation à l'analyse de régression avec Python (2) --Qiita

[Test statistique 2e année / quasi 1e année] Formation à l'analyse régressive avec Python (1)