scipy.stats: 尖度と歪度

# 1.　尖度 `kurtosis`

`kurtosis(a, axis=0, fisher=True, bias=True, nan_policy='propagate')`


```python
import numpy as np
from scipy.stats import kurtosis
x = [1.2, 3.1, 4.6, 2.2, 5.7]
```

## 1.1.　Fisher の定義(正規分布のとき 0 になる)


```python
kurtosis(x) # g2:後述
```

    -1.3684571122281206

関数を使わずに計算すると以下のようになる。

```python
from scipy.stats import zscore
np.mean(zscore(x)**4) - 3
```

    -1.3684571122281206

## 1.2.　Pearson の定義(正規分布のとき 3 になる) `fisher=False`

```python
kurtosis(x, fisher=False)
```

    1.6315428877718794

関数を使わずに計算すると以下のようになる。

```python
from scipy.stats import zscore
np.mean(zscore(x)**4)
```

    1.6315428877718794

## 1.3.　不偏推定量 `bias=False`

### 1.3.1.　`fisher=True`

```python
kurtosis(x, bias=False) # G2:後述
```

    -1.4738284489124838

関数を使わずに計算すると以下のようになる。


```python
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`

```python
kurtosis(x, fisher=False, bias=False)
```

    1.5261715510875162

関数を使わずに計算すると以下のようになる。

```python
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$ になる）。

[^1]: D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189.

```python
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$

[^2]: `fisher=True` のとき


```python
from scipy.stats import kurtosistest, norm, t
np.random.seed(123)
```

## 2.1.　例1.

正規乱数について，両側検定を行う

```python
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

```python
kurtosis(y)
```

    -0.8826271976238784

```python
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$ 分布すなわちコーシー分布からの乱数について，両側検定を行う。

```python
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

```python
kurtosis(z)
```

    18.137971157631473

```python
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')`

```python
import numpy as np
from scipy.stats import skew
x = [1.2, 3.1, 4.6, 2.2, 5.7]
```

## 3.1.　不偏ではない推定量

何のオプションも指定しないときは，古典的な定義による歪度を求める。

```python
skew(x) # g1 後述
```

    0.14460191499270095

関数を使わずに計算すると以下のようになる。

```python
from scipy.stats import zscore
np.mean(zscore(x)**3)
```

    0.14460191499270109

## 3.2.　不偏推定量

不偏推定量を求めるためには `bias=False` を指定する。Excel, SAS, SPSS が採用している方法である。

```python
skew(x, bias=False) # G1 後述
```

    0.2155598077335502

関数を使わずに計算すると以下のようになる。

```python
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$ になる）。

[^3]: D. N. Joanes and C. A. Gill (1998), Comparing measures of sample skewness and kurtosis. The Statistician, 47, 183–189.

```python
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')`

```python
from scipy.stats import skewtest, norm, t
np.random.seed(123)
```

## 4.1.　例1.

正規乱数について，両側検定を行う。

```python
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

```python
skew(y)
```

    -0.009096852799368121

```python
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$ 分布すなわちコーシー分布からの乱数について，両側検定を行う。

```python
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

```python
skew(z)
```

    4.365811733349633

```python
z_value, p_value = skewtest(z)
(z_value, p_value)
```

    (5.825572775476851, 5.691704942463112e-09)

$p \lt 0.05$ なので，帰無仮説は棄却される（$p \lt 0.001$）。すなわち，正規分布からの標本ではないと結論する。