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番ポートを利用しますので、ネットワーク環境に応じて該当ポートを開放してください。

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

一変量解析

二変量解析

確率分布

一変量解析の演習

基本統計量

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

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

各データを変数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')

基準値

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

基準化した値＝基準値

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]

二変量解析の演習

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

量的変数×量的変数

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

相関係数の計算

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

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

決定係数

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

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

質的変数×量的変数

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

質的変数×質的変数

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

確率分布

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()