統計

pythonによる統計学演算実習

準備

pythonでは、あらかじめnumpy, scipyモジュール等をインポートします。

import math
import numpy as np
from scipy import stats

pythonのコマンドラインで実行していくこともできますが、jupyter notebookを利用するとグラフ表示ができ、プログラム修正も簡単に行うことができます。

sudo pip install jupyter # インストールされていない場合はインストール
jupyter notebook & # バックグラウンドで実行

jupyter notebookはTCP9999番ポートを利用しますので、ネットワーク環境に応じて該当ポートを開放してください。

各言語毎の統計量に関する関数

一変量解析

Excel R Python
合計 SUM(値1,値2,…) sum(x) sum(x)
算術平均 AVERAGE(値1,値2,…) mean(x) np.mean(x)
幾何平均 GEOMEAN(値1,値2,…) prod(x)^(1/length(x)) stats.gmean(x)
調和平均 HARMEAN(値1,値2,…) 1/sum(1/x)*length(x) stats.hmean(x)
中央値 MEDIAN(値1,値2,…) median(x) np.median(x)
最頻値 MODE(値1,値2,…)  - stats.mode(x)
偏差平方和 DEVSQ(値1,値2,…) sum((x-mean(x))^2)
標本分散
(Nで割る)
VARP(値1,値2,…) var(x)*(length(x)-1)/length(x) np.var(x)
不偏分散
(N-1で割る)
VAR(値1,値2,…) var(x) np.var(x, ddof=1)
標準偏差
(Nで割る)
STDEVP(値1,値2,…)  - np.std(x)
標準偏差
(N-1で割る)
STDEV(値1,値2,…) sd(x) np.std(x, ddof=1)
最大値 MAX(値1,値2,…) max(x) max(x)
最小値 MIN(値1,値2,…) min(x) min(x)
レンジ MAX-MIN range(x)
パーセンタイル PARCENTILE(配列,戻り値)
四分位点 QUARTILE(配列,戻り値) quantile(x) np.percentile(x,value)
四分位偏差 IQR(x)
5数要約 fivenum(x)
要約統計量 summary(x)
尖度 KURT(値1,値2,…) stats.kurtosis(x)
歪度 SKEW(値1,値2,…) stats.skew(x)
度数 FREQUENCY(データ配列,区間配列)

二変量解析

Excel R Python
相関係数 CORREL(配列1,配列2) np.corrcoef(x, y)
単回帰式の傾きa SLOPE(目的変数,説明変数)
単回帰式の定数項b INTERCEPT(目的変数,説明変数)
決定係数 RSQ(目的変数,説明変数)
相関比
連関係数

確率分布

Excel R Python
二項分布 BINOMDIST(X,n,p,定数)
X:0からnまでの整数で確率Pの事象が起こる回数
n:試行回数
p:ある事象が起こる確率
定数:0または1
(0はP(x)、1はP(0)からP(X)までの累積確率)
ポアソン分布 POISSON(X,m,定数)
X:発生回数
m:平均
定数:0または1
正規分布 NORMDIST(X,n,p,定数)
x:特定の値
m:平均値
s:標準偏差
定数:0または1
標準正規分布の累積確率 NORMSDIST(z)
z:基準値
標準正規分布の基準値 NORMSINV(p)
p:累積確率
t分布における確率
(母平均の検定のp値をとる確率)
TDIST(X,自由度,尾部)
X:検定統計量t
自由度:n-1
定数2
t分布におけるt値 TINV(確率,自由度)
カイ2乗分布における確率 CHIDIST(X,自由度)
X:検定統計量T
自由度:n-1
カイ2乗分布におけるX2値 CHINV(確率,自由度)

一変量解析の演習

基本統計量

【例題】ある会社のセールスマンの1ヶ月間の商品販売台数

名前 販売台数 偏差 偏差の平方
A 5 0 0
B 3 -2 4
C 4 -1 1
D 7 2 4
E 6 1 1
合計 25 0 10
平均 5
分散 2.5
標準偏差 1.6
■標準関数、標準モジュールでの計算

各データを変数xに代入し、各統計量を表示してみます。

x = [5, 3, 4, 7, 6]

print('平均値:', np.mean(x))
print('中央値:', np.median(x))
print('四分位点(25%):', np.percentile(x, 25))
print('四分位点(50%):', np.percentile(x, 50))
print('四分位点(75%):', np.percentile(x, 75))
print('最頻値:', stats.mode(x))
print('合計:', sum(x))
ave =  np.mean(x)
print('偏差:', [d - ave for d in x])
sq = [(d - ave) ** 2 for d in x]
print('偏差の平方:', sq)
print('偏差平方和:', sum(sq))
print('標本分散:', np.var(x))
print('不偏分散:', np.var(x, ddof=1))
print('標準偏差(Nで割る):', np.std(x))
print('標準偏差(N-1で割る):', np.std(x, ddof=1)) # ddof = Delta Degrees of Freedom: 自由度
print('尖度:', stats.kurtosis(x))
print('歪度:', stats.skew(x))

