사분위범위와 강건 측도¶

import numpy as np
import pandas as pd

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
loans_data = pd.read_csv(url)

# 범위 = 최댓값 - 최솟값. 자료 전체에서 딱 두 점만 쓴다.
data_range = loans_data['x'].max() - loans_data['x'].min()
print(f"{data_range = }")

# 그 두 점이 무엇인지도 함께 보자. 범위가 왜 취약한지 바로 드러난다.
print(f"최솟값 {loans_data['x'].min():,}  최댓값 {loans_data['x'].max():,}")
print(f"관측값 {len(loans_data):,}개 중 단 2개가 이 값을 정한다")

x = loans_data['x'].astype(float).values
srt = np.sort(x)
N = len(x)

# 끝점 하나를 지우면 범위가 어떻게 되는가. 최솟값과 최댓값의 사정이 다르다.
print(f"\n최솟값 {srt[0]:,.0f} 은 {int((x == srt[0]).sum())}개, "
      f"최댓값 {srt[-1]:,.0f} 은 {int((x == srt[-1]).sum())}개")
print(f"  최솟값 하나를 지우면 범위 {srt[-1] - srt[1]:,.0f}  ({srt[-1] - srt[1] - data_range:+,.0f})")
print(f"  최댓값 하나를 지우면 범위 {srt[-1] - srt[0]:,.0f}  (변화 없다 — 같은 값이 둘이다)")

# (1) 부분표본 크기 m 을 키우면 범위가 자란다. IQR 은 자라지 않는다.
rng = np.random.default_rng(0)
B = 200
print(f"\n부분표본 {B}회 추출 (비복원)")
print(f"{'m':>7}{'범위 평균':>13}{'범위 sd':>11}{'IQR 평균':>12}{'IQR sd':>10}")
for m in (100, 1000, 10000, 50000):
    R, I = [], []
    for _ in range(B):
        s = rng.choice(x, m, replace=False)
        R.append(s.max() - s.min())
        I.append(np.subtract(*np.percentile(s, [75, 25])))
    print(f"{m:>7}{np.mean(R):>13,.0f}{np.std(R, ddof=1):>11,.0f}"
          f"{np.mean(I):>12,.0f}{np.std(I, ddof=1):>10,.0f}")

# m=100 에서 범위가 전체 범위와 같아질 확률을 정확히 센다.
from math import comb
m = 100
p_min = m / N                                   # 최솟값 1개가 뽑힐 확률
p_max = 1 - comb(N - 2, m) / comb(N, m)         # 최댓값 2개 중 하나라도 뽑힐 확률
print(f"\nm=100 에서  P(최솟값 포함) = {p_min:.5f},  P(최댓값 중 하나 포함) = {p_max:.5f}")
print(f"  P(범위 = 195,000) = 두 사건의 곱에 가까운 {p_min * p_max:.3e}")

data_range = 195000
최솟값 4,000  최댓값 199,000
관측값 50,000개 중 단 2개가 이 값을 정한다

최솟값 4,000 은 1개, 최댓값 199,000 은 2개
  최솟값 하나를 지우면 범위 192,100  (-2,900)
  최댓값 하나를 지우면 범위 195,000  (변화 없다 — 같은 값이 둘이다)

부분표본 200회 추출 (비복원)
      m        범위 평균      범위 sd      IQR 평균    IQR sd
    100      157,818     15,374      40,331     5,112
   1000      183,628      4,335      41,029     1,720
  10000      191,544      1,806      40,320       528
  50000      195,000          0      40,000         0

m=100 에서  P(최솟값 포함) = 0.00200,  P(최댓값 중 하나 포함) = 0.00400
  P(범위 = 195,000) = 두 사건의 곱에 가까운 7.992e-06

import matplotlib
matplotlib.use("Agg")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
df = pd.read_csv(url)

# 같은 자료를 두 짝의 측도로 요약한다.
#   비강건한 짝: 평균 ± 표준편차   (모든 관측값을 다 쓴다)
#   강건한 짝  : 중앙값, Q1, Q3    (순위만 쓴다)
mean_income = df['x'].mean()
median_income = df['x'].median()
std_dev = df['x'].std()
q1 = df['x'].quantile(0.25)
q3 = df['x'].quantile(0.75)
iqr = stats.iqr(df['x'])

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(12, 4))

