왜도와 첨도¶

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

def generate_and_plot_mixed_distribution(seed: int = 0):
    """중심은 같고 퍼짐만 다른 정규분포 셋을 섞는다.

    loc(중심)은 모두 0으로 두고 scale(표준편차)만 1, 2, 4로 키운다.
    좌우가 똑같이 늘어나므로 **대칭이면서 꼬리만 두꺼운** 분포가 된다.
    즉 왜도는 0 근처, 첨도는 정규분포보다 크게 나온다.
    """
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)             # 본체 1000개, 표준편차 1
    minor_1 = stats.norm(scale=2).rvs(200)          # 조금 넓게 200개
    minor_2 = stats.norm(scale=4).rvs(100)          # 아주 넓게 100개
    combined = np.concatenate((main_data, minor_1, minor_2))

    fig, ax = plt.subplots(figsize=(12, 3))
    ax.hist(combined, bins=30)
    plt.show()

    # 눈으로 본 것을 숫자로 확인한다.
    # fisher=True(기본)이면 정규분포의 첨도가 0이 되도록 3을 뺀 초과첨도다.
    print(f"왜도 {stats.skew(combined):+.3f}  (0에 가까움 = 대칭)")
    print(f"첨도 {stats.kurtosis(combined):+.3f}  (0보다 큼 = 정규분포보다 꼬리가 두껍다)")

if __name__ == "__main__":
    generate_and_plot_mixed_distribution()

왜도 -0.008  (0에 가까움 = 대칭)
첨도 +7.150  (0보다 큼 = 정규분포보다 꼬리가 두껍다)

import numpy as np
from scipy import stats

w = np.array([10, 2, 1]) / 13
s = np.array([1.0, 2.0, 4.0])
mu2 = np.sum(w * s ** 2)
mu4 = np.sum(w * 3 * s ** 4)        # 정규분포의 4차 적률은 3 sigma^4
print(f"mu2 = {mu2:.6f}  (= 34/13)")
print(f"mu4 = {mu4:.6f}  (= 894/13)")
print(f"beta2 = {mu4 / mu2 ** 2:.6f},  gamma2 = {mu4 / mu2 ** 2 - 3:.6f}  (= 4077/578)")
print(f"sigma = {np.sqrt(mu2):.6f}")

# (2) 의 공식 gamma2 = 3 Var(V)/(E V)^2 로도 같은 값이 나오는가. V = sigma_i^2
V = s ** 2
EV = np.sum(w * V)
VarV = np.sum(w * (V - EV) ** 2)
print(f"3 * Var(V) / (E V)^2 = {3 * VarV / EV ** 2:.6f}")

# 표본을 키우면 모집단 값에 닿는가
rng = np.random.default_rng(0)
print(f"\n{'n':>12}{'g1':>10}{'g2':>10}")
for N in (13_000, 130_000, 1_300_000, 13_000_000):
    k = np.array([10, 2, 1]) * (N // 13)
    x = np.concatenate([rng.normal(0, s[i], k[i]) for i in range(3)])
    print(f"{len(x):>12,}{stats.skew(x):>+10.4f}{stats.kurtosis(x):>+10.4f}")
print(f"{'모집단':>12}{0.0:>+10.4f}{mu4 / mu2 ** 2 - 3:>+10.4f}")

mu2 = 2.615385  (= 34/13)
mu4 = 68.769231  (= 894/13)
beta2 = 10.053633,  gamma2 = 7.053633  (= 4077/578)
sigma = 1.617215
3 * Var(V) / (E V)^2 = 7.053633

           n        g1        g2
      13,000   -0.1426   +6.1992
     130,000   +0.0102   +7.0013
   1,300,000   -0.0018   +6.9807
  13,000,000   +0.0035   +7.0533
         모집단   +0.0000   +7.0536

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def generate_and_plot_right_skewed_distribution(seed: int = 0):
    """중심을 오른쪽으로만 옮긴 덩어리를 덧붙인다.

    앞 보기와 달리 scale이 아니라 loc를 바꾼다.
    0, +2, +4 로 오른쪽에만 덩어리를 놓으므로 오른쪽 꼬리가 길어진다.
    """
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)             # 본체는 0 중심
    right_1 = stats.norm(loc=2).rvs(200)            # 오른쪽 어깨
    right_2 = stats.norm(loc=4).rvs(100)            # 오른쪽 꼬리
    combined = np.concatenate((main_data, right_1, right_2))

    fig, ax = plt.subplots(figsize=(12, 3))
    ax.hist(combined, bins=30)
    plt.show()

    # 오른쪽 치우침의 두 가지 신호를 확인한다
    print(f"왜도 {stats.skew(combined):+.3f}  (양수 = 오른쪽 치우침)")
    print(f"평균 {combined.mean():+.3f} > 중앙값 {np.median(combined):+.3f}")

if __name__ == "__main__":
    generate_and_plot_right_skewed_distribution()

왜도 +0.848  (양수 = 오른쪽 치우침)
평균 +0.595 > 중앙값 +0.314

import numpy as np
from scipy import stats
from scipy.optimize import brentq

w = np.array([10, 2, 1]) / 13
locs = np.array([0.0, 2.0, 4.0])

