情報基礎「Pythonプログラミング」（ステップ７・統計計算・ファイル入力と相関係数）

Pythonプログラミング（ステップ７・統計計算・ファイル入力と相関係数）

このページでは、Pythonを使って基本的な統計量である相関係数を求める方法について考える。

ファイルからデータを読み込んでリストにセットする

コンピュータに統計処理を行わせる場合、プログラムの中に数値データを書き込んでしまうよりは、データを別ファイルに保存しておいて、プログラムはその内容を読んで、結果を出力する、といった流れのほうが何かと便利だ。例えば、１万件のデータを処理しなければならない場合、手入力で１件ずつデータを入力するやり方は明らかに非現実的である。そこで、ファイルに格納されたデータをPythonのプログラムに読み込む方法について確認しておこう。

ここでは、データファイルの形式として、データ１件分（１レコード）が１行ごとに、以下のように並んでいる場合を考えよう。

身長と体重データの例

ダウンロード

このデータが、height-vs-weight.txtというファイル名で、Pythonプログラムと同じディレクトリ（フォルダ）に保存されていると想定する。以下は、そこからデータを読み出して、値をリストにセットするプログラムの例である：

# coding: utf-8

file=open("height-vs-weight.txt")
lines=file.readlines()
file.close()

x_vs_y=[ ]
for i in range(0,len(lines),1):
    x,y=map(float,lines[i].split())
    print("i=",i," x=",x," y=",y)
    x_vs_y.append([x,y])

このプログラム中に登場した新しい記法や関数について、以下で説明する：

file=open("height-vs-weight.txt"): open("ファイル名")でファイルを開き、「ファイルについての情報」を変数fileにセットする。

lines=file.readlines(): ファイルからデータを読み込み、リストを作成し、linesにセットする。リストの0番目がファイルの1行目, 1番目が2行目・・・、にそれぞれ対応している。

file.close(): データを読み込んだので、ファイルを「閉じる。」

x_vs_y=[ ]: 数値データを保存するため、空のリストを用意する。

x,y=map(float,lines[i].split()): ファイルのi番目の行を実数値に変換して、それぞれ変数x,yにセットする。キーボード入力 x,y=map(float,input().split()) と同じパターン。

x_vs_y.append([x,y]): リストの末尾に [x,y] というリストを追加する。

