統計量
平均値、分散などの１変量統計量、中央値、四分位数などの順序統計量、相関係数、共分散などの２変量統計量、偏相関係数を求めるPythonスクリプト例を以下のように用意した。
Pythonスクリプトファイルは圧縮ファイルScriptStat.zipとした。圧縮ファイルScriptStat.zipには、Jupyter Notebook用のファイルも含まれている。ダウンロード・解凍して自由に利用されたい。
なお、Excelなどの表計算ソフトで作成されるCSV形式のファイルから読み込む例を、本ウェブサイトの後半に用意した。CSVファイルの作成と読み込みの簡単な説明のウェブサイトも用意した。
日本気象協会のウェブサイトで横浜市の10日間最高気温予測を調べると、表A.1のようであった。
|
表A.1 横浜市の10日間最高気温予測より（日本気象協会ウェブサイト：2019.10.28） |
||||||||||
|
日 |
10月29日 |
30日 |
31日 |
11月1日 |
2日 |
3日 |
4日 |
5日 |
6日 |
7日 |
|
温度 |
16℃ |
23℃ |
22℃ |
24℃ |
22℃ |
20℃ |
18℃ |
19℃ |
20℃ |
19℃ |
平均値、分散、標準偏差を求めるスクリプトを次のように作成した。
import numpy as np
data = [16, 23, 22, 24, 22, 20, 18,
19, 20, 19]
m = np.mean(data)
v = np.var(data)
uv = np.var(data, ddof=1)
sd = np.std(data)
ssd = np.std(data, ddof=1)
print('平均値
= ', m)
print('分散 = ', v)
print('不偏分散 = ', uv)
print('標準偏差 = ', sd)
print('サンプル標準偏差 = ', ssd)
実行すると、以下の出力を得る。
平均値 = 20.3
分散 = 5.41
不偏分散 = 6.011111111111111
標準偏差 = 2.3259406699226015
サンプル標準偏差 = 2.4517567397911058
上の表A.1のデータの中央値、四分位数、最小値、最大値を得るスクリプトを、次のように作成した。
import numpy as np
# リストをnumpyの配列ndarray型として用意する
data = np.array([16, 23, 22, 24, 22,
20, 18, 19, 20, 19])
data_ascend = np.sort(data)
idx_d = np.argsort(-data) # dataの値の符号を反転してから昇順に添字を並べる
data_descend = data[idx_d]
v_min = np.amin(data)
q1 = np.percentile(data, 25)
med = np.median(data)
q3 = np.percentile(data, 75)
v_max = np.amax(data)
print('データ = ', data)
print('順序統計量
= ', data_ascend)
print('降順
= ', data_descend)
print('最小値 = ',
v_min)
print('第１四分位数 = ', q1)
print('中央値 = ', med)
print('第３四分位数 = ', q3)
print('最大値 = ',
v_max)
上のスクリプトを実行すると、以下の出力を得る。
データ = [16 23 22 24 22 20 18 19 20 19]
順序統計量 = [16 18 19 19 20 20 22 22 23 24]
降順
= [24 23 22 22 20 20 19 19
18 16]
最小値 = 16
第１四分位数 = 19.0
中央値 = 20.0
第３四分位数 = 22.0
最大値 = 24
まず、並べ換えられたものが順序統計量として出力され、その後、最小値、第１四分位数、中央値、第３四分位数、最大値が出力されている。
2教科、Sub. AとSub. B、のスコアが表F.1のようであるとする。
|
表F.1 2教科[Sub.A,
Sub.B]のスコア |
|
[79, 72]
[77, 70] [74,
75] [53, 44] [76, 52] [68, 60]
[96, 96] [85, 71] [74, 49] [93, 80] [54, 55]
[86, 70] [89,
96] [87, 87] [81, 81] [84, 85]
[78, 64] [76,
53] [54, 57] [75, 60] |
表F.1の２教科のスコアの相関係数、共分散を求めるプログラムを、次のように用意した。
import numpy as np
import scipy.stats as ss
data = [
[79,
72], [77, 70], [74, 75], [53, 44], [76, 52],
[68,
60], [96, 96], [85, 71], [74, 49], [93, 80],
[54,
55], [86, 70], [89, 96], [87, 87], [81, 81],
[84,
85], [78, 64], [76, 53], [54, 57], [75, 60]
]
A = []
B = []
for v in data:
A.append(v[0])
B.append(v[1])
r = np.corrcoef(A, B) # ピアソンの積率相関係数
print(r)
print('相関係数 = ', r[0][1])
C = np.cov(A, B)
print(C)
print('推定共分散 = ', C[0][1])
n = len(data)
C1 = C * (n-1) / n
print('標本共分散 = ', C1[0][1])
sprmn_r, p_spr = ss.spearmanr(A,
B) # スピアマンの順位相関係数
print('Spearman = ',sprmn_r, ' p_value = {0:.5f}'.format(p_spr))
kndl_tau, p_kndl = ss.kendalltau(A,
B) # ケンドールの順位相関係数
print("Kendall's tau = ",
kndl_tau, ' p-value =
{0:.5f}'.format(p_kndl))
リストAとBに対応する順番に値を設定して、相関行列を算出する関数を
r = np.corrcoef(A, B)
と呼び出している。この関数は、相関行列を返すので、相関係数は行列の成分r[0][1]の値として得られる。
共分散も同様であるが、次の関数呼び出し
C = np.cov(A, B)
では、共分散行列は分母を「データ数―１」で割ったもの、したがって分散は不偏分散が算出されるので、次のように
n = len(data)
C1 = C * (n-1) / n
標本数で調整すれば標本分散が得られる。
上のスクリプトを実行すると、以下の出力が得られる。
[[1.
0.77718556]
[0.77718556 1. ]]
相関係数 = 0.7771855633037875
[[148.15526316 144.78157895]
[144.78157895 234.23947368]]
推定共分散 = 144.7815789473684
標本共分散 = 137.54249999999996
Spearman = 0.8131122833458929 p_value = 0.00001
Kendall's tau = 0.6524064171122994 p-value = 0.00007
２変量
と
の相関係数
から変量
の影響を単回帰モデルで除いた相関係数、すなわち偏相関係数
は次式で与えられる。
ここで、
はとの相関係数、
はとの相関係数である。
偏相関係数を求めるスクリプトを次のように用意した。
Rxy = float(input('r_xy = '))
Rxz = float(input('r_xz = '))
Ryz = float(input('r_yz = '))
Rxy_z = (Rxy - Rxz * Ryz) / (((1 -
Rxz**2) * (1 - Ryz**2)) ** 0.5)
print('\nr_xy.z = ', Rxy_z)
上のスクリプトを実行すると、相関係数、、の値を変数r_xy、r_xz、r_yzの値として入力が求められる。３つの相関係数を入力すると、偏相関係数の値がr_xy.zの値として出力される。実行例を以下に示す。
r_xy = 0.45
r_xz = 0.8
r_yz = 0.75
r_xy.z = -0.3779644730092275
CSV形式のファイル（表計算ソフトExcelなどで作成）から読み込む場合
CSVファイルとして用意されたデータを読み込むサンプルスクリプトを用意した。CSVファイルは、Excelなどの表計算ソフトで作成したシート上のデータを、ファイル拡張子「.csv」を指定して保存すると、CSV形式のファイルとして保存される。１変量の統計量、２変量の統計量、偏相関係数の計算例を用意したが、ファイルにおける変量の数は任意である。例えば、１変量の統計量の計算スクリプトにおいて、データは多変量データのファイルであってもよく、スクリプトの実行時に統計量を計算する変数を１つ選ぶ。
スクリプトファイルなどは、圧縮ファイルForCSV.zipにまとめた。圧縮ファイルForCSV.zipには、Jupyter Notebook用のファイルも含まれている。
CSVファイルからデータを読み込み、１変量統計量の計算を行うスクリプトをリストCSV1に示す。
リストCSV1のスクリプトを実行すると、入力データファイル名の設定が求められる。
File Name
(*.csv) =
CSVファイルは、図CSV1に示す形式で用意する。
図CSV1 横浜市の10日間最高気温予測より（日本気象協会ウェブサイト：2019.10.28）
１行目に変数名を書き、その下に値を書いていく。図CSV１では、変数は日と気温の2つであるが、３つ以上であってもよい。
データを記入し終わったら、ファイル名の拡張子を「.csv」と指定して保存する。
次の例では、入力ファイルとして図ＣＳＶ１のファイルOneVar_jp.csvを指定している。
File Name
(*.csv) = OneVar_jp.csv
ファイル名を設定すると、データが読み込まれ、変数名のリストが表示される。リストの中から、統計量を計算する変数名の指定が求められる。
File Name (*.csv)
= OneVar_jp.csv
['日' '気温']
Variable name
=
変数名のリストは、[と]の間に、変数名を’と’で挟んだ形で表示される。変数名の設定は、’と’を除いて設定する。
Variable name
= 気温
変数名を設定すると、計算が始まる。
計算が終了すると、ヒストグラムが表示される（図CSV2）。
図CSV2
図CSV2のWindowを閉じると、端末には以下の出力が表示されている。
平均値 = 20.3
分散 = 5.41
不偏分散 = 6.011111111111111
標準偏差 = 2.3259406699226015
サンプル標準偏差 = 2.4517567397911058
データ = [16 23 22 24 22 20 18 19 20 19]
順序統計量 = [16 18 19 19 20 20 22 22 23 24]
降順
= [24 23 22 22 20 20 19 19
18 16]
最小値 = 16
第１四分位数 = 19.0
中央値 = 20.0
第３四分位数 = 22.0
最大値 = 24
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
in_fn = input('File Name (*.csv) = ')
d_frame = pd.read_csv(in_fn,
encoding='shift_jis')
print(d_frame.keys().values)
var_name = input('Variable name = ')
data = d_frame[var_name].values
m = np.mean(data)
v = np.var(data)
uv = np.var(data, ddof=1)
sd = np.std(data)
ssd = np.std(data, ddof=1)
print('平均値
= ', m)
print('分散 = ', v)
print('不偏分散 = ', uv)
print('標準偏差 = ', sd)
print('サンプル標準偏差 = ', ssd)
data_ascend = np.sort(data)
idx_d = np.argsort(-data) # dataの値の符号を反転してから昇順に添字を並べる
data_descend = data[idx_d]
v_min = np.amin(data)
q1 = np.percentile(data, 25)
med = np.median(data)
q3 = np.percentile(data, 75)
v_max = np.amax(data)
print('データ = ', data)
print('順序統計量
= ', data_ascend)
print('降順
= ', data_descend)
print('最小値 = ',
v_min)
print('第１四分位数 = ', q1)
print('中央値 = ', med)
print('第３四分位数 = ', q3)
print('最大値 = ',
v_max)
plt.hist(data)
plt.xlabel(var_name, fontsize=20)
plt.ylabel('度数', fontsize=20)
plt.tight_layout()
plt.show()
CSVファイルからデータを読み込んで、２変量統計量の計算を行うスクリプトをリストCSV２に示す。
リストCSV２のスクリプトを実行すると、入力データファイル名の設定が求められる。
File Name
(*.csv) =
CSVファイルは、図CSV３に示す形式で用意する。
図CSV３ 仮想データ
１行目に変数名を書き、その下に値を書いていく。図CSV２では、変数はデータ番号、科目Aと科目Bの３つであるが、４つ以上であってもよい。
データを記入し終わったら、ファイル名の拡張子を「.csv」と指定して保存する。
次の例では、入力ファイルとして図ＣＳＶ３のファイルTwoVars_jp.csvを指定している。
File Name (*.csv) = TwoVars_jp.csv
ファイル名を設定すると、データが読み込まれ、変数名のリストが表示される。リストの中から、統計量を計算する変数名の指定が求められる。
File Name
(*.csv) = TwoVars_jp.csv
['データ番号' '科目A' '科目B']
Variable
Name_1 =
変数名のリストは、[と]の間に、変数名を’と’で挟んだ形で表示される。変数名の設定は、’と’を除いて設定する。
まず、第１変数を設定する。
Variable
Name_1 = 科目A
Variable
Name_2 =
第１変数の名前を設定すると、第２変数の名前を聞いてくる。
Variable Name_1 = 科目A
Variable Name_2 = 科目B
第２変数の名前を設定すると、計算が始まる。
計算が終わると、散布図が表示される（図CSV４）。
図CSV４
図CSV４のWindowを閉じると、端末には以下の出力が表示されている。
相関係数 = 0.7771855633037875
[[148.15526316 144.78157895]
[144.78157895 234.23947368]]
推定共分散 = 144.7815789473684
標本共分散 = 137.54249999999996
Spearman = 0.8131122833458929 p_value = 0.00001
Kendall's tau = 0.6524064171122994 p-value = 0.00007
import numpy as np
import scipy.stats as ss
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
in_fn = input('File Name (*.csv) = ')
d_frame = pd.read_csv(in_fn,
encoding='shift_jis')
print(d_frame.keys().values)
nm_A = input('Variable Name_1 = ')
nm_B = input('Variable Name_2 = ')
A = d_frame[nm_A].values
B = d_frame[nm_B].values
r = np.corrcoef(A, B) # ピアソンの積率相関係数
print('相関係数 = ', r[0][1])
C = np.cov(A, B)
print(C)
print('推定共分散 = ', C[0][1])
n = len(A)
C1 = C * (n-1) / n
print('標本共分散 = ', C1[0][1])
sprmn_r, p_spr = ss.spearmanr(A,
B) # スピアマンの順位相関係数
print('Spearman = ',sprmn_r, f' p_value = {p_spr:.5f}')
kndl_tau, p_kndl = ss.kendalltau(A,
B) # ケンドールの順位相関係数
print("Kendall's tau = ",
kndl_tau, f' p-value =
{p_kndl:.5f}')
plt.scatter(A, B)
plt.title(f'r =
{r[0][1]:.3f}'.format(r), fontsize=20)
plt.xlabel(nm_A, fontsize=16)
plt.ylabel(nm_B, fontsize=16)
plt.tight_layout()
plt.show()
CSVファイルからデータを読み込んで、偏相関係数（partial coefficient correlation）の計算を行うスクリプトをリストCSV３に示す。
リストCSV３のスクリプトを実行すると、入力データファイル名の設定が求められる。
File Name
(*.csv) =
CSVファイルは、図CSV５に示す形式で用意する。
図CSV５ 仮想データ（リストCSV４のスクリプトで作成）
１行目に変数名を書き、その下に値を書いていく。図CSV５では、変数はデータN、科目X、科目Yと科目Zの４つであるが、５個以上であってもよい。
データを記入し終わったら、ファイル名の拡張子を「.csv」と指定して保存する。
次の例では、入力ファイルとして図ＣＳＶ５のファイルThreeVars_jp.csvを指定している。
File Name (*.csv) = ThreeVars_jp.csv
ファイル名を設定すると、データが読み込まれ、変数名のリストが表示される。リストの中から、統計量を計算する変数名の指定が求められる。
File Name
(*.csv) = ThreeVars_jp.csv
['データN' '科目X' '科目Y' '科目Z']
Variable
Name_X =
変数名のリストは、[と]の間に、変数名を’と’で挟んだ形で表示される。変数名の設定は、’と’を除いて設定する。
まず、第１変数を設定する。
Variable
Name_X = 科目X
Variable Name_Y
=
第１変数の名前を設定すると、第２変数の名前を、第２変数の名前を設定すると第３変数の名前の指定が求められる。
Variable Name_X = 科目X
Variable Name_Y = 科目Y
Variable Name_Z = 科目Z
第３変数の名前を設定すると、計算が始まる。
計算が終了すると、散布図が表示される（図CSV６）。
図CSV６
図CSV６のWindowを閉じると、端末には以下の出力が表示されている。
r_xy.z = -0.48095159236812796
Rxy_z_reg = -0.48095159236812685
偏相関係数r_xy.zは、公式（岡本、2022、p.41）によるものであり、Rxy_z_regは定義に基づき残差の相関係数として算出されたものである。
import numpy as np
import scipy.stats as ss
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
in_fn = input('File Name (*.csv) = ')
d_frame = pd.read_csv(in_fn,
encoding='shift_jis')
print(d_frame.keys().values)
nm_X = input('Variable Name_X = ')
nm_Y = input('Variable Name_Y = ')
nm_Z = input('Variable Name_Z = ')
X = d_frame[nm_X].values
Y = d_frame[nm_Y].values
Z = d_frame[nm_Z].values
XYZ = np.vstack((X, Y, Z))
Rs = np.corrcoef(XYZ) # Rsは相関行列を表す
Rxy = Rs[0][1]
Rxz = Rs[0][2]
Ryz = Rs[1][2]
#
# 偏相関係数の公式による計算（岡本、2022、p.41）
#
Rxy_z = (Rxy - Rxz * Ryz) / (((1 -
Rxz**2) * (1 - Ryz**2)) ** 0.5)
print('\nr_xy.z = ', Rxy_z)
n = len(X)
Axz = np.hstack((np.ones((n,1)),
np.reshape(Z,(n,1))))
Bxz = np.linalg.pinv(Axz) @
np.reshape(X, (n,1))
X_Z = X - np.reshape(Axz @ Bxz, n)
Ayz = np.hstack((np.ones((n,1)),
np.reshape(Z,(n,1))))
Byz = np.linalg.pinv(Ayz) @
np.reshape(Y, (n,1))
Y_Z = Y - np.reshape(Ayz @ Byz, n)
#
# 変数Ｚのリニアな影響を除いた残差による偏相関係数の計算（定義）
#
Rxy_z_reg = np.corrcoef(X_Z,
Y_Z)[0][1]
print('Rxy_z_reg =', Rxy_z_reg)
plt.figure(figsize=(7,7)) # ４枚の散布図を描く画面
plt.subplot(221)
plt.scatter(X, Y)
plt.title(f'Cor({nm_X},{nm_Y}) =
{Rxy:.3f}', fontsize=14)
plt.xlabel(nm_X, fontsize=14)
plt.ylabel(nm_Y, fontsize=14)
plt.subplot(222)
plt.scatter(X_Z, Y_Z)
plt.title(f'P.Cor({nm_X},{nm_Y}|{nm_Z})
= {Rxy_z:.3f}', fontsize=12)
plt.xlabel(f'{nm_X}|{nm_Z}',
fontsize=14)
plt.ylabel(f'{nm_Y}|{nm_Z}',
fontsize=14)
plt.subplot(223)
plt.scatter(Z, X)
plt.title(f'Cor({nm_Z},{nm_X}) =
{Rxz:.3f}', fontsize=14)
plt.xlabel(nm_Z, fontsize=14)
plt.ylabel(nm_X, fontsize=14)
plt.subplot(224)
plt.scatter(Z, Y)
plt.title(f'Cor({nm_Z},{nm_Y}) =
{Ryz:.3f}', fontsize=14)
plt.xlabel(nm_Z, fontsize=14)
plt.ylabel(nm_Y, fontsize=14)
plt.tight_layout()
plt.show()
import numpy as np
import matplotlib.pyplot as plt
rngen = np.random.default_rng(12345)
xy = rngen.multivariate_normal([50,
50],
[[10**2, -0.5*10**2],
[-0.5*10**2, 10**2]],
size=100)
z = rngen.normal(50, 10, size=100)
x = xy.T[0] + z - 50
y = xy.T[1] + z - 50
with open('ThreeVars_jp.csv', 'w') as
f:
f.write('"DataN","科目X","科目Y","科目Z"')
for i, v in enumerate(zip(x,
y,z)):
f.write(f'\n{i+1},{int(v[0])},{int(v[1])},{int(v[2])}')
plt.figure(figsize=(7, 7))
plt.subplot(221)
plt.scatter(x, y)
plt.xlabel('x', fontsize=18)
plt.ylabel('y', fontsize=18)
plt.subplot(222)
plt.scatter(x, z)
plt.xlabel('x', fontsize=18)
plt.ylabel('z', fontsize=18)
plt.subplot(223)
plt.scatter(y,z)
plt.xlabel('y', fontsize=18)
plt.ylabel('z', fontsize=18)
plt.tight_layout()
plt.show()