出力結果

平均値: 5.0
中央値: 5.0
四分位点(25%): 4.0
四分位点(50%): 5.0
四分位点(75%): 6.0
最頻値: ModeResult(mode=array([3]), count=array([1]))
合計: 25
偏差: [0.0, -2.0, -1.0, 2.0, 1.0]
偏差の平方: [0.0, 4.0, 1.0, 4.0, 1.0]
偏差平方和: 10.0
標本分散: 2.0
不偏分散: 2.5
標準偏差(Nで割る): 1.41421356237
標準偏差(N-1で割る): 1.58113883008
尖度: -1.3
歪度: 0.0
■自作関数での計算

平均値、分散、標準偏差などの統計値はnumpyモジュールの関数で簡単に出力することができますが、計算過程の理解を深めるために、関数を自作して実行し、上記と同じ結果になることを確認します。

def mean(a):
    return sum(a) / len(a)

def var(a, ddof=0):
    m = mean(a)
    v = 0.0
    for x in a:
        v += (x - m) ** 2
    return v / (len(a) - ddof)

def std(a, ddof=0):
    m = mean(a)
    v = 0.0
    for x in a:
        v += (x - m) ** 2
    return math.sqrt(v / (len(a) - ddof))

print('標本分散:', var(x))
print('不偏分散:', var(x, ddof=1))
print('標準偏差(Nで割る):', std(x))
print('標準偏差(N-1で割る):', std(x, ddof=1))

出力結果

標本分散: 2.0
不偏分散: 2.5
標準偏差(Nで割る): 1.4142135623730951
標準偏差(N-1で割る): 1.5811388300841898

今後のために各種統計量を出力する関数を作成しておきます。

def print_statistics_values(x):
    print('平均値:', np.mean(x))
    print('中央値:', np.median(x))
    print('四分位点(25%):', np.percentile(x, 25))
    print('四分位点(50%):', np.percentile(x, 50))
    print('四分位点(75%):', np.percentile(x, 75))
    print('最頻値:', stats.mode(x))
    print('合計:', sum(x))
    ave =  np.mean(x)
    print('偏差:', [d - ave for d in x])
    sq = [(d - ave) ** 2 for d in x]
    print('偏差の平方:', sq)
    print('偏差平方和:', sum(sq))
    print('標本分散:', np.var(x))
    print('不偏分散:', np.var(x, ddof=1))
    print('標準偏差(Nで割る):', np.std(x))
    print('標準偏差(N-1で割る):', np.std(x, ddof=1))
    print('尖度:', stats.kurtosis(x))
    print('歪度:', stats.skew(x))
    print('\n')

基準値

【例題】あるクラスの国語と数学の点数

No. 国語 数学
1 90 93
2 57 90
3 56 80
4 54 63
5 53 55
6 52 45
7 50 40
8 45 27
9 40 20
10 33 17
平均値 53.0 53.0
標準偏差 15.0 28.1

基準化した値=基準値

x = [90, 57, 56, 54, 53, 52, 50, 45, 40, 33]
y = [93, 90, 80, 63, 55, 45, 40, 27, 20, 17]

ave_x =  np.mean(x)
ave_y =  np.mean(y)
std_x = np.std(x, ddof=1)
std_y = np.std(y, ddof=1)
print('国語の平均値:', ave_x)
print('数学の平均値:', ave_y)
print('国語の標準偏差:', std_x)
print('数学の標準偏差:', std_y)

出力結果

国語の平均値: 53.0
数学の平均値: 53.0
国語の標準偏差: 15.0480711204
数学の標準偏差: 28.0792529182
国語の基準値: [2.458786890623796, 0.26581479898635629, 0.19936109923976725, 0.066453699746589073, 0.0, -0.066453699746589073, -0.19936109923976725, -0.53162959797271259, -0.86389809670565798, -1.3290739949317816]
数学の基準値: [1.4245393250497873, 1.3176988756710533, 0.96156404440860632, 0.35613483126244683, 0.071226966252489365, -0.28490786500995746, -0.46297528064118082, -0.92595056128236164, -1.1752449431660745, -1.2820853925448086]

基準値の平均は0、標準偏差は1になります。

print('国語の基準値の平均:', np.mean(norm_x))
print('数学の基準値の平均:', np.mean(norm_y))
print('国語の基準値の標準偏差:', np.std(norm_x, ddof=1))
print('数学の基準値の標準偏差:', np.std(norm_y, ddof=1))

