중앙값 절대편차 (MAD)¶

import numpy as np
import pandas as pd
from statsmodels import robust

# 미국 50개 주의 인구와 살인율. 오른쪽으로 크게 치우친 전형적인 자료다.
url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)

# 표준편차는 제곱을 쓰므로 멀리 떨어진 값 하나에 크게 흔들린다.
std_dev = state['Population'].std()
print(f"표준편차     : {std_dev:,.0f}")

# statsmodels 의 mad 는 정규분포에서 표준편차와 눈금이 맞도록 이미 보정해 준다.
mad = robust.scale.mad(state['Population'])
print(f"MAD (보정)   : {mad:,.0f}")

# 정의대로 직접 구해 본다. 중앙값에서의 절대편차, 그 중앙값이다.
median_pop = state['Population'].median()
abs_deviations = abs(state['Population'] - median_pop)
mad_manual = abs_deviations.median()
# 0.6745 는 표준정규의 0.75 분위점이다. 이 값으로 나누면 위 mad 와 눈금이 맞는다.
mad_standardized = mad_manual / 0.6744897501960817
print(f"MAD (직접)   : {mad_standardized:,.0f}")
print(f"보정 전 MAD  : {mad_manual:,.0f}")
print(f"s / MAD(보정) = {std_dev / mad:.4f}   <- 정규자료라면 1 쯤이어야 한다")

# --- 최댓값 하나를 빼면 두 척도가 각각 얼마나 움직이는가 ---
pop = state['Population'].astype(float)
n = len(pop)
xbar = pop.mean()
SS = ((pop - xbar) ** 2).sum()
k = pop.idxmax()
xk = pop[k]
print(f"\n최댓값 {state.loc[k, 'State']} = {xk:,.0f},  평균 {xbar:,.1f}")
print(f"  그 한 점이 제곱합에서 차지하는 몫 = {n / (n - 1) * (xk - xbar) ** 2 / SS:.4f}")

# 제곱합 갱신식 SS' = SS - n/(n-1) (x_k - xbar)^2 으로 s' 를 미리 계산한다.
SS_new = SS - n / (n - 1) * (xk - xbar) ** 2
s_pred = np.sqrt(SS_new / (n - 2))
pop2 = pop.drop(k)
print(f"  SS  = {SS:.6e} -> SS' 예측 {SS_new:.6e},  실제 {((pop2 - pop2.mean()) ** 2).sum():.6e}")
print(f"  s   = {std_dev:,.2f} -> s'  예측 {s_pred:,.2f},  실제 {pop2.std():,.2f}"
      f"   ({100 * (s_pred / std_dev - 1):+.2f}%)")
mad2 = robust.scale.mad(pop2)
print(f"  MAD = {mad:,.2f} -> MAD' 실제 {mad2:,.2f}   ({100 * (mad2 / mad - 1):+.2f}%)")

# MAD 가 왜 거의 안 움직이는지는 순위만 보면 안다.
srt, srt2 = np.sort(pop.values), np.sort(pop2.values)
d1 = np.sort(np.abs(pop.values - np.median(pop.values)))
d2 = np.sort(np.abs(pop2.values - np.median(pop2.values)))
print(f"\n  n=50 중앙값 = (x_(25)+x_(26))/2 = ({srt[24]:,.0f}+{srt[25]:,.0f})/2 = {np.median(pop):,.1f}")
print(f"  n=49 중앙값 = x_(25)             = {srt2[24]:,.0f}")
print(f"  n=50 MAD(보정 전) = (d_(25)+d_(26))/2 = "
      f"({d1[24]:,.0f}+{d1[25]:,.0f})/2 = {(d1[24] + d1[25]) / 2:,.1f}")
print(f"  n=49 MAD(보정 전) = d_(25)            = {d2[24]:,.0f}")

표준편차     : 6,848,235
MAD (보정)   : 3,849,876
MAD (직접)   : 3,849,876
보정 전 MAD  : 2,596,702
s / MAD(보정) = 1.7788   <- 정규자료라면 1 쯤이어야 한다