# 왼쪽: 평균 ± 표준편차.
# 소득은 오른쪽으로 치우쳐 있어 평균이 봉우리보다 오른쪽에 놓이고,
# "평균 - 표준편차"가 자료가 별로 없는 곳을 가리킨다.
ax1.hist(df['x'], bins=30, density=True, color="#DCEBFB", edgecolor="white")
ax1.axvline(mean_income, color="#1565C0", linestyle='--', lw=2, label="평균")
ax1.axvline(mean_income - std_dev, color="#D32F2F", linestyle='--', lw=2,
            label="평균 - 표준편차")
ax1.axvline(mean_income + std_dev, color="#E65100", linestyle='--', lw=2,
            label="평균 + 표준편차")
ax1.legend(fontsize=8)
ax1.set_title("평균과 표준편차")
ax1.set_xlabel("소득 (달러)")
ax1.set_ylabel("밀도")

# 가운데: 중앙값과 사분위수.
# Q1과 Q3 사이가 정확히 자료의 가운데 50%이며, 치우침에 흔들리지 않는다.
ax2.hist(df['x'], bins=30, density=True, color="#DCEBFB", edgecolor="white")
ax2.axvline(median_income, color="#1565C0", linestyle='--', lw=2, label="중앙값")
ax2.axvline(q1, color="#D32F2F", linestyle='--', lw=2, label="$Q_1$")
ax2.axvline(q3, color="#E65100", linestyle='--', lw=2, label="$Q_3$")
ax2.legend(fontsize=8)
ax2.set_title("중앙값과 사분위수")
ax2.set_xlabel("소득 (달러)")
ax2.set_ylabel("밀도")

# 오른쪽: 같은 사분위수를 상자그림으로 옮긴 것
ax3.boxplot(df['x'], vert=True, patch_artist=True, labels=["소득"])
ax3.set_title("상자그림")
ax3.set_ylabel("소득 (달러)")

fig.tight_layout()
fig.savefig("robust_63.png", dpi=170, facecolor="white", bbox_inches="tight")

# 두 짝의 숫자를 나란히 찍어 비교한다
print(f"평균   {mean_income:>9,.0f}   표준편차 {std_dev:>9,.0f}")
print(f"중앙값 {median_income:>9,.0f}   IQR      {iqr:>9,.0f}")
print(f"평균 - 표준편차 = {mean_income - std_dev:>9,.0f}   "
      f"(최솟값 {df['x'].min():,}보다 큰가? "
      f"{'예' if mean_income - std_dev > df['x'].min() else '아니오'})")

# 그림의 눈금이 실제로 어디를 가리키는지 센다.
v = df['x'].astype(float).values
print(f"\n{'눈금':>12}{'값':>12}{'그 아래 자료 비율':>18}")
for name, val in (("평균 - s", mean_income - std_dev), ("Q1", q1),
                  ("평균", mean_income), ("중앙값", median_income),
                  ("Q3", q3), ("평균 + s", mean_income + std_dev)):
    print(f"{name:>12}{val:>12,.0f}{np.mean(v < val):>18.4f}")
print(f"\n[평균-s, 평균+s] 안의 비율 {np.mean((v >= mean_income - std_dev) & (v <= mean_income + std_dev)):.4f}"
      f"  (정규라면 0.6827)")
print(f"[Q1, Q3]          안의 비율 {np.mean((v >= q1) & (v <= q3)):.4f}  (정의상 0.5)")
print(f"\n왜도 {stats.skew(v):.4f}  초과첨도 {stats.kurtosis(v):.4f}")
print(f"IQR/1.349 = {iqr / 1.349:,.0f}  대  s = {std_dev:,.0f}   비 {std_dev / (iqr / 1.349):.4f}")
print(f"45,000 인 관측값 {int((v == 45000).sum()):,}개, 85,000 인 관측값 {int((v == 85000).sum()):,}개")

평균      68,761   표준편차    32,872
중앙값    62,000   IQR         40,000
평균 - 표준편차 =    35,888   (최솟값 4,000보다 큰가? 예)

          눈금           값        그 아래 자료 비율
      평균 - s      35,888            0.1275
          Q1      45,000            0.2423
          평균      68,761            0.5762
         중앙값      62,000            0.4976
          Q3      85,000            0.7350
      평균 + s     101,633            0.8544

[평균-s, 평균+s] 안의 비율 0.7269  (정규라면 0.6827)
[Q1, Q3]          안의 비율 0.5100  (정의상 0.5)

왜도 1.0488  초과첨도 1.0808
IQR/1.349 = 29,652  대  s = 32,872   비 1.1086
45,000 인 관측값 1,300개, 85,000 인 관측값 863개

import numpy as np
import pandas as pd

url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
df = pd.read_csv(url)