mu = np.sum(w * locs)
d = locs - mu                       # 각 성분이 혼합평균에서 떨어진 거리
m2 = np.sum(w * (d ** 2 + 1))
m3 = np.sum(w * (d ** 3 + 3 * d))
m4 = np.sum(w * (d ** 4 + 6 * d ** 2 + 3))
print(f"mu  = {mu:.6f}  (= 8/13)")
print(f"d   = {np.round(d, 6)}")
print(f"mu2 = {m2:.6f}   mu3 = {m3:.6f}   mu4 = {m4:.6f}")
print(f"gamma1 = {m3 / m2 ** 1.5:.6f}   gamma2 = {m4 / m2 ** 2 - 3:.6f}")

# 중앙값은 F(x) = 0.5 를 수치적으로 푼다.
F = lambda x: float(np.sum(w * stats.norm.cdf(x, locs, 1.0)))
med = brentq(lambda x: F(x) - 0.5, -6, 12, xtol=1e-13)
print(f"\n모집단 중앙값 = {med:.6f}")
print(f"평균 - 중앙값 = {mu - med:.6f}  (> 0 이면 어림대로다)")
print(f"(평균 - 중앙값) / sigma = {(mu - med) / np.sqrt(m2):.6f}")

print(f"\n{'':>10}{'모집단':>12}{'표본 n=1300':>14}")
print(f"{'왜도':>10}{m3 / m2 ** 1.5:>+12.4f}{0.848:>+14.4f}")
print(f"{'평균':>10}{mu:>+12.4f}{0.595:>+14.4f}")
print(f"{'중앙값':>10}{med:>+12.4f}{0.314:>+14.4f}")

mu  = 0.615385  (= 8/13)
d   = [-0.615385  1.384615  3.384615]
mu2 = 2.467456   mu3 = 3.211652   mu4 = 22.575225
gamma1 = 0.828618   gamma2 = 0.707946

모집단 중앙값 = 0.358251
평균 - 중앙값 = 0.257134  (> 0 이면 어림대로다)
(평균 - 중앙값) / sigma = 0.163695

                   모집단     표본 n=1300
        왜도     +0.8286       +0.8480
        평균     +0.6154       +0.5950
       중앙값     +0.3583       +0.3140

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

def generate_and_plot_left_skewed_distribution(seed: int = 0):
    """앞 함수의 부호만 뒤집었다. loc가 -2, -4 로 왼쪽에 놓인다."""
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)
    left_1 = stats.norm(loc=-2).rvs(200)            # 왼쪽 어깨
    left_2 = stats.norm(loc=-4).rvs(100)            # 왼쪽 꼬리
    combined = np.concatenate((main_data, left_1, left_2))

    fig, ax = plt.subplots(figsize=(12, 3))
    ax.hist(combined, bins=30)
    plt.show()

    # 부호가 정확히 반대로 나온다
    print(f"왜도 {stats.skew(combined):+.3f}  (음수 = 왼쪽 치우침)")
    print(f"평균 {combined.mean():+.3f} < 중앙값 {np.median(combined):+.3f}")

if __name__ == "__main__":
    generate_and_plot_left_skewed_distribution()

왜도 -0.853  (음수 = 왼쪽 치우침)
평균 -0.636 < 중앙값 -0.396

import numpy as np
import scipy.stats as stats


def mixture(sign):
    """loc 의 부호만 바꾸어 같은 씨앗으로 표본을 만든다."""
    np.random.seed(0)
    main = stats.norm().rvs(1_000)
    lump1 = stats.norm(loc=sign * 2).rvs(200)
    lump2 = stats.norm(loc=sign * 4).rvs(100)
    return np.concatenate((main, lump1, lump2))


right, left = mixture(+1), mixture(-1)
print(f"{'':>8}{'왜도':>12}{'초과첨도':>12}{'평균':>12}{'중앙값':>12}")
for name, x in (("오른쪽", right), ("왼쪽", left)):
    print(f"{name:>8}{stats.skew(x):>+12.6f}{stats.kurtosis(x):>+12.6f}"
          f"{x.mean():>+12.6f}{np.median(x):>+12.6f}")

print(f"\n두 표본이 서로 거울상인가: "
      f"{np.allclose(np.sort(right), -np.sort(left)[::-1])}")
print(f"왜도의 합 = {stats.skew(right) + stats.skew(left):+.6f}   (거울상이면 0)")
print(f"평균의 합 = {right.mean() + left.mean():+.6f}")

# 진짜 거울상을 만들어 보면 부호가 정확히 뒤집힌다.
mirror = -right
print(f"\n-right 의 왜도 {stats.skew(mirror):+.6f}  "
      f"(= -{stats.skew(right):.6f}),  초과첨도 {stats.kurtosis(mirror):+.6f}")
print(f"  right 의 초과첨도 {stats.kurtosis(right):+.6f}  <- 네제곱이라 그대로다")

                  왜도        초과첨도          평균         중앙값
     오른쪽   +0.847545   +0.680438   +0.595182   +0.313942
      왼쪽   -0.852891   +0.924781   -0.635588   -0.396250

두 표본이 서로 거울상인가: False
왜도의 합 = -0.005345   (거울상이면 0)
평균의 합 = -0.040406

-right 의 왜도 -0.847545  (= -0.847545),  초과첨도 +0.680438
  right 의 초과첨도 +0.680438  <- 네제곱이라 그대로다

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