出力結果
※pythonのfloat形式の計算の過程上、平均値は完全には0になっていません。

国語の基準値の平均: 4.4408920985e-17
数学の基準値の平均: -2.22044604925e-17
国語の基準値の標準偏差: 1.0
数学の基準値の標準偏差: 1.0

偏差値に変換してみます。(見やすくするために小数点以下を切り捨てます。)

devval_x = [int(50 + (d * 10)) for d in norm_x]
devval_y = [int(50 + (d * 10)) for d in norm_y]
print('国語の偏差値:', devval_x)
print('数学の偏差値:', devval_y)

出力結果

国語の偏差値: [74, 52, 51, 50, 50, 49, 48, 44, 41, 36]
数学の偏差値: [64, 63, 59, 53, 50, 47, 45, 40, 38, 37]

二変量解析の演習

■標準関数、標準モジュールでの計算

量的変数×量的変数

【例題】身長と体重の相関

No. 身長x 体重y xの
偏差
yの
偏差
xの
偏差2乗
yの
偏差2乗
偏差の積
A 146 45 -4 -5 16 25 20
B 145 46 -5 -4 25 16 20
C 147 47 -3 -3 9 9 9
D 149 49 -1 -1 1 1 1
E 151 48 1 -2 1 4 -2
F 149 51 -1 1 1 1 -1
G 151 52 1 2 1 4 2
H 154 53 4 3 16 9 12
I 153 54 3 4 9 16 12
J 155 55 5 5 25 25 25
平均値 150 50
合計 104 110 98

相関係数の計算

基本統計量を確認してみます。

x = [146, 145, 147, 149, 151, 149, 151, 154, 153, 155]
y = [45, 46, 47, 49, 48, 51, 52, 53, 54, 55]
print_statistics_values(x)
print_statistics_values(y)

出力結果

平均値: 150.0
中央値: 150.0
四分位点(25%): 147.5
四分位点(50%): 150.0
四分位点(75%): 152.5
最頻値: ModeResult(mode=array([149]), count=array([2]))
合計: 1500
偏差: [-4.0, -5.0, -3.0, -1.0, 1.0, -1.0, 1.0, 4.0, 3.0, 5.0]
偏差の平方: [16.0, 25.0, 9.0, 1.0, 1.0, 1.0, 1.0, 16.0, 9.0, 25.0]
偏差平方和: 104.0
標本分散: 10.4
不偏分散: 11.5555555556
標準偏差(Nで割る): 3.22490309932
標準偏差(N-1で割る): 3.3993463424
尖度: -1.2174556213017753
歪度: 0.0

平均値: 50.0
中央値: 50.0
四分位点(25%): 47.25
四分位点(50%): 50.0
四分位点(75%): 52.75
最頻値: ModeResult(mode=array([45]), count=array([1]))
合計: 500
偏差: [-5.0, -4.0, -3.0, -1.0, -2.0, 1.0, 2.0, 3.0, 4.0, 5.0]
偏差の平方: [25.0, 16.0, 9.0, 1.0, 4.0, 1.0, 4.0, 9.0, 16.0, 25.0]
偏差平方和: 110.0
標本分散: 11.0
不偏分散: 12.2222222222
標準偏差(Nで割る): 3.31662479036
標準偏差(N-1で割る): 3.4960294939
尖度: -1.3818181818181818
歪度: 0.0

相関係数を計算します。

ave_x =  np.mean(x)
ave_y =  np.mean(y)
sq_x = [(d - ave_x) ** 2 for d in x]
sq_y = [(d - ave_y) ** 2 for d in y]
product_xy = [(dx - ave_x) * (dy - ave_y) for dx, dy in zip(x, y)]
print('xの偏差2乗:', sq_x)
print('yの偏差2乗:', sq_y)
print('偏差の積:', product_xy)

s_xx = sum(sq_x)
s_yy = sum(sq_y)
s_xy = sum(product_xy)
print('Sxx:', s_xx)
print('Syy:', s_yy)
print('Sxy:', s_xy)
r = s_xy / math.sqrt(s_xx * s_yy)
print('相関係数', r)

出力結果

xの偏差2乗: [16.0, 25.0, 9.0, 1.0, 1.0, 1.0, 1.0, 16.0, 9.0, 25.0]
yの偏差2乗: [25.0, 16.0, 9.0, 1.0, 4.0, 1.0, 4.0, 9.0, 16.0, 25.0]
偏差の積: [20.0, 20.0, 9.0, 1.0, -2.0, -1.0, 2.0, 12.0, 12.0, 25.0]
Sxx: 104.0
Syy: 110.0
Sxy: 98.0
相関係数 0.916248050208

単回帰分析

