２ 統計量
本書第２章に関係する統計量、平均値・分散・標準偏差、歪度、尖度、四分位数、箱ひげ図、共分散と相関係数、回帰直線、偏相関係数、順位相関係数の計算Pythonスクリプト例を以下のように書いてみた。スクリプトファイルなどは圧縮ファイルchap2_files.zipにまとめた。ダウンロード・展開して自由に利用されたい。
Jupyter Notebook用にまとめたものも用意して、圧縮ファイルchap2_prgs.zipとした。こちらもダウンロード・展開して自由に利用されたい。
平均値、分散、標準偏差を定義に基づいて計算する場合と、Pythonの関数を用いる場合のコードを、リスト2.1のように書いてみた。
リスト2.1 平均値・分散・標準偏差の計算スクリプト例（s2_mn_sd_a.py）
import numpy as np
data = [5, 6, 6, 8]
v = 0.0
for t in data:
v += t
v_mean = v / len(data)
print('平均値 =', v_mean)
print('平均値 =', np.mean(data))
v = 0.0
for t in data:
v += (t - v_mean) ** 2
v_S2 = v / len(data)
v_s2 = v / (len(data) - 1)
print('分散 =', v_S2)
print('不偏分散 =', v_s2)
print(np.var(data, ddof = 0))
print(np.var(data, ddof = 1))
v_S = v_S2 ** 0.5
v_s = v_s2 ** 0.5
print('S =', v_S)
print('s =', v_s)
print(np.std(data, ddof = 0))
print(np.std(data, ddof = 1))
パラメータddofは、和を（データの個数―ddof）で割ることを示す。
リスト2.1を実行すると、以下の出力が得られる。
平均値 = 6.25
平均値 = 6.25
分散 = 1.1875
不偏分散 = 1.5833333333333333
1.1875
1.5833333333333333
S = 1.0897247358851685
s = 1.2583057392117916
1.0897247358851685
1.2583057392117916
歪度を本書（2.1.16）によって算出する場合とPythonの関数を用いる場合のスクリプト例をリスト2.2に示す。
リスト2.2 歪度の計算（s2_skewness.py）
import scipy.stats as ss
import numpy as np
data = [1]*10 +[2,3,4,5,6,7,8,9,10]
print(data)
data = np.array(data) # ndarray型に変換
m = np.mean(data)
S = np.std(data, ddof = 0)
g = np.mean((((data - m) / S) ** 3))
print(g)
print(ss.skew(data))
リスト2.2を実行すると、以下の出力を得る。
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3,
4, 5, 6, 7, 8, 9, 10]
0.9233534786348041
0.9233534786348037
浮動小数点の計算精度で同じ値である。
尖度を本書（2.1.17）式によって計算するコードと、Pythonの関数によって求める場合のスクリプト例をリスト２.3に示す。
リスト2.3 尖度の計算（s2_kurtosis.py）
import numpy as np
import scipy.stats as ss
data = [1,2,3,4,5,6,7,8,9,10] * 2
print(data)
data = np.array(data)
m = np.mean(data)
S = np.std(data, ddof = 0)
g2 = np.mean((((data - m) / S) ** 4))
print(g2)
print(ss.kurtosis(data, fisher =
False))
print(g2-3)
print(ss.kurtosis(data, fisher =
True))
パラメータfisherの値がFalseであれば、式（2.1.17）による値が、Trueであれば式（2.1.17）から３を引いた値が返される。
リスト2.3を実行すると以下の出力が得られる。
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1, 2,
3, 4, 5, 6, 7, 8, 9, 10]
1.7757575757575759
1.7757575757575756
-1.2242424242424241
-1.2242424242424244
データ「10, 20, 30, 40」の四分位数を求めるスクリプトをリスト2.4のように作成した。
リスト2.4 四分位数を求めるスクリプト（s2_percentiles.py）
import numpy as np
data = [10, 20, 30, 40]
print(data)
Qs = np.percentile(data, [25, 50, 75])
print(Qs)
リスト2.4を実行すると以下の出力を得る。
[10, 20, 30, 40]
[17.5 25. 32.5]
中央値が25となっているが、第１四分位数と第3四分位数は17.5と32.5で、それぞれ10と20、および30と40の中間ではない。
これらが中間の値になるように求める関数例が、岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版、p.11にMyPercentileとして載っている。これを用いたスクリプト例をリスト2.5に示す。
リスト2.5 ウエイトを用いて四分位数を求めるスクリプト（s2_percentiles_m.py）
import numpy as np
data = [10, 20, 30, 40]
print(data)
def MyPercentile( d, p ):
x = sorted(d)
n = len(x)
pos = n * (p / 100.0) + 0.5
if pos <= 1:
return x[0]
elif pos >= n:
return x[n - 1]
else:
i =
int(pos)
v =
(pos - i) * x[i] + (i + 1 - pos) * x[i - 1]
return v
Qs = [MyPercentile(data, p) for p in
[25, 50, 75]]
print(Qs)
リスト2.5を実行すると、以下の出力を得る。
[10, 20, 30, 40]
[15.0, 25.0, 35.0]
第１四分位数および第２四分位数が、10と20および30と40の中間地として算出されている。
Numpyのpercentile関数では、いろいろな計算法が用意されている。上の例に対してこれらの計算法による結果を比較するスクリプトをリスト2.5aのように作成してみた。
リスト2.5a Percentile計算法の比較（check_percentiles.py）
import numpy as np
a = [10, 20, 30, 40]
methods = ['inverted_cdf',
'averaged_inverted_cdf', 'closest_observation',
'interpolated_inverted_cdf', 'hazen', 'weibull',
'linear', # (default)
'median_unbiased', 'normal_unbiased', 'lower',
'higher', 'midpoint', 'nearest']
print('a =', a)
for m in methods:
Qs = np.percentile(a, [25,
50, 75], method = m)
print(f'{m:>25s}', Qs)
リスト2.5aを実行すると以下の結果を得る。
a = [10, 20, 30, 40]
inverted_cdf [10 20 30]
averaged_inverted_cdf [15.
25. 35.]
closest_observation [10 20 30]
interpolated_inverted_cdf [10. 20.
30.]
hazen [15. 25. 35.]
weibull
[12.5 25. 37.5]
linear [17.5 25. 32.5]
median_unbiased [14.16666667 25.
35.83333333]
normal_unbiased [14.375 25. 35.625]
lower [10 20 30]
higher [20 30 40]
midpoint [15. 25. 35.]
nearest [20 30 30]
本書では、箱ひげ図はヒンジに基づいて説明されている。しかし、Pythonで描かれる箱図は、ヒンジの代わりに四分位数が用いられている。本書の表1.1.1のデータの箱ひげ図を描くスクリプトをリスト2.6のように用意した。
リスト2.6 箱ひげ図描画スクリプト（s2_boxplot.py）
import matplotlib.pyplot as plt
data = [
66,
59, 62, 64, 63, 68, 65, 59, 68, 64,
65,
51, 67, 64, 83, 59, 61, 62, 57, 72,
65,
64, 54, 60, 53, 65, 67, 60, 53, 79,
74,
53, 61, 68, 75, 50, 57, 55, 66, 56,
55,
61, 70, 71, 49, 69, 70, 80, 73, 72]
# 日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP',
'Hiragino Maru Gothic Pro']
plt.boxplot(data, labels = ['テスト'])
plt.xticks(fontsize = 16)
plt.ylabel('点数', fontsize = 16)
plt.show()
リスト2.6を実行すると、図2.1の箱ひげ図が描画される。
図2.1
箱ひげ図の説明は、このウェブサイトでも行っている。
共分散を本書p. 29の式、相関件数をp. 30の式によって算出するスクリプトをリスト2.7のように用意した。
リスト2.7 共分散と相関係数の計算スクリプト（s2_fcov_cor.py）
import matplotlib.pyplot as plt
import numpy as np
data = [
(70,
168), (68, 165), (62, 159), (56, 156), (70, 170),
(64,
162), (59, 162), (70, 169), (68, 168), (69, 170),
(67,
168), (63, 166), (58, 163), (74, 170), (65, 168),
(60,
160), (65, 165), (72, 175), (61, 164), (77, 174),
(68,
164), (64, 161), (70, 166), (47, 147), (66, 165),
(68,
170), (73, 174), (74, 172), (63, 164), (83, 182),
(58,
159), (64, 166), (67, 167), (75, 174), (67, 163),
(68,
166), (68, 167), (66, 165), (63, 165), (55, 154),
(61,
163), (60, 160), (62, 164), (62, 160), (59, 160),
(51,
155), (59, 163), (68, 170), (56, 157), (79, 180)]
data = np.array(data)
print(data)
mx = np.mean(data.T[0])
my = np.mean(data.T[1])
Cxy = np.mean((data.T[0] - mx) *
(data.T[1] - my))
print('共分散 =', Cxy)
r = Cxy / \
((np.mean((data.T[0] - mx)
** 2) * np.mean((data.T[1] - my) ** 2)) ** 0.5)
print('相関係数 =', r)
リスト2.7を実行すると以下の出力を得る。
共分散 = 40.867999999999995
相関係数 = 0.9428052005583535
リスト2.7をPythonの関数を用いて、リスト2.8のように書き直した。
リスト2.8 Pythonの関数による共分散と相関係数の計算（s2_fcov_cor_by_func.py）
import matplotlib.pyplot as plt
import numpy as np
data = [
(70,
168), (68, 165), (62, 159), (56, 156), (70, 170),
(64,
162), (59, 162), (70, 169), (68, 168), (69, 170),
(67,
168), (63, 166), (58, 163), (74, 170), (65, 168),
(60,
160), (65, 165), (72, 175), (61, 164), (77, 174),
(68,
164), (64, 161), (70, 166), (47, 147), (66, 165),
(68,
170), (73, 174), (74, 172), (63, 164), (83, 182),
(58,
159), (64, 166), (67, 167), (75, 174), (67, 163),
(68,
166), (68, 167), (66, 165), (63, 165), (55, 154),
(61,
163), (60, 160), (62, 164), (62, 160), (59, 160),
(51,
155), (59, 163), (68, 170), (56, 157), (79, 180)]
data = np.array(data)
print(data)
Cxy = np.cov(data.T, ddof = 0)
print('共分散行列：\n', Cxy)
#print(np.cov(data.T, ddof = 0))
r = np.corrcoef(data.T)
print('相関行列：\n', r)
リスト2.8を実行すると、以下の出力を得る。
共分散行列：
[[46.8224 40.868 ]
[40.868 40.13 ]]
相関行列：
[[1.
0.9428052]
[0.9428052 1. ]]
関数の値は、それぞれ共分散行列および相関行列である。
回帰直線の係数を、本書（2.2.5）式により求めるスクリプトをリスト2.9のように作成した。
リスト2.9 単回帰直線を求めるスクリプト（s2_regression.py）
import matplotlib.pyplot as plt
import numpy as np
data = [
(70,
168), (68, 165), (62, 159), (56, 156), (70, 170),
(64,
162), (59, 162), (70, 169), (68, 168), (69, 170),
(67,
168), (63, 166), (58, 163), (74, 170), (65, 168),
(60,
160), (65, 165), (72, 175), (61, 164), (77, 174),
(68,
164), (64, 161), (70, 166), (47, 147), (66, 165),
(68,
170), (73, 174), (74, 172), (63, 164), (83, 182),
(58,
159), (64, 166), (67, 167), (75, 174), (67, 163),
(68,
166), (68, 167), (66, 165), (63, 165), (55, 154),
(61,
163), (60, 160), (62, 164), (62, 160), (59, 160),
(51,
155), (59, 163), (68, 170), (56, 157), (79, 180)]
data = np.array(data)
print(data)
N = len(data)
print('N =', N)
a = (N * np.sum(data.T[0] * data.T[1])
-
np.sum(data.T[0]) *
np.sum(data.T[1])) / \
(N *
np.sum(data.T[0]**2) - np.sum(data.T[0])**2)
b = (np.sum(data.T[0]**2) *
np.sum(data.T[1]) -
np.sum(data.T[0])
* np.sum(data.T[0] * data.T[1])) / \
(N *
np.sum(data.T[0]**2) - np.sum(data.T[0])**2)
print('a =', a)
print('b =', b)
# 日本語表記の設定
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP',
'Hiragino
Maru Gothic Pro']
plt.plot(data.T[0], data.T[1], 'o')
x_min = np.min(data.T[0])
x_max = np.max(data.T[0])
plt.plot([x_min, x_max], [a*x_min + b,
a*x_max + b],
label = 'a = {0:.3f}\nb = {1:.3f}'.format(a, b))
plt.title('回帰式 a*x + b', fontsize = 24)
plt.legend(fontsize = 12)
plt.show()
リスト2.9を実行すると図2.2のグラフが描かれる。
図2.2
図2.2を閉じると、端末には以下の出力が表示されている。
N = 50
a = 0.8728300984144342
b = 108.35656437944232
本書（2.2.13）式に基づいて偏相関係数を求めるスクリプトを、リスト2.10のように作成した。
リスト2.10 偏相関係数を求めるスクリプト（s2_partial_corr.py）
r_xy = float(input('r_xy = '))
r_xz = float(input('r_xz = '))
r_yz = float(input('r_yz = '))
r_xy_z = (r_xy - r_xz * r_yz) / \
(((1 - r_xz**2) * (1 - r_yz**2)) ** 0.5)
print('\nr_xy_z =', r_xy_z)
リスト2.10を実行すると、まず
の値の入力が求められる。
r_xy =
次のように、例えば「0.56」と設定する。
r_xy = 0.56
次に、
の値の設定が求められる。
r_xy = 0.56
r_xz =
値を、例えば「0.8」と設定する。
r_xy = 0.56
r_xz = 0.8
次に、
の値の設定が求められる。
r_xy = 0.56
r_xz = 0.8
r_yz = 0.7
例えば、「0.7」と設定して、Enterキーを押すと計算が始まり、以下のように偏相関係数
が表示される。
r_xy = 0.56
r_xz = 0.8
r_yz = 0.7
r_xy_z = 2.5910386503532426e-16
Spearmanの順位相関係数とKendallの相関係数をPearsonの相関係数とともに求めるスクリプトをリスト2.11のように用意してみた。
リスト2.11 順位相関係数を求めるスクリプト（r2_spearman_corr.py）
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
data = np.array([[1, 1], [2, 3], [5,
4], [6, 6]])
r_Pearson, p_v_pearson =
ss.pearsonr(data.T[0], data.T[1])
print('r_Pearson =', r_Pearson, ' p_value =', p_v_pearson)
r_Spearman, p_v_spearman =
ss.spearmanr(data)
print('r_Spearman =', r_Spearman,
' p_value =', p_v_spearman)
r_Kendall, p_v_kendall =
ss.kendalltau(data.T[0], data.T[1])
print('r_Kendall =', r_Kendall, ' p_value =', p_v_kendall)
plt.plot(data.T[0], data.T[1], 'o',
label = ('Pearson = {0:.3f}\n' +
'Spearman = {1:.3f}\nKendall = {2:.3f}').format(
r_Pearson, r_Spearman, r_Kendall))
plt.legend(fontsize = 12)
plt.show()
リスト2.11を実行すると、図2.3のグラフが表示される。
図2.3
図2.3のグラフを閉じると、端末には以下の出力が表示されている。
r_Pearson = 0.9417419115948377 p_value =
0.058258088405162334
r_Spearman = 1.0 p_value = 0.0
r_Kendall = 1.0 p_value = 0.08333333333333333
各相関係数の値と、そのｐ値が出力されている。