def generate_and_plot_histogram_and_box_plot_mixed_distribution(seed: int = 0):
    """같은 자료를 히스토그램과 상자그림으로 나란히 본다.

    자료는 첫 보기와 같은 대칭·두꺼운꼬리 혼합분포다.
    두 그림을 위아래로 붙여 x축을 눈으로 맞추면,
    히스토그램의 꼬리가 상자그림에서 어떻게 "이상치 점"으로 바뀌는지 보인다.
    """
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)
    minor_1 = stats.norm(scale=2).rvs(200)
    minor_2 = stats.norm(scale=4).rvs(100)
    combined = np.concatenate((main_data, minor_1, minor_2))

    fig, (ax_hist, ax_box) = plt.subplots(2, 1, figsize=(12, 6))
    ax_hist.hist(combined, density=True, bins=30)
    ax_hist.set_title("합친 자료의 히스토그램 (밀도)")
    ax_box.boxplot(combined, vert=False)          # vert=False 로 눕혀 위 그림과 축을 맞춘다
    ax_box.set_title("합친 자료의 상자그림")
    plt.tight_layout()
    plt.show()

    # 상자그림이 이상치로 찍는 점이 몇 개인지 세어 본다.
    # 꼬리가 두꺼우면 1.5*IQR 울타리 밖의 점이 많아진다.
    q1, q3 = np.percentile(combined, [25, 75])
    iqr = q3 - q1
    out = ((combined < q1 - 1.5*iqr) | (combined > q3 + 1.5*iqr)).sum()
    print(f"IQR = {iqr:.3f},  울타리 밖 점 {out}개 / {len(combined)}개 "
          f"({out/len(combined):.1%})")
    print("정규분포라면 약 0.7% 이므로, 이보다 많으면 꼬리가 두꺼운 것이다.")

if __name__ == "__main__":
    generate_and_plot_histogram_and_box_plot_mixed_distribution()

IQR = 1.565,  울타리 밖 점 63개 / 1300개 (4.8%)
정규분포라면 약 0.7% 이므로, 이보다 많으면 꼬리가 두꺼운 것이다.

import numpy as np
from scipy import stats
from scipy.optimize import brentq

w = np.array([10, 2, 1]) / 13
s = np.array([1.0, 2.0, 4.0])

G = lambda x: float(np.sum(w * stats.norm.cdf(x / s)))
q3 = brentq(lambda x: G(x) - 0.75, 0.01, 20, xtol=1e-13)
iqr = 2 * q3
fence = 4 * q3
p = 2 * float(np.sum(w * stats.norm.cdf(-fence / s)))
print(f"Q3 = {q3:.6f},  IQR = 2*Q3 = {iqr:.6f},  울타리 = 4*Q3 = {fence:.6f}")
print(f"울타리 밖 확률 p = {p:.6f}  ({p * 100:.3f}%)")

n = 1300
print(f"\nn = {n} 에서 기대 개수 np = {n * p:.2f},  표준편차 {np.sqrt(n * p * (1 - p)):.2f}")
print(f"  코드가 센 개수 63 개는 {(63 - n * p) / np.sqrt(n * p * (1 - p)):+.2f} 표준편차")

qn = stats.norm.ppf(0.75)
pn = 2 * stats.norm.cdf(-4 * qn)
print(f"\n정규분포: Q3 = {qn:.6f},  울타리 {4 * qn:.6f},  밖 확률 {pn:.6f} ({pn * 100:.3f}%)")
print(f"두꺼운꼬리 혼합은 정규의 {p / pn:.2f} 배")
print(f"표본 IQR 1.565 대 모집단 IQR {iqr:.3f}")

Q3 = 0.793607,  IQR = 2*Q3 = 1.587213,  울타리 = 4*Q3 = 3.174427
울타리 밖 확률 p = 0.051336  (5.134%)

n = 1300 에서 기대 개수 np = 66.74,  표준편차 7.96
  코드가 센 개수 63 개는 -0.47 표준편차

정규분포: Q3 = 0.674490,  울타리 2.697959,  밖 확률 0.006977 (0.698%)
두꺼운꼬리 혼합은 정규의 7.36 배
표본 IQR 1.565 대 모집단 IQR 1.587

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

def generate_and_plot_histogram_and_box_plot_right_skewed(seed: int = 0):
    """오른쪽으로 치우친 자료를 히스토그램과 상자그림으로 나란히 본다.

    눕힌 상자그림에서 치우침은 세 가지로 드러난다. 중앙값 선이 상자의
    왼쪽(작은 값 쪽)으로 밀리고, 오른쪽 수염이 왼쪽보다 길어지며,
    이상치 점이 오른쪽에만 찍힌다.
    """
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)
    right_1 = stats.norm(loc=2).rvs(200)
    right_2 = stats.norm(loc=4).rvs(100)
    combined = np.concatenate((main_data, right_1, right_2))

    fig, (ax_hist, ax_box) = plt.subplots(2, 1, figsize=(12, 6))
    ax_hist.hist(combined, density=True, bins=30)
    ax_hist.set_title("합친 자료의 히스토그램 (밀도)")
    ax_box.boxplot(combined, vert=False)      # 눕혀야 위 히스토그램과 x축이 맞는다
    ax_box.set_title("합친 자료의 상자그림")
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    generate_and_plot_histogram_and_box_plot_right_skewed()

import numpy as np
import scipy.stats as stats
from scipy.optimize import brentq