최댓값 California = 37,253,956,  평균 6,162,876.3
  그 한 점이 제곱합에서 차지하는 몫 = 0.4292
  SS  = 2.298018e+15 -> SS' 예측 1.311635e+15,  실제 1.311635e+15
  s   = 6,848,235.35 -> s'  예측 5,227,402.05,  실제 5,227,402.05   (-23.67%)
  MAD = 3,849,876.15 -> MAD' 실제 3,686,302.13   (-4.25%)

  n=50 중앙값 = (x_(25)+x_(26))/2 = (4,339,367+4,533,372)/2 = 4,436,369.5
  n=49 중앙값 = x_(25)             = 4,339,367
  n=50 MAD(보정 전) = (d_(25)+d_(26))/2 = (2,583,376+2,610,028)/2 = 2,596,702.0
  n=49 MAD(보정 전) = d_(25)            = 2,486,373

import pandas as pd
import numpy as np
from statsmodels import robust

# 먼저 원래 자료에서 두 척도를 재 둔다.
url = ('https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/state.csv')
state = pd.read_csv(url)
original_std = state['Population'].std()
original_mad = robust.scale.mad(state['Population'])

# 여기에 있을 수 없을 만큼 큰 가상의 주 둘을 끼워 넣는다.
population_with_outliers = pd.concat([
    state['Population'],
    pd.Series([100_000_000, 150_000_000])
])

outlier_std = population_with_outliers.std()
outlier_mad = robust.scale.mad(population_with_outliers)

# 자료 52개 중 둘만 바뀌었는데 두 척도가 받는 충격은 전혀 다르다.
print("이상치의 영향:")
print(f"  표준편차: {original_std:,.0f} → {outlier_std:,.0f} ({100 * (outlier_std - original_std) / original_std:.1f}% 증가)")
print(f"  MAD     : {original_mad:,.0f} → {outlier_mad:,.0f} ({100 * (outlier_mad - original_mad) / original_mad:.1f}% 증가)")

# (1) MAD 의 11% 가 이상치의 크기와 무관함을 보인다. 크기만 바꿔 가며 다시 잰다.
print("\n이상치의 크기를 키워 가며:")
print(f"{'이상치 둘':>22}{'s':>14}{'MAD(보정)':>16}")
for a, b in [(1e8, 1.5e8), (1e9, 1e10), (1e12, 1e13), (1e30, 1e31)]:
    q = pd.concat([state['Population'].astype(float), pd.Series([a, b])])
    print(f"{a:>10.1e},{b:>10.1e}{q.std():>14.4e}{robust.scale.mad(q):>16,.0f}")

# 왜 그런가: MAD 는 순서통계량만 본다.
pw = population_with_outliers.astype(float).values
srt = np.sort(pw)
med = np.median(pw)
d = np.sort(np.abs(pw - med))
print(f"\n  n=52 중앙값 = (x_(26)+x_(27))/2 = ({srt[25]:,.0f}+{srt[26]:,.0f})/2 = {med:,.1f}")
print(f"  n=52 MAD(보정 전) = (d_(26)+d_(27))/2 = ({d[25]:,.0f}+{d[26]:,.0f})/2 = {(d[25] + d[26]) / 2:,.1f}")
print(f"  두 이상치의 편차는 순위 51, 52 번: {d[-2]:,.0f}, {d[-1]:,.0f}  (중앙값 계산에 쓰이지 않는다)")

# s 쪽은 두 이상치가 제곱합을 거의 독점한다.
pwm = pw.mean()
SS = ((pw - pwm) ** 2).sum()
share = ((1e8 - pwm) ** 2 + (1.5e8 - pwm) ** 2) / SS
print(f"\n  새 평균 {pwm:,.0f},  두 이상치가 제곱합에서 차지하는 몫 {share:.4f}")
print(f"  두 이상치만으로 어림한 s = {np.sqrt(((1e8 - pwm) ** 2 + (1.5e8 - pwm) ** 2) / 51):,.0f}"
      f"  (실제 {outlier_std:,.0f},  비 {np.sqrt(share):.4f})")

이상치의 영향:
  표준편차: 6,848,235 → 24,537,372 (258.3% 증가)
  MAD     : 3,849,876 → 4,273,462 (11.0% 증가)