# quantile(p)는 자료의 p 비율이 그 아래에 놓이는 값을 돌려준다.
q1 = df['x'].quantile(0.25)      # 아래에서 25%
q2 = df['x'].median()            # 아래에서 50% = 중앙값
q3 = df['x'].quantile(0.75)      # 아래에서 75%

print(f"{q1 = }")
print(f"{q2 = }")  # 중앙값
print(f"{q3 = }")

# IQR은 Q3 - Q1. 자료의 가운데 절반이 차지하는 폭이다.
print(f"IQR = {q3 - q1:,.0f}")

# (1) linear 이 읽는 위치와 그 양옆의 순서통계량
x = df['x'].astype(float).values
srt = np.sort(x)
n = len(x)
for p in (0.25, 0.75):
    h = (n - 1) * p
    lo, hi = int(np.floor(h)), int(np.floor(h)) + 1
    print(f"\np={p}: h=(n-1)p={h}")
    print(f"  x_({lo+1}) = {srt[lo]:,.0f},  x_({hi+1}) = {srt[hi]:,.0f}"
          f"   {'같다 -> 보간해도 그대로' if srt[lo] == srt[hi] else '다르다 -> 보간이 값을 바꾼다'}")

# (2) 보간법을 바꿔 본다
print(f"\n{'method':>10}{'Q1':>12}{'Q3':>12}{'IQR':>12}")
for meth in ("linear", "lower", "higher", "midpoint", "nearest"):
    a, b = np.percentile(x, [25, 75], method=meth)
    print(f"{meth:>10}{a:>12,.0f}{b:>12,.0f}{b - a:>12,.0f}")

print(f"\n고유값 {len(np.unique(x)):,}개 / 관측값 {n:,}개")
print(f"  45,000 인 관측값 {int((x == 45000).sum()):,}개,"
      f"  85,000 인 관측값 {int((x == 85000).sum()):,}개")

q1 = 45000.0
q2 = 62000.0
q3 = 85000.0
IQR = 40,000

p=0.25: h=(n-1)p=12499.75
  x_(12500) = 45,000,  x_(12501) = 45,000   같다 -> 보간해도 그대로

p=0.75: h=(n-1)p=37499.25
  x_(37500) = 85,000,  x_(37501) = 85,000   같다 -> 보간해도 그대로

    method          Q1          Q3         IQR
    linear      45,000      85,000      40,000
     lower      45,000      85,000      40,000
    higher      45,000      85,000      40,000
  midpoint      45,000      85,000      40,000
   nearest      45,000      85,000      40,000

고유값 5,271개 / 관측값 50,000개
  45,000 인 관측값 1,300개,  85,000 인 관측값 863개

import numpy as np

print(np.percentile([2, 4, 5, 7, 8, 9, 11, 13, 15, 80], [25, 75]))

[ 5.5 12.5]

import numpy as np

rng = np.random.default_rng(0)
B, n = 100_000, 20
X = rng.normal(0, 1, (B, n))

estimators = {
    "표본표준편차 s": X.std(axis=1, ddof=1),
    "IQR / 1.349":  np.subtract(*np.percentile(X, [75, 25], axis=1)) / 1.349,
    "MAD x 1.4826": 1.4826 * np.median(
        np.abs(X - np.median(X, axis=1, keepdims=True)), axis=1),
    "범위 / 3.735":  (X.max(1) - X.min(1)) / 3.735,
}

base = None
print(f"정규 N(0,1), n={n}")
print(f"{'추정량':<16}{'평균':>9}{'표준편차':>11}{'MSE':>11}{'상대효율':>11}")
for name, v in estimators.items():
    mse = np.mean((v - 1) ** 2)
    base = base or mse
    print(f"{name:<16}{v.mean():>9.4f}{v.std():>11.4f}{mse:>11.5f}{base / mse:>11.4f}")

정규 N(0,1), n=20
추정량                    평균       표준편차        MSE       상대효율
표본표준편차 s           0.9867     0.1611    0.02614     1.0000
IQR / 1.349        0.9320     0.2383    0.06140     0.4258
MAD x 1.4826       0.9583     0.2495    0.06399     0.4085
범위 / 3.735         1.0002     0.1953    0.03815     0.6853

import numpy as np

rng = np.random.default_rng(0)
d2 = {2: 1.128, 5: 2.326, 10: 3.078, 20: 3.735, 50: 4.498, 100: 5.015, 1000: 6.483}