w = np.array([10, 2, 1]) / 13
locs = np.array([0.0, 2.0, 4.0])
F = lambda x: float(np.sum(w * stats.norm.cdf(x, locs, 1.0)))
Q1, Q2, Q3 = (brentq(lambda x: F(x) - p, -8, 14, xtol=1e-13) for p in (0.25, 0.5, 0.75))
IQR = Q3 - Q1
print(f"모집단  Q1={Q1:.6f}  Q2={Q2:.6f}  Q3={Q3:.6f}  IQR={IQR:.6f}")
print(f"  왼쪽 상자 Q2-Q1 = {Q2 - Q1:.6f},  오른쪽 상자 Q3-Q2 = {Q3 - Q2:.6f}")
print(f"  보울리 사분위왜도 = {(Q3 - 2 * Q2 + Q1) / IQR:.6f}   (gamma1 = 0.828618)")

lo, hi = Q1 - 1.5 * IQR, Q3 + 1.5 * IQR
p_lo, p_hi = F(lo), 1 - F(hi)
print(f"\n울타리 [{lo:.6f}, {hi:.6f}]")
print(f"  아래 확률 {p_lo:.6e} -> n=1300 기대 {p_lo * 1300:.2f} 개")
print(f"  위   확률 {p_hi:.6f} -> n=1300 기대 {p_hi * 1300:.2f} 개")

# 표본에서 실제로
np.random.seed(0)
x = np.concatenate((stats.norm().rvs(1_000), stats.norm(loc=2).rvs(200),
                    stats.norm(loc=4).rvs(100)))
q1, q2, q3 = np.percentile(x, [25, 50, 75])
iqr = q3 - q1
l, h = q1 - 1.5 * iqr, q3 + 1.5 * iqr
inl = x[(x >= l) & (x <= h)]
print(f"\n표본    Q1={q1:.4f}  Q2={q2:.4f}  Q3={q3:.4f}  IQR={iqr:.4f}")
print(f"  왼쪽 상자 {q2 - q1:.4f},  오른쪽 상자 {q3 - q2:.4f},  보울리 {(q3 - 2 * q2 + q1) / iqr:.4f}")
print(f"  이상치: 아래 {np.sum(x < l)} 개, 위 {np.sum(x > h)} 개")
print(f"  수염: 왼쪽 {q1 - inl.min():.4f}, 오른쪽 {inl.max() - q3:.4f}")

모집단  Q1=-0.457652  Q2=0.358251  Q3=1.402022  IQR=1.859674
  왼쪽 상자 Q2-Q1 = 0.815903,  오른쪽 상자 Q3-Q2 = 1.043771
  보울리 사분위왜도 = 0.122531   (gamma1 = 0.828618)

울타리 [-3.247164, 4.191534]
  아래 확률 4.483242e-04 -> n=1300 기대 0.58 개
  위   확률 0.034816 -> n=1300 기대 45.26 개

표본    Q1=-0.4800  Q2=0.3139  Q3=1.4115  IQR=1.8915
  왼쪽 상자 0.7939,  오른쪽 상자 1.0975,  보울리 0.1605
  이상치: 아래 0 개, 위 48 개
  수염: 왼쪽 2.5661, 오른쪽 2.8144

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

def generate_and_plot_histogram_and_box_plot_left_skewed(seed: int = 0):
    """왼쪽으로 치우친 자료를 히스토그램과 상자그림으로 나란히 본다.

    앞 보기를 좌우로 뒤집은 것이다. 중앙값 선이 상자의 오른쪽(큰 값 쪽)으로
    밀리고, 왼쪽 수염이 오른쪽보다 길어지며, 이상치 점이 왼쪽에만 찍힌다.
    """
    np.random.seed(seed)
    main_data = stats.norm().rvs(1_000)
    left_1 = stats.norm(loc=-2).rvs(200)
    left_2 = stats.norm(loc=-4).rvs(100)
    combined = np.concatenate((main_data, left_1, left_2))

    fig, (ax_hist, ax_box) = plt.subplots(2, 1, figsize=(12, 6))
    ax_hist.hist(combined, density=True, bins=30)
    ax_hist.set_title("합친 자료의 히스토그램 (밀도)")
    ax_box.boxplot(combined, vert=False)      # 눕혀야 위 히스토그램과 x축이 맞는다
    ax_box.set_title("합친 자료의 상자그림")
    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    generate_and_plot_histogram_and_box_plot_left_skewed()

import numpy as np
import scipy.stats as stats
from scipy.optimize import brentq

w = np.array([10, 2, 1]) / 13

def box(locs, scales):
    """혼합의 사분위수, 보울리 사분위왜도, 울타리 밖 확률을 돌려준다."""
    F = lambda x: float(np.sum(w * stats.norm.cdf(x, locs, scales)))
    q1, q2, q3 = (brentq(lambda x: F(x) - p, -40, 40, xtol=1e-13)
                  for p in (0.25, 0.5, 0.75))
    iqr = q3 - q1
    lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
    return q1, q2, q3, (q3 - 2 * q2 + q1) / iqr, F(lo), 1 - F(hi)

cases = {
    "오른쪽 치우침": (np.array([0.0, 2.0, 4.0]), np.ones(3)),
    "왼쪽 치우침": (np.array([0.0, -2.0, -4.0]), np.ones(3)),
    "대칭 두꺼운꼬리": (np.zeros(3), np.array([1.0, 2.0, 4.0])),
}
print(f"{'':>14}{'Q1':>10}{'Q2':>10}{'Q3':>10}{'IQR':>9}{'보울리':>10}"
      f"{'아래 밖%':>10}{'위 밖%':>9}")
