왜도와 첨도¶

import numpy as np
from scipy import stats

np.random.seed(0)

# 정규분포는 좌우대칭이므로 왜도가 0 근처로 나온다. 표준편차를 3 으로
# 키워도 마찬가지다 — 왜도는 척도에 휘둘리지 않는 값이다.
data = np.random.normal(0, 3, 1000)
# data = np.random.exponential(1, 1000)

skewness_value = stats.skew(data)
print(f"Skewness: {skewness_value:.4f}")

Skewness: 0.0339

import numpy as np
from scipy import stats

np.random.seed(0)

# 같은 씨앗에서 척도만 바꾸어 뽑으면 자료는 정확히 3배가 된다.
a = np.random.normal(0, 3, 1000)
np.random.seed(0)
b = np.random.normal(0, 1, 1000)
print(f"a 가 b 의 정확히 3배인가: {np.array_equal(a, 3 * b)}")
print(f"skew(a) = {stats.skew(a):.10f}")
print(f"skew(b) = {stats.skew(b):.10f}")

# 표준오차: 보정판 G1 에 대한 정확식과 어림값.
n = 1000
se_G1 = np.sqrt(6 * n * (n - 1) / ((n - 2) * (n + 1) * (n + 3)))
se_g1 = se_G1 * (n - 2) / np.sqrt(n * (n - 1))
print(f"\nSE(G1) 정확식 = {se_G1:.5f},  sqrt(6/n) = {np.sqrt(6 / n):.5f}")
print(f"SE(g1) = SE(G1) * (n-2)/sqrt(n(n-1)) = {se_g1:.5f}")

# 모의실험으로 두 판본의 표준편차를 직접 잰다.
rng = np.random.default_rng(11)
R = 50_000
X = rng.standard_normal((R, n))
sd_g1 = stats.skew(X, axis=1).std(ddof=1)
sd_G1 = stats.skew(X, axis=1, bias=False).std(ddof=1)
print(f"모의 sd(g1) = {sd_g1:.5f},  모의 sd(G1) = {sd_G1:.5f}  (R = {R})")
print(f"관측된 g1 = {stats.skew(b):.4f} 은 SE 의 {stats.skew(b) / sd_g1:.2f} 배")

a 가 b 의 정확히 3배인가: True
skew(a) = 0.0338589532
skew(b) = 0.0338589532

SE(G1) 정확식 = 0.07734,  sqrt(6/n) = 0.07746
SE(g1) = SE(G1) * (n-2)/sqrt(n(n-1)) = 0.07723
모의 sd(g1) = 0.07741,  모의 sd(G1) = 0.07753  (R = 50000)
관측된 g1 = 0.0339 은 SE 의 0.44 배

import numpy as np
from scipy import stats

np.random.seed(0)

# 이번에는 지수분포다. 오른쪽으로 길게 늘어졌으므로 왜도가 양으로 크게 나온다.
# data = np.random.normal(0, 1, 1000)
data = np.random.exponential(1, 1000)

skewness_value = stats.skew(data)
print(f"Skewness: {skewness_value:.4f}")

# scipy 의 kurtosis 는 기본이 초과첨도다(정규가 0). 다른 책에서 3 을 더한
# 값을 쓰는 일이 많으니 견줄 때 주의해야 한다.
kurtosis_value = stats.kurtosis(data)
print(f"Kurtosis: {kurtosis_value:.4}")

Skewness: 2.0526
Kurtosis: 6.476

import numpy as np
from scipy import stats

# 지수분포의 이론 적률: E[X^k] = k! (비율 1)
mu = 1.0
m2 = 2 - 1**2
m3 = 6 - 3 * 2 * 1 + 3 * 1 * 1**2 - 1**3
m4 = 24 - 4 * 6 * 1 + 6 * 2 * 1**2 - 4 * 1 * 1**3 + 1**4
print(f"중심적률  mu2 = {m2:.0f}, mu3 = {m3:.0f}, mu4 = {m4:.0f}")
print(f"이론 왜도   = mu3 / mu2^1.5 = {m3 / m2**1.5:.0f}")
print(f"이론 초과첨도 = mu4 / mu2^2 - 3 = {m4 / m2**2 - 3:.0f}")

