1. 尖度 kurtosis
kurtosis(a, axis=0, fisher=True, bias=True, nan_policy='propagate')
import numpy as np
from scipy.stats import kurtosis
x = [1.2, 3.1, 4.6, 2.2, 5.7]
1.1. Fisher の定義(正規分布のとき 0 になる)
kurtosis(x) # g2:後述
-1.3684571122281206
関数を使わずに計算すると以下のようになる。
from scipy.stats import zscore
np.mean(zscore(x)**4) - 3
-1.3684571122281206
1.2. Pearson の定義(正規分布のとき 3 になる) fisher=False
kurtosis(x, fisher=False)
1.6315428877718794
関数を使わずに計算すると以下のようになる。
from scipy.stats import zscore
np.mean(zscore(x)**4)
1.6315428877718794
1.3. 不偏推定量 bias=False
1.3.1. fisher=True
kurtosis(x, bias=False) # G2:後述
-1.4738284489124838
関数を使わずに計算すると以下のようになる。
n = len(x)
n*(n+1)*sum(zscore(x, ddof=1)**4)/(n-1)/(n-2)/(n-3) - 3*(n-1)**2/(n-2)/(n-3)
-1.4738284489124807
1.3.2. fisher=False
kurtosis(x, fisher=False, bias=False)
1.5261715510875162
関数を使わずに計算すると以下のようになる。
n*(n+1)*sum(zscore(x,ddof=1)**4)/(n-1)/(n-2)/(n-3) - 3*(n-1)**2/(n-2)/(n-3) + 3
1.5261715510875193
1.4. 第 3 の定義 b2
R の e1071::kurtosis は以上の 2 種類の尖度の他に，第 3 の推定値も求めることができる。1
観察値 $x_i$，サンプルサイズ $n$，平均値 $\mu$，標準偏差 $s$ とおき，$r$ 次のモーメントを $m_r = \sum_i(x_i - \mu)^r\ /\ n$ としたとき，以下の 3 通りの定義がある。
- $g_2 = m_4\ /\ m_2^2 - 3$
- $G_2 = {(n+1)g_2 + 6} (n-1) \ /\ {(n-2)(n-3)}$
- $b_2 = m_4\ /\ s^4 - 3 = (g_2 + 3)(1-1/n)^2 - 3$
$g_2$ は kurtosis(x) で求められる，多くの古い教科書に載っているものである。
$G_2$ は kurtosis(x, bias=False) で求められる，SAS や SPSS や Excel でも使われているものである。
$b_2$ は MINITAB， BMDP で使われているもので，e1071::kurtosis でも求めることができる。$g_2$ から求めることもできるが，np.mean(zscore(x, ddof=1)**4) - 3 でも計算できる（ddof=0 ならば $g_2$ になる）。
g2 = kurtosis(x)
G2 = ((n+1)*g2 + 6) * (n-1) / (n-2) / (n-3)
b2 = (g2 + 3)*(1-1/n)**2 - 3
print(g2, G2, b2)
print(np.mean(zscore(x, ddof=1)**4) - 3)
-1.3684571122281206 -1.4738284489124827 -1.955812551825997
-1.9558125518259968
R の moments::kurtosis（x） は kurtosis(x, fisher=False) すなわち $m_4\ /\ m_2^2 = g_2+3$ を求める。
2. 尖度の検定 kurtosistest
kurtosistest(a, axis=0, nan_policy='propagate', alternative='two-sided')
$H_0:\text{尖度} = 0$2
$H_1:\text{尖度} \neq 0$
from scipy.stats import kurtosistest, norm, t
np.random.seed(123)
2.1. 例1.
正規乱数について，両側検定を行う
y = norm.rvs(loc=50, scale=10, size=25).round(1)
print(*y)
39.1 60.0 52.8 34.9 44.2 66.5 25.7 45.7 62.7 41.3 43.2 49.1 64.9 43.6 45.6 45.7 72.1 71.9 60.0 53.9 57.4 64.9 40.6 61.8 37.5
kurtosis(y)
-0.8826271976238784
z_value, p_value = kurtosistest(y)
(z_value, p_value)
(-1.032048194835839, 0.3020495406844782)
$p \gt 0.05$ なので，帰無仮説は棄却されない（$p = 0.302$）。すなわち，正規分布からの標本でないとはいえない。
2.2. 例2.
自由度 1 の $t$ 分布すなわちコーシー分布からの乱数について，両側検定を行う。
z = t.rvs(df=1, loc=50, scale=10, size=25).round(1) # コーシー分布
print(*z)
34.6 47.7 35.7 54.1 46.1 44.3 57.8 70.2 88.0 50.3 61.7 21.9 108.7 671.3 87.3 64.7 44.6 40.6 48.8 64.7 17.8 39.9 -31.0 57.5 42.0
kurtosis(z)
18.137971157631473
z_value, p_value = kurtosistest(z)
(z_value, p_value)
(4.954534755281744, 7.250355695269661e-07)
$p \lt 0.05$ なので，帰無仮説は棄却される（$p \lt 0.001$）。すなわち，正規分布からの標本ではないと結論する。
3. 歪度 skew
- データの分布が正規分布にしたがうとき，歪度はほぼ 0 である。
- 0 より大きいときは分布の重心が右にある（左の裾が長い）。
- 0 より小さいときは分布の重心が左にある（右の裾が長い）。
skew(a, axis=0, bias=True, nan_policy='propagate')
import numpy as np
from scipy.stats import skew
x = [1.2, 3.1, 4.6, 2.2, 5.7]
3.1. 不偏ではない推定量
何のオプションも指定しないときは，古典的な定義による歪度を求める。
skew(x) # g1 後述
0.14460191499270095
関数を使わずに計算すると以下のようになる。
from scipy.stats import zscore
np.mean(zscore(x)**3)
0.14460191499270109
3.2. 不偏推定量
不偏推定量を求めるためには bias=False を指定する。Excel, SAS, SPSS が採用している方法である。
skew(x, bias=False) # G1 後述
0.2155598077335502
関数を使わずに計算すると以下のようになる。
n = len(x)
n*sum(zscore(x, ddof=1)**3)/(n-1)/(n-2)
0.2155598077335504
3.3. 第 3 の推定量
R の e1071::skewness は以上の 2 種類の歪度の他に，第 3 の推定値も求めることができる。1
観察値 $x_i$，サンプルサイズ $n$，平均値 $\mu$，標準偏差 $s$ とおき，$r$ 次のモーメントを $m_r = \sum_i(x_i - \mu)^r\ /\ n$ としたとき，以下の 3 通りの定義がある。
- $g_1 = m_3\ /\ m_2^{3/2}$
- $G_1 = g_1\sqrt{n\ (n-1)}\ /\ (n-2)$
- $b_1 = m_3\ /\ s^3 = g_1\ {(n-1)\ /\ n}^{3/2}$
$g_1$ は kurtosis(x) で求められる，多くの古い教科書に載っているものである。
$G_1$ は kurtosis(x, bias=False) で求められる，SAS や SPSS や Excel でも使われているものである。
$b_1$ は MINITAB， BMDP で使われているもので，e1071::skewness でも求めることができる。$g_1$ から求めることもできるが，np.mean(zscore(x, ddof=1)**3) でも計算できる（ddof=0 ならば $g_1$ になる）。
g1 = skew(x)
G1 = g1 * np.sqrt(n*(n-1)) / (n-2)
b1 = g1*((n-1)/n)**(3/2)
print(g1, G1, b1)
print(np.mean(zscore(x, ddof=1)**3))
0.14460191499270095 0.2155598077335502 0.10346870771210412
0.10346870771210419
R の moments::skewness(x） は skew(x) すなわち $m_3\ /\ m_2^{3/2} = g_1$ を求める。
4. 歪度の検定 skewtest
$H_0:\text{歪度} = 0$
$H_1:\text{歪度} \neq 0$
skewtest(a, axis=0, nan_policy='propagate', alternative='two-sided')
from scipy.stats import skewtest, norm, t
np.random.seed(123)
4.1. 例1.
正規乱数について，両側検定を行う。
y = norm.rvs(loc=50, scale=10, size=25).round(1)
print(*y)
39.1 60.0 52.8 34.9 44.2 66.5 25.7 45.7 62.7 41.3 43.2 49.1 64.9 43.6 45.6 45.7 72.1 71.9 60.0 53.9 57.4 64.9 40.6 61.8 37.5
skew(y)
-0.009096852799368121
z_value, p_value = skewtest(y)
(z_value, p_value)
(-0.022210278084285384, 0.982280218888247)
$p \gt 0.05$ なので，帰無仮説は棄却されない（$p = 0.982$）。すなわち，正規分布からの標本でないとはいえない。
4.2. 例2.
自由度 1 の $t$ 分布すなわちコーシー分布からの乱数について，両側検定を行う。
z = t.rvs(df=1, loc=50, scale=10, size=25).round(1) # コーシー分布
print(*z)
34.6 47.7 35.7 54.1 46.1 44.3 57.8 70.2 88.0 50.3 61.7 21.9 108.7 671.3 87.3 64.7 44.6 40.6 48.8 64.7 17.8 39.9 -31.0 57.5 42.0
skew(z)
4.365811733349633
z_value, p_value = skewtest(z)
(z_value, p_value)
(5.825572775476851, 5.691704942463112e-09)
$p \lt 0.05$ なので，帰無仮説は棄却される（$p \lt 0.001$）。すなわち，正規分布からの標本ではないと結論する。