print(f"{'n':>6}{'E[범위]':>11}{'범위 추정 MSE':>16}{'s 의 MSE':>12}{'범위의 상대효율':>17}")
for n, c in d2.items():
    X = rng.normal(0, 1, (60_000, n))
    r = (X.max(1) - X.min(1)) / c
    s = X.std(axis=1, ddof=1)
    print(f"{n:>6}{(X.max(1) - X.min(1)).mean():>11.4f}{np.mean((r - 1) ** 2):>16.5f}"
          f"{np.mean((s - 1) ** 2):>12.5f}{np.mean((s - 1) ** 2) / np.mean((r - 1) ** 2):>17.4f}")

     n      E[범위]       범위 추정 MSE     s 의 MSE         범위의 상대효율
     2     1.1334         0.57777     0.40699           0.7044
     5     2.3293         0.13812     0.11956           0.8657
    10     3.0792         0.06727     0.05507           0.8187
    20     3.7332         0.03785     0.02604           0.6881
    50     4.4999         0.02106     0.01019           0.4837
   100     5.0139         0.01466     0.00509           0.3471
  1000     6.4842         0.00591     0.00051           0.0857

import numpy as np

rng = np.random.default_rng(1)
n = 100

def trimmed_sd(x, p):
    k = int(len(x) * p)
    return np.sort(x)[k:len(x) - k].std(ddof=1)

def winsorized_sd(x, p):
    lo, hi = np.quantile(x, [p, 1 - p])
    return np.clip(x, lo, hi).std(ddof=1)

clean = rng.normal(0, 1, (4000, n))
dirty = clean.copy()
dirty[:, :10] = rng.normal(0, 6, (4000, 10))       # 10% 오염

print(f"{'추정량':<15}{'정규 자료':>12}{'(표준편차)':>12}{'10% 오염':>12}")
for name, f in [("s", lambda x: x.std(ddof=1)),
                ("10% 절단 sd", lambda x: trimmed_sd(x, 0.10)),
                ("10% 윈저 sd", lambda x: winsorized_sd(x, 0.10)),
                ("MAD x 1.4826", lambda x: 1.4826 * np.median(np.abs(x - np.median(x))))]:
    a = np.array([f(r) for r in clean])
    b = np.array([f(r) for r in dirty])
    print(f"{name:<15}{a.mean():>12.4f}{a.std():>12.4f}{b.mean():>12.4f}")

추정량                   정규 자료      (표준편차)      10% 오염
s                    0.9954      0.0707      2.0875
10% 절단 sd            0.6649      0.0594      0.7517
10% 윈저 sd            0.8186      0.0692      0.9435
MAD x 1.4826         0.9886      0.1164      1.0951

import numpy as np

rng = np.random.default_rng(1)
for s in (0.5, 1.0, 1.5):
    x = rng.lognormal(0, s, 2_000_000)
    q1, q2, q3 = np.quantile(x, [0.25, 0.5, 0.75])
    print(f"lognormal(0, {s}):")
    print(f"  평균 {x.mean():>8.4f}   표준편차 {x.std():>8.4f}   평균 ± sd = "
          f"[{x.mean() - x.std():>7.4f}, {x.mean() + x.std():>7.4f}]")
    print(f"  중앙값 {q2:>6.4f}   IQR {q3 - q1:>8.4f}   [Q1, Q3] = "
          f"[{q1:>7.4f}, {q3:>7.4f}]")
    print(f"  실제로 [평균-sd, 평균+sd] 안에 있는 비율 "
          f"{np.mean((x > x.mean() - x.std()) & (x < x.mean() + x.std())):.4f}")

lognormal(0, 0.5):
  평균   1.1333   표준편차   0.6036   평균 ± sd = [ 0.5297,  1.7369]
  중앙값 1.0005   IQR   0.6861   [Q1, Q3] = [ 0.7142,  1.4003]
  실제로 [평균-sd, 평균+sd] 안에 있는 비율 0.7639
lognormal(0, 1.0):
  평균   1.6513   표준편차   2.1658   평균 ± sd = [-0.5145,  3.8171]
  중앙값 1.0011   IQR   1.4572   [Q1, Q3] = [ 0.5102,  1.9674]
  실제로 [평균-sd, 평균+sd] 안에 있는 비율 0.9095
lognormal(0, 1.5):
  평균   3.0806   표준편차   8.9507   평균 ± sd = [-5.8700, 12.0313]
  중앙값 0.9989   IQR   2.3850   [Q1, Q3] = [ 0.3637,  2.7488]
  실제로 [평균-sd, 평균+sd] 안에 있는 비율 0.9514