이상치의 크기를 키워 가며:
                 이상치 둘             s         MAD(보정)
   1.0e+08,   1.5e+08    2.4537e+07       4,273,462
   1.0e+09,   1.0e+10    1.3901e+09       4,273,462
   1.0e+12,   1.0e+13    1.3910e+12       4,273,462
   1.0e+30,   1.0e+31    1.3910e+30       4,273,462

  n=52 중앙값 = (x_(26)+x_(27))/2 = (4,533,372+4,625,364)/2 = 4,579,368.0
  n=52 MAD(보정 전) = (d_(26)+d_(27))/2 = (2,753,027+3,011,786)/2 = 2,882,406.5
  두 이상치의 편차는 순위 51, 52 번: 95,420,632, 145,420,632  (중앙값 계산에 쓰이지 않는다)

  새 평균 10,733,535,  두 이상치가 제곱합에서 차지하는 몫 0.8911
  두 이상치만으로 어림한 s = 23,163,380  (실제 24,537,372,  비 0.9440)

import numpy as np
import pandas as pd
from statsmodels import robust

# 어느 주식의 일별 수익률이라고 하자. 마지막 하루가 폭락일(-50%)이다.
returns = pd.Series([0.01, 0.02, -0.01, 0.015, -0.005, 0.03, -0.02,
                      0.01, -0.01, 0.005, -0.015, 0.02, -0.50])

print(f"표준편차   : {returns.std():.4f}")
print(f"MAD (보정) : {robust.scale.mad(returns):.4f}")
print(f"s / MAD    : {returns.std() / robust.scale.mad(returns):.4f}")