for name, (lo_, sc) in cases.items():
    q1, q2, q3, bw, pl, ph = box(lo_, sc)
    print(f"{name:>14}{q1:>10.4f}{q2:>10.4f}{q3:>10.4f}{q3 - q1:>9.4f}"
          f"{bw:>+10.4f}{pl * 100:>10.3f}{ph * 100:>9.3f}")

# 표본에서 좌우 신호가 뒤집히는가
np.random.seed(0)
y = np.concatenate((stats.norm().rvs(1_000), stats.norm(loc=-2).rvs(200),
                    stats.norm(loc=-4).rvs(100)))
a1, a2, a3 = np.percentile(y, [25, 50, 75])
ir = a3 - a1
l, h = a1 - 1.5 * ir, a3 + 1.5 * ir
print(f"\n왼쪽 치우침 표본: 왼쪽 상자 {a2 - a1:.4f},  오른쪽 상자 {a3 - a2:.4f},"
      f"  보울리 {(a3 - 2 * a2 + a1) / ir:+.4f}")
print(f"  이상치: 아래 {np.sum(y < l)} 개, 위 {np.sum(y > h)} 개")

                      Q1        Q2        Q3      IQR       보울리     아래 밖%     위 밖%
       오른쪽 치우침   -0.4577    0.3583    1.4020   1.8597   +0.1225     0.045    3.482
        왼쪽 치우침   -1.4020   -0.3583    0.4577   1.8597   -0.1225     3.482    0.045
      대칭 두꺼운꼬리   -0.7936    0.0000    0.7936   1.5872   -0.0000     2.567    2.567

왼쪽 치우침 표본: 왼쪽 상자 1.0108,  오른쪽 상자 0.7800,  보울리 -0.1289
  이상치: 아래 44 개, 위 0 개

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

def generate_samples(main_size, right_size, left_size):
    """왼쪽·오른쪽 덩어리의 **개수 차이**로 치우침을 만든다.

    right_size > left_size 이면 오른쪽이 무거워져 양의 왜도가 되고,
    두 값이 같으면 대칭이 된다. 아래 main()에서 개수를 바꿔 가며
    왜도가 어떻게 움직이는지 확인할 수 있다.
    """
    main_sample = np.random.normal(0, 1, main_size)      # 중앙 덩어리
    right_sample = np.random.normal(2, 1, right_size)    # 오른쪽 덩어리
    left_sample = np.random.normal(-2, 1, left_size)     # 왼쪽 덩어리
    return np.concatenate([main_sample, right_sample, left_sample])

def calculate_statistics(data):
    """평균, 표준편차, 왜도를 구한다.

    여기서는 n으로 나누는 모집단 표준편차를 쓴다(정규 밀도를 겹쳐 그리기 위함).
    표본표준편차가 필요하면 n-1로 나눠야 한다.
    """
    n = data.shape[0]
    mean = data.sum() / n
    std_dev = np.sqrt(np.sum((data - mean) ** 2) / n)
    skewness = stats.describe(data).skewness
    return mean, std_dev, skewness

def plot_distribution_with_normal_fit(data, mean, std_dev, skewness, title):
    """히스토그램 위에 같은 평균·표준편차의 정규분포를 겹쳐 그린다.

    두 곡선이 어긋나는 방식이 곧 왜도(또는 첨도)의 시각적 정체다.
    치우친 자료에서는 정규곡선이 봉우리를 지나치고 꼬리 쪽에서 벌어진다.
    """
    fig, ax = plt.subplots(figsize=(12, 3))
    _, bins, _ = ax.hist(data, density=True, bins=100, label="표본")
    normal_pdf = stats.norm(loc=mean, scale=std_dev).pdf(bins)
    ax.plot(bins, normal_pdf, "--r", label="정규분포 밀도")
    ax.set_title(f"{title}\n왜도 = {skewness:.4f}")
    ax.legend()
    ax.spines["left"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)
    plt.show()

def main():
    np.random.seed(0)
    main_size = 10_000
    right_size = 3_000
    left_size = 3_000

    samples = generate_samples(main_size, right_size, left_size)
    mean, std_dev, skewness = calculate_statistics(samples)

    if right_size > left_size:
        title = "오른쪽으로 치우친 분포"
    elif right_size < left_size:
        title = "왼쪽으로 치우친 분포"
    else:
        title = "대칭 분포"

    plot_distribution_with_normal_fit(samples, mean, std_dev, skewness, title)

    print(f"{title}")
    print(f"  평균   {mean:+.4f}")
    print(f"  표준편차 {std_dev:.4f}")
    print(f"  왜도   {skewness:+.4f}")
    print(f"  중앙값 {np.median(samples):+.4f}  (대칭이면 평균과 같아진다)")

if __name__ == "__main__":
    main()

대칭 분포
  평균   -0.0093
  표준편차 1.5661
  왜도   +0.0456
  중앙값 -0.0273  (대칭이면 평균과 같아진다)

import numpy as np
from scipy import stats

