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