以上の手順でデータを読み込むと、(0番目から始めて）i番目のデータのxはx_vs_y[i][0]、yはx_vs_y[i][1] で参照できる。身長(x)と体重(y)の平均をそれぞれ計算するコードの例を以下に示す。

# coding: utf-8

file=open("height-vs-weight.txt")
lines=file.readlines()
file.close()

x_vs_y=[ ]
for line in lines:
    x,y=map(float,line.split())
    x_vs_y.append([x,y])

n=len(x_vs_y)
s_x=0
s_y=0
for i in range(0,n,1):
    s_x = s_x + x_vs_y[i][0]
    s_y = s_y + x_vs_y[i][1]

mu_x=s_x/n
mu_y=s_y/n

print("身長の平均={0:3.1f} [cm]".format(mu_x))
print("体重の平均={0:2.1f} [kg]".format(mu_y))

データ間の関係性を調べる：ピアソンの相関係数

二つの変量$X$と$Y$の間の関係性を特徴づける量として、相関係数（さらに詳しくは、ここで述べるのはピアソンの相関係数）が広く使われている。その具体的な計算手順は以下のとおりである：

$n$ 個のデータのペア $(x_i,y_i)$ （ただし、$n=0,1,\cdots,n-1$）が与えられたとして、それぞれの平均 $$ E(X) = \frac{1}{n} \sum_{i=0}^{n-1} x_i,\ \ \ E(Y) = \frac{1}{n} \sum_{i=0}^{n-1} y_i $$ および、二乗平均 $$ E(XX) = \frac{1}{n} \sum_{i=0}^{n-1} {x_i}^2,\ \ \ E(YY) = \frac{1}{n} \sum_{i=0}^{n-1} {y_i}^2 $$ を求める。これらを用いると、$X$と$Y$の分散は $$ V_X = E(XX) - E(X)^2 ,\ \ \ V_Y = E(YY) - E(Y)^2 , $$ となる。

同様に、$X$と$Y$の積の平均 $$ E(XY) = \frac{1}{n} \sum_{i=0}^{n-1} x_i y_i $$ から、共分散(covariance) $$ V_{XY} = E(XY) - E(X) E(Y) $$ が得られる。共分散を分散で規格化した量がピアソンの相関係数で、 $$ r = \frac{V_{XY}}{\sqrt{V_X V_Y}} $$ で与えられる。相関係数は $-1 \le r \le +1$ の値を取る。

横軸と縦軸のスケールをデータのばらつきに合わせ調整した上で、$X$, $Y$を散布図で表現した際に、点が斜め上45度方向の直線に沿って分布すると$r$は1に近く、右下がり45度方向の直線に沿って分布すると$r$は-1に近い値を取る。反対に、$r$が0に近い場合は、$X$が増えても$Y$は同じようには（直線的には）増えない、という状況に対応しているが、このことは、必ずしも$X$と$Y$の間に関係性がないことを意味するわけではない。

　練習：相関係数

プロ野球球団の選手の身長(cm)、体重(kg)、年齢、推定年俸(万円）などのデータがプロ野球データFreakで公開されている。それを元に、2015年のある球団について、

データファイル: rakuten.txt

身長(cm)  体重(kg)  年齢  推定年俸（万円）

の数値を並べたファイルをこちらに用意した。このファイル（rakuten.txt）を読み込み、 W:身長 X:体重 Y:年齢 Z:年俸としたとき、年俸との相関係数（$R_{WZ}, R_{XZ}, R_{YZ}$）を計算するプログラムを作成せよ。

なお、年俸は、他の変量と較べて、選手（働き）によって桁が違ってくるような量である。すなわち、金額そのものよりも、その「桁」のほうに注目するのが自然であろう。そのような観点から、年俸については、math.log10()関数を使い、金額の対数を変量とした解析を行うこと。

年俸と最も相関が強いと考えられるのは、身長か、体重か、年齢か？

　ヒント

データを4件ずつ読み込むためには、このページで示したサンプルプログラムの中で

import math

file=open("rakuten.txt")
lines=file.readlines()
file.close()

wxyz=[ ]
for i in range(0,len(lines),1):
    w,x,y,z=map(float,lines[i].split())
    wxyz.append([w,x,y,math.log10(z)])

といった書き換えが必要となる。上記のコードでは、4番目の変量については、元データの常用対数(math.log10())を取っている。

　解説: 疑似相関

$X$と$Y$の相関係数が大きいからと言って、$X$と$Y$の間に直接的な因果関係があるとは限らない。たとえば、Spurious Correlationsのサイトには、疑似相関の例がいくつも紹介されている。

例えば、第三の（隠れた）要因$Z$があって、$X, Y$が共に$Z$に影響されているような場合、$X$と$Y$の間に強い相関が見られる可能性がある。

データの分布形を仮定できない場合：スピアマンの順位相関

取り得る値が制限されていたり、桁違いに大きな（小さな）値が含まれるようなデータの相関関係をみる場合、データを一旦変換しておいて、変換されたほうの値で考えるのが自然な場合も多い。その中でも、「順位」は、元データの確率密度関数には依らない、一般的な指標となる。例えば、試験の点数を考えてみると、試験の内容毎に点数の分布は異なるが、点数の順位は、受験者の実力レベルを（安定して）反映しているはずである。（ここでは当然、全員が満点となったりしないような、弁別性の高い試験を想定している）。

$n$件の対となるデータ $(X_i, Y_i)$ の組が与えられており、$\{X_j \,|\, j=0,\cdots,n-1 \}$ の中での $X_i$の順位が$r_X(X_i)$、$\{Y_j \,|\, j=0,\cdots,n-1\}$ の中での $Y_i$の順位が $r_Y(Y_i)$ であったとしよう。

式が煩雑になるので、以下では、カテゴリー別の順位を $x_i = r_X(X_i)$, $y_i = r_Y(Y_i)$ と表すことにする。

これら順位に対して、ピアソンの相関係数と全く同様に、共分散$V_{XY}$、分散$V_X$, $V_Y$ から計算される $$ r_s = \frac{V_{XY}}{\sqrt{V_X V_Y}} $$ を、スピアマンの順位相関係数と呼ぶ（定義から $-1 \le r_s \le +1$）。

ここでは簡単のため、同順位の(tied）データは含まれないものとする。

ここで、順位は$1$から$n$までの値を取るとすると $$ E(X) = E(Y) = \frac{1}{n} \sum_{k=1}^{n} k = \frac{n+1}{2} $$ $$ E(XX) = E(YY) = \frac{1}{n} \sum_{k=1}^{n} k^2 = \frac{(n+1)(2n+1)}{6} $$ $$ V_X = V_Y = \frac{(n+1)(2n+1)}{6} -\left( \frac{n+1}{2} \right)^2 = \frac{n^2 - 1}{12} $$ $$ V_{XY} = \frac{1}{n} \sum_{i=0}^{n-1} x_i y_i - \left( \frac{n+1}{2} \right)^2 \\ $$ ここで、共分散 $V_{XY}$ に $$ \frac{1}{n} \sum_{i=0}^{n-1} x_i y_i = \frac{(n+1)(2n+1)}{6} - \frac{1}{2n} \sum_{i=0}^{n-1} (x_i - y_i)^2 $$ を代入すると、 $$ V_{XY} = \sqrt{V_X V_Y} - \frac{1}{2n} \sum_{i=0}^{n-1} (x_i - y_i)^2 $$ となるから、順位相関 $r_s$ は、順位の差の自乗和 $\sum_i (x_i - y_i)^2$ を使って $$ r_s = \frac{V_{XY}}{\sqrt{V_X V_Y}} = 1 - \frac{6}{n^3-n} \sum_{i=0}^{n-1} (x_i - y_i)^2 $$ で計算できる。

　解説:データのランキングの取得

リスト中のデータの順位は、SciPyライブラリを使って

from scipy.stats import rankdata

x = [3.1, 1.2, 6.1, 9.8, 2.8]
rank = rankdata(x)

のように記述することによって得ることができる（小さいデータから順に順位づけられる）。特に指定しない場合、同じ値のデータが含まれる場合、順位の平均値が返される（２番めに大きいデータが２つあれば、順位 2.5、等）。詳しくはscipy.stats.rankdataのマニュアルを参照のこと。

　練習：順位相関

前節の課題で分析したプロ野球選手のデータについて、スピアマンの順位相関係数を計算するプログラムを作成してみなさい。相関を計算するペアは、年俸と身長、年俸と体重、年俸と年齢とする。