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 の推定値も求めることができる。¹

観察値 $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$²

$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 の推定値も求めることができる。¹

観察値 $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$）。すなわち，正規分布からの標本ではないと結論する。

D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189. ↩ ↩²
fisher=True のとき ↩