# 손계산을 따라가 본다. n = 13 이라 중앙값은 7번째 값 하나다.
n = len(returns)
srt = np.sort(returns.values)
med = srt[n // 2]
d = np.sort(np.abs(returns.values - med))
print(f"\nn = {n},  정렬한 수익률 {srt}")
print(f"중앙값 = x_(7) = {med}")
print(f"절대편차 정렬 {d}")
print(f"MAD(보정 전) = d_(7) = {d[n // 2]},  보정 후 {d[n // 2] / 0.6744897501960817:.6f}")

# (1) 한 점이 압도할 때 s = |x_0|/sqrt(n) 인가. 폭락 깊이를 키워 가며 본다.
print(f"\n{'폭락일':>10}{'s 실제':>12}{'|x0|/sqrt(n)':>15}{'비':>9}{'MAD(보정)':>12}")
for crash in (-0.50, -5.00, -50.00):
    q = pd.Series(list(returns[:-1]) + [crash])
    approx = abs(crash) / np.sqrt(n)
    print(f"{crash:>10.2f}{q.std():>12.4f}{approx:>15.4f}"
          f"{q.std() / approx:>9.4f}{robust.scale.mad(q):>12.4f}")

# 폭락일을 아예 빼면?
q0 = returns[:-1]
print(f"\n폭락일 제외 (n={len(q0)}):  s = {q0.std():.4f},  MAD(보정) = {robust.scale.mad(q0):.4f}")
print(f"  폭락일이 s 를 {returns.std() / q0.std():.2f}배, MAD 를 "
      f"{robust.scale.mad(returns) / robust.scale.mad(q0):.2f}배로 만들었다")

표준편차   : 0.1407
MAD (보정) : 0.0222
s / MAD    : 6.3249

n = 13,  정렬한 수익률 [-0.5   -0.02  -0.015 -0.01  -0.01  -0.005  0.005  0.01   0.01   0.015
  0.02   0.02   0.03 ]
중앙값 = x_(7) = 0.005
절대편차 정렬 [0.    0.005 0.005 0.01  0.01  0.015 0.015 0.015 0.015 0.02  0.025 0.025
 0.505]
MAD(보정 전) = d_(7) = 0.015,  보정 후 0.022239

       폭락일        s 실제   |x0|/sqrt(n)        비     MAD(보정)
     -0.50      0.1407         0.1387   1.0143      0.0222
     -5.00      1.3880         1.3868   1.0009      0.0222
    -50.00     13.8687        13.8675   1.0001      0.0222

폭락일 제외 (n=12):  s = 0.0159,  MAD(보정) = 0.0185
  폭락일이 s 를 8.83배, MAD 를 1.20배로 만들었다

from statsmodels import robust
import numpy as np
import pandas as pd

# 마지막 100 이 이상치다. 나머지 다섯 값은 1부터 5까지 고르게 놓여 있다.
data = pd.Series([1, 2, 3, 4, 5, 100])
mad = robust.scale.mad(data)
print(f"MAD: {mad:.2f}")

# 손계산을 따라간다. n = 6 이라 중앙값은 3번째와 4번째의 평균이다.
srt = np.sort(data.values.astype(float))
med = (srt[2] + srt[3]) / 2
d = np.sort(np.abs(srt - med))
print(f"\n정렬 {srt},  중앙값 = ({srt[2]}+{srt[3]})/2 = {med}")
print(f"절대편차 정렬 {d}")
print(f"MAD(보정 전) = (d_(3)+d_(4))/2 = ({d[2]}+{d[3]})/2 = {(d[2] + d[3]) / 2}")
print(f"보정 후 = {(d[2] + d[3]) / 2 / 0.6744897501960817:.6f}")

# 이상치만 키워 가며 두 척도를 본다. s 는 보기 3 의 |x_0|/sqrt(n) 을 따라야 한다.
n = len(data)
print(f"\n{'최댓값':>10}{'MAD(보정)':>12}{'s 실제':>14}{'|x0|/sqrt(6)':>15}{'비':>9}")
for big in (100.0, 1e4, 1e6, 1e12):
    dd = pd.Series([1, 2, 3, 4, 5, big])
    approx = big / np.sqrt(n)
    print(f"{big:>10.0e}{robust.scale.mad(dd):>12.4f}{dd.std():>14.4e}"
          f"{approx:>15.4e}{dd.std() / approx:>9.4f}")

MAD: 2.22

정렬 [  1.   2.   3.   4.   5. 100.],  중앙값 = (3.0+4.0)/2 = 3.5
절대편차 정렬 [ 0.5  0.5  1.5  1.5  2.5 96.5]
MAD(보정 전) = (d_(3)+d_(4))/2 = (1.5+1.5)/2 = 1.5
보정 후 = 2.223903

       최댓값     MAD(보정)          s 실제   |x0|/sqrt(6)        비
     1e+02      2.2239    3.9625e+01     4.0825e+01   0.9706
     1e+04      2.2239    4.0813e+03     4.0825e+03   0.9997
     1e+06      2.2239    4.0825e+05     4.0825e+05   1.0000
     1e+12      2.2239    4.0825e+11     4.0825e+11   1.0000

import numpy as np
import pandas as pd
from statsmodels import robust

data = pd.Series([1, 2, 3, 4, 5, 100])

# 정의를 세 줄로 그대로 옮긴 것이다: 중앙값 → 절대편차 → 그 중앙값.
median = data.median()
abs_dev = abs(data - median)
mad = abs_dev.median()

# 정규분포에서 표준편차와 눈금을 맞추기 위한 보정상수다.
mad_standardized = mad / 0.6744897501960817
print(f"MAD (보정): {mad_standardized:.2f}")
print(f"statsmodels: {robust.scale.mad(data):.2f}  "
      f"(차 {mad_standardized - robust.scale.mad(data):.1e})")

# 두 답이 같은 것은 세 가지가 모두 일치했기 때문이다.
# 그 가운데 분위수 보간법만 바꾸면 답이 달라진다.
x = data.values.astype(float)
print(f"\n{'method':>10}{'중앙값':>9}{'MAD(보정 전)':>14}{'MAD(보정)':>12}")
for meth in ("linear", "lower", "higher", "midpoint", "nearest"):
    m = np.percentile(x, 50, method=meth)
    v = np.percentile(np.abs(x - m), 50, method=meth)
    print(f"{meth:>10}{m:>9.2f}{v:>14.2f}{v / 0.6744897501960817:>12.4f}")

MAD (보정): 2.22
statsmodels: 2.22  (차 0.0e+00)

    method      중앙값     MAD(보정 전)     MAD(보정)
    linear     3.50          1.50      2.2239
     lower     3.00          1.00      1.4826
    higher     4.00          2.00      2.9652
  midpoint     3.50          1.50      2.2239
   nearest     3.00          1.00      1.4826

import numpy as np

rng = np.random.default_rng(0)

def Qn(x):
    x = np.sort(x); k = len(x)
    v = np.abs(x[:, None] - x[None, :])[np.triu_indices(k, 1)]
    h = k // 2 + 1
    return 2.2219 * np.sort(v)[h * (h - 1) // 2 - 1]

def Sn(x):
    k = len(x)
    d = np.abs(x[:, None] - x[None, :])
    return 1.1926 * np.median([np.median(np.delete(d[i], i)) for i in range(k)])

n, B = 40, 8000
X = rng.normal(0, 1, (B, n))
estimators = {
    "s": X.std(axis=1, ddof=1),
    "MAD x 1.4826": 1.4826 * np.median(
        np.abs(X - np.median(X, axis=1, keepdims=True)), axis=1),
    "Qn": np.array([Qn(r) for r in X]),
    "Sn": np.array([Sn(r) for r in X]),
}

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

정규 N(0,1), n=40
추정량                  평균        MSE       상대효율
s                0.9954    0.01246     1.0000
MAD x 1.4826     0.9809    0.03204     0.3890
Qn               1.0963    0.02971     0.4194
Sn               1.0135    0.02187     0.5699

import numpy as np

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

print(f"{'분포':>10}{'MAD x 1.4826':>15}{'SD':>10}{'아래쪽 편차':>14}{'위쪽 편차':>13}{'비':>8}")
for label, x in [("정규", rng.normal(0, 1, n)),
                 ("지수", rng.exponential(1, n)),
                 ("로그정규", rng.lognormal(0, 1, n))]:
    m = np.median(x)
    mad = np.median(np.abs(x - m))
    lower = np.median(m - x[x < m])       # 아래쪽 편차의 중앙값
    upper = np.median(x[x > m] - m)       # 위쪽 편차의 중앙값
    print(f"{label:>10}{1.4826 * mad:>15.4f}{x.std():>10.4f}"
          f"{lower:>14.4f}{upper:>13.4f}{upper / lower:>8.3f}")

        분포   MAD x 1.4826        SD        아래쪽 편차        위쪽 편차       비
        정규         1.0023    1.0012        0.6741       0.6776   1.005
        지수         0.7134    0.9997        0.4064       0.6928   1.704
      로그정규         0.8852    2.1585        0.4896       0.9543   1.949

import numpy as np

def Qn(x):
    x = np.sort(np.asarray(x, float)); k = len(x)
    v = np.abs(x[:, None] - x[None, :])[np.triu_indices(k, 1)]
    h = k // 2 + 1
    return 2.2219 * np.sort(v)[h * (h - 1) // 2 - 1]

d = np.array([5, 5, 5, 5, 5, 5, 5, 6, 7, 40], float)
med = np.median(d)
mad = np.median(np.abs(d - med))

print(f"자료 {d.astype(int)}")
print(f"  중앙값 {med}   MAD {mad}")
print(f"  → 수정 Z 점수 = 0.6745(x - 5)/0  : 0 으로 나누게 되어 계산 불가")
print(f"\n  IQR {np.subtract(*np.percentile(d, [75, 25])):.4f}")
print(f"  Qn  {Qn(d):.4f}")
print(f"  표본표준편차 {d.std(ddof=1):.4f}")

자료 [ 5  5  5  5  5  5  5  6  7 40]
  중앙값 5.0   MAD 0.0
  → 수정 Z 점수 = 0.6745(x - 5)/0  : 0 으로 나누게 되어 계산 불가

  IQR 0.7500
  Qn  0.0000
  표본표준편차 10.9828

import numpy as np

rng = np.random.default_rng(11)

def summarize(x, name):
    med = np.median(x)
    mad = 1.4826 * np.median(np.abs(x - med))
    q1, q3 = np.percentile(x, [25, 75])
    print(f"{name:<22}{x.std(ddof=1):>10.3f}{mad:>10.3f}{(q3 - q1) / 1.349:>10.3f}"
          f"{np.mean(np.abs(x - med) > 3 * mad) * 100:>12.2f}%")

print(f"{'자료':<22}{'s':>10}{'MAD*c':>10}{'IQR/1.349':>10}{'|MZ|>3 비율':>13}")
summarize(rng.normal(100, 15, 5000), "정규")
summarize(np.r_[rng.normal(100, 15, 4750), rng.normal(100, 90, 250)], "5% 오염")
summarize(rng.lognormal(4.6, 0.8, 5000), "로그정규(치우침)")
summarize(np.r_[np.full(3500, 5.0), rng.integers(6, 12, 1500)], "70% 동점")

자료                             s     MAD*c IQR/1.349    |MZ|>3 비율
정규                        14.992    14.782    14.697        0.30%
5% 오염                     24.122    16.049    16.034        3.04%
로그정규(치우침)                125.395    73.574    83.564        7.30%
70% 동점                     1.840     0.000     0.741       30.00%