w = np.array([10, 3, 3]) / 16
d = np.array([0.0, 2.0, -2.0])          # 평균이 0 이라 중심이 곧 편차다
m2 = np.sum(w * (d ** 2 + 1))
m3 = np.sum(w * (d ** 3 + 3 * d))
m4 = np.sum(w * (d ** 4 + 6 * d ** 2 + 3))
print(f"mu2 = {m2:.6f}  mu3 = {m3:.6f}  mu4 = {m4:.6f}")
print(f"sigma = {np.sqrt(m2):.6f},  gamma1 = {m3 / m2 ** 1.5:.6f},"
      f"  gamma2 = {m4 / m2 ** 2 - 3:.6f}")

# 같은 분포에서 표본을 되풀이 뽑아 g1 의 표집분포를 잰다.
rng = np.random.default_rng(11)
R, n = 2000, 16_000
g = np.empty(R)
for b in range(R):
    x = np.concatenate([rng.normal(0, 1, 10_000), rng.normal(2, 1, 3_000),
                        rng.normal(-2, 1, 3_000)])
    g[b] = stats.skew(x)
se = g.std()
print(f"\n반복 {R}회,  n = {n}")
print(f"  g1 의 평균 {g.mean():+.5f}  (참값 0,  몬테카를로 오차 {se / np.sqrt(R):.5f})")
print(f"  g1 의 표준편차(= 표준오차) {se:.5f}   정규 어림 sqrt(6/n) = {np.sqrt(6 / n):.5f}")
print(f"  관측값 +0.0456 은 {0.0456 / se:.2f} 표준오차")
print(f"  |g1| >= 0.0456 일 확률 {np.mean(np.abs(g) >= 0.0456):.4f}")

# 왜 양수가 나왔는가 -- 덩어리별 표본평균을 본다.
np.random.seed(0)
main = np.random.normal(0, 1, 10_000)
right = np.random.normal(2, 1, 3_000)
left = np.random.normal(-2, 1, 3_000)
print(f"\n덩어리별 표본평균: 본체 {main.mean():+.4f}, 오른쪽 {right.mean():+.4f}"
      f" (참값 +2), 왼쪽 {left.mean():+.4f} (참값 -2)")
print(f"  두 덩어리의 중심이 참값보다 각각 {right.mean() - 2:+.4f}, {left.mean() + 2:+.4f}"
      f" 만큼 치우쳐 있다")

mu2 = 2.500000  mu3 = 0.000000  mu4 = 18.000000
sigma = 1.581139,  gamma1 = 0.000000,  gamma2 = -0.120000

반복 2000회,  n = 16000
  g1 의 평균 +0.00046  (참값 0,  몬테카를로 오차 0.00037)
  g1 의 표준편차(= 표준오차) 0.01646   정규 어림 sqrt(6/n) = 0.01936
  관측값 +0.0456 은 2.77 표준오차
  |g1| >= 0.0456 일 확률 0.0045

덩어리별 표본평균: 본체 -0.0184, 오른쪽 +1.9881 (참값 +2), 왼쪽 -1.9763 (참값 -2)
  두 덩어리의 중심이 참값보다 각각 -0.0119, +0.0237 만큼 치우쳐 있다

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

def generate_samples(main_size, peak_size):
    """중앙에 아주 좁은(표준편차 0.2) 덩어리를 얹어 봉우리를 뾰족하게 만든다.

    중심은 둘 다 0이므로 대칭은 유지되고, 봉우리만 솟는다.
    이것이 급첨(leptokurtic) 분포를 만드는 가장 간단한 방법이다.
    """
    main_sample = np.random.normal(0, 1, main_size)      # 넓은 본체
    peak_sample = np.random.normal(0, 0.2, peak_size)    # 좁고 뾰족한 봉우리
    return np.concatenate([main_sample, peak_sample])

def calculate_statistics(data):
    """첨도를 정의 그대로 계산한다.

    표준화한 값의 네제곱 평균이 첨도다. 네제곱이므로 중심에서 멀리 떨어진
    값이 압도적으로 큰 기여를 한다. 첨도가 사실상 "꼬리의 무게"를 재는 이유다.
    """
    mean = np.mean(data)
    std_dev = np.std(data)
    skewness = stats.describe(data).skewness
    kurtosis = np.mean(((data - mean) / std_dev) ** 4)
    excess_kurtosis = kurtosis - 3     # 정규분포의 첨도 3을 빼면 초과첨도
    return mean, std_dev, skewness, kurtosis, excess_kurtosis