単回帰式の傾きaと切片bを求め、式を導出します。

a = s_xy / s_xx
b = ave_y - a * ave_x
print('傾きa:', a)
print('切片b:', b)
print('単回帰式:y=' + str(a) + 'x' +str(b))

出力結果

傾きa: 0.942307692308
切片b: -91.3461538462
単回帰式:y=0.942307692308x-91.3461538462

決定係数

回帰式による予測値と実測値との差を計算します。

No. 身長x 体重y
実測値
体重y
予測値
残差 残差の
平方
A 146 45 46.2 1.2 1.44
B 145 46 45.3 -0.7 0.49
C 147 47 47.2 0.2 0.04
D 149 49 49.1 0.1 0.01
E 151 48 50.9 2.9 8.41
F 149 51 49.1 -1.9 3.61
G 151 52 50.9 -1.1 1.21
H 154 53 53.8 0.8 0.64
I 153 54 52.8 -1.2 1.44
J 155 55 54.7 -0.3 0.09
合計 0 17.38
exp = [a * d + b for d in x]
diff = [de - dy for de, dy in zip(exp, y)]
diff_sq = [d ** 2 for d in diff]
print('予測値:', exp)
print('残差:', diff)
print('残差の平方:', diff_sq)
s_e = sum(diff_sq)
r2 = 1 - s_e / s_yy
print('残差平方和Se:', s_e)
print('目的変数yの偏差平方和Syy:', s_yy)
print('決定係数r2:', r2)

出力結果

予測値: [46.230769230769226, 45.288461538461547, 47.173076923076934, 49.057692307692321, 50.942307692307708, 49.057692307692321, 50.942307692307708, 53.769230769230774, 52.826923076923094, 54.711538461538453]
残差: [1.2307692307692264, -0.71153846153845279, 0.17307692307693401, 0.05769230769232081, 2.9423076923077076, -1.9423076923076792, -1.0576923076922924, 0.7692307692307736, -1.1730769230769056, -0.28846153846154721]
残差の平方: [1.5147928994082733, 0.50628698224850821, 0.029955621301778933, 0.003328402366865419, 8.6571745562131071, 3.7725591715975821, 1.1187130177514468, 0.59171597633136763, 1.3761094674555803, 0.083210059171602679]
残差平方和Se: 17.6538461538
目的変数yの偏差平方和Syy: 110.0
決定係数r2: 0.83951048951

決定係数r2の平方根は相関係数rに等しいことを確認します。

print('決定係数r2の平方根:', math.sqrt(r2), '=相関係数r')

出力結果

決定係数r2の平方根: 0.9162480502082883 =相関係数r

質的変数×量的変数

【例題】年齢と好きな商品の相関

No. 年齢(歳) 好きな商品
1 29 A
2 32 A
3 35 A
4 36 A
5 38 B
6 40 B
7 41 B
8 43 B
9 48 B
10 20 C
11 22 C
12 24 C
13 29 C
14 35 C
15 38 C

質的変数×質的変数

【例題】中学生の学年とゲーム機購入予定の相関

No. 学年 好きな商品
1 1年 ある
2 1年 ある
3 1年 ない
4 2年 ある
5 2年 ある
6 2年 ある
7 2年 ない
8 3年 ある
9 3年 ない
10 3年 ない

確率分布

matplotlibをグラフ描画に利用します。

import matplotlib.pyplot as plt

二項分布

成功確率pの試行をn回行った場合の二項分布から生成される乱数を生成します。

x = np.random.binomial(3, 0.5, 100) # n=3, p=0.5の二項分布、データ数100個
plt.hist(x, bins=10) # 描く棒の数10個
plt.show()

x = np.random.binomial(10, 0.5, 300) # n=10, p=0.5の二項分布、データ数300個
plt.hist(x, bins=20)
plt.show()

x = np.random.binomial(100, 0.5, 3000) # n=100, p=0.5の二項分布、データ数3000個
plt.hist(x, bins=200)
plt.show()

試行回数nを増やせば増やすほど、正規分布の形に近づいていくことがわかります。
成功確率pの値を変えると、グラフの頂点部分が移動することがわかります。

ポアソン分布

単位時間当たり平均λ回発生する事象について、発生回数の乱数を生成します。

x = np.random.poisson(10, 100) # λ=10のポアソン分布、データ数100個
plt.hist(x, bins=20)
plt.show()

たとえば1時間当たり10回クリックされるWEB広告において、それぞれクリックされる回数毎の確率がわかります。

正規分布

標準正規分布(平均0、標準偏差1)の乱数を生成します。

R = np.random.randn(10000) # 標準正規分布で乱数を1万個生成
plt.hist(R, bins=100)
plt.show()

 

 

カテゴリー