# 표본값이 이론값에서 얼마나 떨어질 수 있는가 — 표집분포를 직접 잰다.
rng = np.random.default_rng(3)
R, n = 20_000, 1000
Y = rng.exponential(1.0, size=(R, n))
g1 = stats.skew(Y, axis=1)
g2 = stats.kurtosis(Y, axis=1)
print(f"\n지수 자료 n = {n} 의 표집분포 (R = {R})")
print(f"  g1: 평균 {g1.mean():.3f}  표준편차 {g1.std(ddof=1):.3f}")
print(f"  g2: 평균 {g2.mean():.3f}  표준편차 {g2.std(ddof=1):.3f}")
print(f"  관측 g1 = 2.0526 은 평균에서 {(2.0526 - g1.mean()) / g1.std(ddof=1):+.2f} 표준편차")
print(f"  관측 g2 = 6.476  은 평균에서 {(6.476 - g2.mean()) / g2.std(ddof=1):+.2f} 표준편차")

중심적률  mu2 = 1, mu3 = 2, mu4 = 9
이론 왜도   = mu3 / mu2^1.5 = 2
이론 초과첨도 = mu4 / mu2^2 - 3 = 6

지수 자료 n = 1000 의 표집분포 (R = 20000)
  g1: 평균 1.974  표준편차 0.251
  g2: 평균 5.741  표준편차 2.470
  관측 g1 = 2.0526 은 평균에서 +0.31 표준편차
  관측 g2 = 6.476  은 평균에서 +0.30 표준편차

import numpy as np
from scipy import stats
d = np.array([1, 2, 2, 3, 3, 3, 4, 4, 5, 100.0])
print(stats.skew(d))       # 2.6596
print(stats.kurtosis(d))   # 5.0891

2.659594282898852
5.089118705274686

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
n = 300_000

print(f"{'분포':>8}{'왜도':>10}{'초과첨도':>12}{'왜도^2':>10}{'초과첨도+2':>13}")
for label, x in [("정규", rng.normal(0, 1, n)), ("지수", rng.exponential(1, n)),
                 ("로그정규", rng.lognormal(0, 1, n)),
                 ("균등", rng.uniform(-1, 1, n))]:
    s, k = stats.skew(x), stats.kurtosis(x)
    print(f"{label:>8}{s:>10.3f}{k:>12.3f}{s * s:>10.3f}{k + 2:>13.3f}")

      분포        왜도        초과첨도      왜도^2       초과첨도+2
      정규    -0.005       0.006     0.000        2.006
      지수     2.002       5.966     4.007        7.966
    로그정규     5.978      91.554    35.735       93.554
      균등    -0.001      -1.200     0.000        0.800

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
print(f"{'n':>6}{'SE(왜도)':>22}{'SE(첨도)':>24}")
print(f"{'':>6}{'기본값':>11}{'bias=False':>11}{'이론':>10}{'bias=False':>13}{'이론':>10}")
for n in (20, 50, 200, 1000):
    s_def = np.std([stats.skew(rng.normal(0, 1, n)) for _ in range(20_000)])
    s_unb = np.std([stats.skew(rng.normal(0, 1, n), bias=False) for _ in range(20_000)])
    k_unb = np.std([stats.kurtosis(rng.normal(0, 1, n), bias=False) for _ in range(20_000)])
    se_s = np.sqrt(6 * n * (n - 1) / ((n - 2) * (n + 1) * (n + 3)))
    se_k = np.sqrt(24 * n * (n - 1) ** 2 / ((n - 3) * (n - 2) * (n + 3) * (n + 5)))
    print(f"{n:>6}{s_def:>11.4f}{s_unb:>11.4f}{se_s:>10.4f}{k_unb:>13.4f}{se_k:>10.4f}")

     n                SE(왜도)                  SE(첨도)
              기본값 bias=False        이론   bias=False        이론
    20     0.4712     0.5143    0.5121       0.9969    0.9924
    50     0.3246     0.3375    0.3366       0.6609    0.6619
   200     0.1707     0.1719    0.1719       0.3461    0.3422
  1000     0.0773     0.0767    0.0773       0.1541    0.1545