def plot_distribution_with_normal_fit(data, mean, std_dev, excess_kurtosis, title):
    fig, ax = plt.subplots(figsize=(12, 3))
    _, bins, _ = ax.hist(data, density=True, bins=100, label="표본")
    normal_pdf = stats.norm(loc=mean, scale=std_dev).pdf(bins)
    ax.plot(bins, normal_pdf, "--r", label="정규분포 밀도")
    ax.set_title(f"{title}\n초과첨도 = {excess_kurtosis:.4f}")
    ax.legend()
    ax.spines["left"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.spines["top"].set_visible(False)
    plt.show()

def main():
    np.random.seed(0)
    main_size = 10_000
    peak_size = 500

    data = generate_samples(main_size, peak_size)
    mean, std_dev, skewness, kurtosis, excess_kurtosis = calculate_statistics(data)

    if excess_kurtosis > 0:
        title = "급첨 분포"
    elif excess_kurtosis < 0:
        title = "평첨 분포"
    else:
        title = "중첨 분포"

    plot_distribution_with_normal_fit(data, mean, std_dev, excess_kurtosis, title)

    print(f"{title}")
    print(f"  왜도       {skewness:+.4f}  (좌우 대칭이므로 0 근처)")
    print(f"  첨도       {kurtosis:.4f}")
    print(f"  초과첨도   {excess_kurtosis:+.4f}  (정규분포는 0)")

if __name__ == "__main__":
    main()

급첨 분포
  왜도       +0.0231  (좌우 대칭이므로 0 근처)
  첨도       3.1053
  초과첨도   +0.1053  (정규분포는 0)

import numpy as np


def mixture_kurtosis(weights, scales):
    """중심이 같은 정규혼합의 초과첨도. 보기 1 에서 유도한 식 그대로다."""
    w = np.asarray(weights, dtype=float)
    w = w / w.sum()
    V = np.asarray(scales, dtype=float) ** 2       # 성분 분산
    EV = np.sum(w * V)
    VarV = np.sum(w * (V - EV) ** 2)
    return 3 * VarV / EV ** 2


print(f"보기 8 (10000 @ 1, 500 @ 0.2) : gamma2 = {mixture_kurtosis([10000, 500], [1, 0.2]):.6f}")
print(f"보기 1 (1000, 200, 100 @ 1,2,4): gamma2 = {mixture_kurtosis([1000, 200, 100], [1, 2, 4]):.6f}")

print(f"\n봉우리를 좁게 해도 첨도는 멈춘다")
print(f"{'봉우리 sd':>10}{'gamma2':>12}")
for sd in (0.5, 0.2, 0.1, 0.01, 0.0):
    print(f"{sd:>10.2f}{mixture_kurtosis([10000, 500], [1, sd]):>12.6f}")

print(f"\n같은 가중치로 바깥 성분을 넓히면")
print(f"{'바깥 sd':>10}{'gamma2':>12}")
for sd in (2, 4, 10, 50, 200):
    print(f"{sd:>10}{mixture_kurtosis([10000, 500], [1, sd]):>12.4f}")

보기 8 (10000 @ 1, 500 @ 0.2) : gamma2 = 0.137689
보기 1 (1000, 200, 100 @ 1,2,4): gamma2 = 7.053633

봉우리를 좁게 해도 첨도는 멈춘다
    봉우리 sd      gamma2
      0.50    0.082305
      0.20    0.137689
      0.10    0.146868
      0.01    0.149969
      0.00    0.150000

같은 가중치로 바깥 성분을 넓히면
     바깥 sd      gamma2
         2      0.9375
         4     10.4167
        10     40.8375
        50     59.0042
       200     59.9370

from scipy import stats
import numpy as np

np.random.seed(0)                       # 시드를 고정해야 아래 출력이 재현된다
data = np.random.normal(0, 1, 10000)    # 정규분포이므로 초과첨도의 참값은 0

# 두 함수 모두 "초과첨도"(첨도 - 3)를 돌려준다.
# 즉 정규분포에서 0이 나오도록 이미 3을 빼 놓았다.
# 3을 빼지 않은 값이 필요하면 fisher=False 를 준다.
print(stats.kurtosis(data))
print(stats.describe(data).kurtosis)
print(stats.kurtosis(data, fisher=False), "  <- 3을 빼지 않은 값")

-0.03095451095565238
-0.03095451095565238
2.9690454890443476   <- 3을 빼지 않은 값

import numpy as np
from scipy import stats

np.random.seed(0)
data = np.random.normal(0, 1, 10000)
n = len(data)

g2 = stats.kurtosis(data)                       # 기본: fisher=True, bias=True
G2 = stats.kurtosis(data, bias=False)           # 불편보정
print(f"g2 (bias=True ) = {g2:.10f}")
print(f"G2 (bias=False) = {G2:.10f}")
print(f"공식 (n-1)/((n-2)(n-3)) * ((n+1)g2 + 6) = "
      f"{(n - 1) / ((n - 2) * (n - 3)) * ((n + 1) * g2 + 6):.10f}")
print(f"왜도도 같다: g1 = {stats.skew(data):.10f}, G1 = {stats.skew(data, bias=False):.10f},"
      f" 공식 {np.sqrt(n * (n - 1)) / (n - 2) * stats.skew(data):.10f}")

# 표본첨도의 천장. 한 점만 멀리 둔 배열이 천장을 이룬다.
print(f"\n{'n':>5}{'천장 b2':>12}{'달성한 b2':>12}{'천장(초과)':>12}{'달성(초과)':>12}")
for m in (5, 8, 12, 20, 100):
    x = np.zeros(m)
    x[-1] = 1e9
    ceil = (m * m - 3 * m + 3) / (m - 1)
    print(f"{m:>5}{ceil:>12.6f}{np.mean(((x - x.mean()) / x.std()) ** 4):>12.6f}"
          f"{ceil - 3:>12.6f}{stats.kurtosis(x):>12.6f}")

print(f"\n보기 1 의 참 초과첨도 7.053633 을 재려면 천장이 그보다 커야 한다")
for m in (8, 10, 11, 12):
    print(f"  n = {m:>3}:  천장(초과) {(m * m - 3 * m + 3) / (m - 1) - 3:>8.4f}"
          f"   7.053633 을 담는가: {(m * m - 3 * m + 3) / (m - 1) - 3 >= 7.053633}")

g2 (bias=True ) = -0.0309545110
G2 (bias=False) = -0.0303697537
공식 (n-1)/((n-2)(n-3)) * ((n+1)g2 + 6) = -0.0303697537
왜도도 같다: g1 = 0.0266346167, G1 = 0.0266386127, 공식 0.0266386127

    n       천장 b2      달성한 b2      천장(초과)      달성(초과)
    5    3.250000    3.250000    0.250000    0.250000
    8    6.142857    6.142857    3.142857    3.142857
   12   10.090909   10.090909    7.090909    7.090909
   20   18.052632   18.052632   15.052632   15.052632
  100   98.010101   98.010101   95.010101   95.010101

보기 1 의 참 초과첨도 7.053633 을 재려면 천장이 그보다 커야 한다
  n =   8:  천장(초과)   3.1429   7.053633 을 담는가: False
  n =  10:  천장(초과)   5.1111   7.053633 을 담는가: False
  n =  11:  천장(초과)   6.1000   7.053633 을 담는가: False
  n =  12:  천장(초과)   7.0909   7.053633 을 담는가: True

import numpy as np
from scipy import stats

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

normal = rng.normal(0, 1, n)
mixed = np.where(rng.random(n) < 0.95, rng.normal(0, 0.6, n), rng.normal(0, 2.6, n))
mixed = mixed / mixed.std()                      # 분산을 1로 맞춘다
uniform = rng.uniform(-np.sqrt(3), np.sqrt(3), n)

print(f"{'분포':>10}{'분산':>9}{'초과첨도':>11}{'|X|<0.5 비율':>15}{'|X|>3 비율':>14}")
for label, x in [("정규", normal), ("혼합", mixed), ("균등", uniform)]:
    print(f"{label:>10}{x.var():>9.4f}{stats.kurtosis(x):>11.4f}"
          f"{np.mean(np.abs(x) < 0.5):>15.4f}{np.mean(np.abs(x) > 3):>14.5f}")

        분포       분산       초과첨도     |X|<0.5 비율      |X|>3 비율
        정규   1.0028    -0.0013         0.3816       0.00280
        혼합   1.0000    12.9937         0.4908       0.01710
        균등   0.9986    -1.1986         0.2891       0.00000

import numpy as np
from scipy import stats

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

def bowley(x):
    q1, q2, q3 = np.quantile(x, [0.25, 0.5, 0.75])
    return (q3 + q1 - 2 * q2) / (q3 - q1)

def moors(x):
    e = np.quantile(x, [0.125, 0.25, 0.375, 0.625, 0.75, 0.875])
    return ((e[5] - e[3]) + (e[2] - e[0])) / (e[4] - e[1])

print(f"{'분포':>8}{'고전 왜도':>12}{'보울리':>10}{'초과첨도':>12}{'무어스':>10}")
for label, x in [("정규", rng.normal(0, 1, n)), ("지수", rng.exponential(1, n)),
                 ("t(5)", rng.standard_t(5, n)), ("균등", rng.uniform(-1, 1, n))]:
    print(f"{label:>8}{stats.skew(x):>12.4f}{bowley(x):>10.4f}"
          f"{stats.kurtosis(x):>12.4f}{moors(x):>10.4f}")
print("\n(정규분포의 무어스 첨도 기준값은 1.2331)")

print("\n이상치 하나에 대한 민감도")
base = rng.normal(0, 1, 1000)
for v in (0, 10, 30):
    x = np.append(base, v) if v else base
    print(f"  극단값 {v:>3} 추가: 고전 왜도 {stats.skew(x):>8.4f}  보울리 {bowley(x):>7.4f}"
          f"   초과첨도 {stats.kurtosis(x):>9.4f}  무어스 {moors(x):.4f}")

      분포       고전 왜도       보울리        초과첨도       무어스
      정규     -0.0053    0.0031      0.0061    1.2293
      지수      2.0017    0.2618      5.9660    1.2995
    t(5)      0.0248   -0.0030      9.7070    1.3262
      균등      0.0012   -0.0001     -1.1998    0.9989

(정규분포의 무어스 첨도 기준값은 1.2331)

이상치 하나에 대한 민감도
  극단값   0 추가: 고전 왜도  -0.1205  보울리  0.0119   초과첨도   -0.0097  무어스 1.2893
  극단값  10 추가: 고전 왜도   0.7166  보울리  0.0115   초과첨도    7.5040  무어스 1.2875
  극단값  30 추가: 고전 왜도  10.1092  보울리  0.0115   초과첨도  218.8426  무어스 1.2875

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
cases = {
    '정규':          rng.normal(0, 1, 2_000_000),
    '지수':          rng.exponential(1, 2_000_000),
    '균등':          rng.uniform(0, 1, 2_000_000),
    '베르누이(0.5)': rng.binomial(1, 0.5, 2_000_000).astype(float),
    '베르누이(0.1)': rng.binomial(1, 0.1, 2_000_000).astype(float),
    '로그정규':      rng.lognormal(0, 1, 2_000_000),
}
print(f"{'분포':>16}{'왜도':>10}{'첨도':>10}{'왜도^2+1':>12}{'여유':>10}")
for k, v in cases.items():
    s = stats.skew(v)
    ku = stats.kurtosis(v, fisher=False)
    print(f"{k:>16}{s:>10.4f}{ku:>10.4f}{s**2+1:>12.4f}{ku-(s**2+1):>10.4f}")