경험적 분포와 분위수-분위수 그림¶

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

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

np.random.seed(1)
x = 4 + np.random.normal(0, 1.5, 100)     # 평균 4, 표준편차 1.5의 정규 표본 100개

# 표본에서 모수를 추정한다. 참값(4, 1.5)이 아니라 자료에서 잰 값을 쓴다.
loc = x.mean()
scale = x.std()

# 이론적 CDF는 x가 정렬되어 있어야 선으로 이어 그릴 수 있다
x.sort()
cdf = stats.norm(loc=loc, scale=scale).cdf(x)

fig, ax = plt.subplots(figsize=(12, 3))

# ax.ecdf 가 경험적 누적분포함수를 그린다.
# 자료점마다 1/n 씩 올라가는 계단함수이며, 여기서는 n=100 이라 계단이 촘촘해
# 매끄러운 곡선처럼 보인다.
ax.ecdf(x, ls="-", c="#D32F2F", label="경험적 누적분포함수")

# 같은 자료에 적합한 정규분포의 이론적 CDF를 겹친다.
# 두 곡선의 벌어짐이 곧 "정규분포 가정이 얼마나 맞는가"이다.
ax.plot(x, cdf, "-", c="#1565C0", label="적합된 정규분포의 누적분포함수")
ax.set_xlabel("관측값")
ax.set_ylabel("누적 비율")
ax.set_title("경험적 누적분포함수와 이론적 누적분포함수")
ax.legend()
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("ecdf_25.png", dpi=170, facecolor="white", bbox_inches="tight")

# --- (1) 뜀의 크기가 정말 1/n 인가 ---
n = len(x)
print(f"n = {n},  서로 다른 값의 개수 = {len(np.unique(x))}")
print(f"뜀의 크기로 나타난 값 = {np.unique(np.round(np.diff(np.arange(1, n + 1) / n), 12))}")

# --- (1) D 를 순서통계량만으로 계산한다 ---
# [x_(i), x_(i+1)) 에서 ECDF 는 i/n 으로 납작하고 F 는 올라가므로
# 위로 벌어지는 최대는 왼쪽 끝, 아래로 벌어지는 최대는 왼쪽 끝 바로 앞에서 난다.
i = np.arange(1, n + 1)
Dplus = (i / n - cdf).max()
Dminus = (cdf - (i - 1) / n).max()
print(f"D+ = {Dplus:.6f}  (i = {(i / n - cdf).argmax() + 1}번째 순서통계량)")
print(f"D- = {Dminus:.6f}  (i = {(cdf - (i - 1) / n).argmax() + 1}번째 순서통계량)")

# 두 곡선의 최대 수직거리가 콜모고로프-스미르노프 통계량 D 다.
# 이 눈대중을 형식적 검정으로 만든 것이 KS 검정이다.
D, pval = stats.kstest(x, stats.norm(loc=loc, scale=scale).cdf)
print(f"최대 수직거리 D = {D:.4f}   (손으로 구한 max(D+, D-) = {max(Dplus, Dminus):.4f})")
print(f"KS 검정 p값     = {pval:.4f}")

# --- (2) 두 경로가 주는 신뢰띠 폭 ---
alpha = 0.05
eps_dkw = np.sqrt(np.log(2 / alpha) / (2 * n))          # DKW 부등식
t95 = stats.kstwobign.ppf(0.95)                          # 콜모고로프 극한분포의 95분위
print(f"\nDKW 띠        eps = {eps_dkw:.7f}")
print(f"콜모고로프 띠 eps = {t95 / np.sqrt(n):.7f}   (t_0.95 = {t95:.6f})")
# DKW 의 오른변은 콜모고로프 꼬리급수의 첫 항이다. 둘째 항의 크기를 재 보면
# 왜 두 수가 소수 여섯째 자리까지 같은지 알 수 있다.
print(f"  첫 항 2exp(-2t^2)   = {2 * np.exp(-2 * t95 ** 2):.8f}")
print(f"  둘째 항 -2exp(-8t^2) = {-2 * np.exp(-8 * t95 ** 2):.2e}")
print(f"  콜모고로프 꼬리 정확값 = {stats.kstwobign.sf(t95):.8f}")

# 유한한 n 에서의 정확한 포함률은 kstwo 가 준다 (극한이 아니라 n=100 의 분포다).
print(f"n = {n} 에서 정확한 포함률 = {stats.kstwo.cdf(eps_dkw, n):.6f}")

# --- (2) 모의실험으로 포함률을 재고, 분포에 무관함도 함께 본다 ---
rng = np.random.default_rng(0)
reps = 2000
for name, draw, F in (
        ("정규", lambda m: rng.normal(4, 1.5, size=(m, n)), stats.norm(4, 1.5).cdf),
        ("지수", lambda m: rng.exponential(1.0, size=(m, n)), stats.expon(scale=1.0).cdf),
):
    s = np.sort(draw(reps), axis=1)
    Fv = F(s)
    sup = np.maximum((i / n - Fv).max(axis=1), (Fv - (i - 1) / n).max(axis=1))
    cov = (sup <= eps_dkw).mean()
    print(f"{name}모집단 포함률 = {cov:.4f} +- {np.sqrt(cov * (1 - cov) / reps):.4f}  "
          f"(평균 sup = {sup.mean():.4f})")

# --- 신뢰띠 그림 ---
ecdf_y = i / n
fig, ax = plt.subplots(figsize=(12, 3.2))
ax.step(x, ecdf_y, where="post", c="#D32F2F", label="경험적 누적분포함수")
ax.fill_between(x, np.clip(ecdf_y - eps_dkw, 0, 1), np.clip(ecdf_y + eps_dkw, 0, 1),
                step="post", color="#DCEBFB", label="DKW 신뢰띠 (95%)")
ax.plot(x, cdf, "-", c="#1565C0", label="적합된 정규분포의 누적분포함수")
ax.set_xlabel("관측값")
ax.set_ylabel("누적 비율")
ax.set_title("경험적 누적분포함수의 DKW 신뢰띠")
ax.legend(loc="upper left")
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("ecdf_dkw.png", dpi=170, facecolor="white", bbox_inches="tight")

n = 100,  서로 다른 값의 개수 = 100
뜀의 크기로 나타난 값 = [0.01]
D+ = 0.043818  (i = 44번째 순서통계량)
D- = 0.037042  (i = 6번째 순서통계량)
최대 수직거리 D = 0.0438   (손으로 구한 max(D+, D-) = 0.0438)
KS 검정 p값     = 0.9863

DKW 띠        eps = 0.1358102
콜모고로프 띠 eps = 0.1358099   (t_0.95 = 1.358099)
  첫 항 2exp(-2t^2)   = 0.05000078
  둘째 항 -2exp(-8t^2) = -7.81e-07
  콜모고로프 꼬리 정확값 = 0.05000000
n = 100 에서 정확한 포함률 = 0.954666
정규모집단 포함률 = 0.9480 +- 0.0050  (평균 sup = 0.0854)
지수모집단 포함률 = 0.9560 +- 0.0046  (평균 sup = 0.0850)

import matplotlib
matplotlib.use("Agg")
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

loc = 1        # 평균
scale = 2      # 표준편차
normal = stats.norm(loc=loc, scale=scale)

# 평균에서 좌우 3 표준편차까지를 촘촘히 훑는다.
# 정규분포는 이 범위에 확률의 99.7%가 들어 있다.
x = np.linspace(loc - 3 * scale, loc + 3 * scale, 1_000)
pdf = normal.pdf(x)     # 밀도함수: 각 점에서의 "빽빽함"
cdf = normal.cdf(x)     # 분포함수: 그 점까지 누적된 확률

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, pdf, "-", c="#1565C0", label="확률밀도함수")
ax.plot(x, cdf, "-", c="#D32F2F", label="누적분포함수")
ax.set_xlabel("$x$")
ax.set_ylabel("밀도 / 누적 확률")
ax.set_title("정규분포의 확률밀도함수와 누적분포함수")
ax.legend()
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("ecdf_50.png", dpi=170, facecolor="white", bbox_inches="tight")

# 두 함수의 관계를 숫자로 확인한다.
#   CDF는 PDF를 적분한 것이므로 평균에서 정확히 0.5,
#   PDF가 최대인 곳(평균)에서 CDF의 기울기가 가장 가파르다.
print(f"CDF(평균)     = {normal.cdf(loc):.4f}")
print(f"PDF 최댓값    = {pdf.max():.4f}  (x = {x[pdf.argmax()]:.2f})")
print(f"P(|X-mu|<3s)  = {normal.cdf(loc+3*scale) - normal.cdf(loc-3*scale):.4f}")

# 격자가 평균을 비켜 간다. 봉우리의 "자리"는 격자만큼만 정확하다.
print(f"  이론 최댓값 1/(s*sqrt(2pi)) = {1 / (scale * np.sqrt(2 * np.pi)):.6f}")
print(f"  격자 간격 = {x[1] - x[0]:.6f},  평균 1 에 가장 가까운 격자점 = {x[np.abs(x - loc).argmin()]:.6f}")

# --- 같은 모집단에서 n=100 을 뽑아 ECDF 의 한 점을 재 본다 ---
# 점 t 를 고정하면 n*Fhat(t) ~ Bin(n, F(t)) 이므로
#   E[Fhat(t)] = F(t),  Var[Fhat(t)] = F(t)(1-F(t))/n
# 이고 분산은 F = 1/2, 곧 t = 평균에서 가장 크다.
n, reps = 100, 20_000
rng = np.random.default_rng(2)
sample = rng.normal(loc, scale, size=(reps, n))
print("\n  t     F(t)     E[Fhat] 실측   sd 이론    sd 실측")
for t in (-1.0, 1.0, 3.0):
    Ft = normal.cdf(t)
    hat = (sample <= t).mean(axis=1)
    print(f"{t:5.1f}  {Ft:.6f}  {hat.mean():.6f}    "
          f"{np.sqrt(Ft * (1 - Ft) / n):.6f}  {hat.std(ddof=1):.6f}")

CDF(평균)     = 0.5000
PDF 최댓값    = 0.1995  (x = 0.99)
P(|X-mu|<3s)  = 0.9973
  이론 최댓값 1/(s*sqrt(2pi)) = 0.199471
  격자 간격 = 0.012012,  평균 1 에 가장 가까운 격자점 = 0.993994

  t     F(t)     E[Fhat] 실측   sd 이론    sd 실측
 -1.0  0.158655  0.158431    0.036535  0.036318
  1.0  0.500000  0.499498    0.050000  0.050071
  3.0  0.841345  0.841218    0.036535  0.036537

import pandas as pd
import numpy as np
from scipy import stats

data = {'x': [4, 4, 6, 7, 10, 11, 12, 14, 15]}
df = pd.DataFrame(data)

# 같은 75번째 백분위수를 세 라이브러리로 구한다. 인자의 단위가 서로 다르다.
# 기본 보간법이 셋 다 선형이라 값은 일치한다.
print(f"{df.x.quantile(0.75) = }")                   # pandas: 비율 [0, 1]
print(f"{np.percentile(df.x.values, 75) = }")        # numpy: 백분율 [0, 100]
print(f"{stats.scoreatpercentile(df.x.values, 75) = }")   # scipy: 백분율 [0, 100]

# 셋이 같은 이유는 모두 같은 보간법("linear", Hyndman-Fan 7형)을 쓰기 때문이다.
# 가상 지표 h = (n-1)p + 1 이 정수로 떨어지면 보간이 일어나지 않는다.
s = np.sort(df.x.values)
n = len(s)
for p in (0.25, 0.50, 0.75):
    h = (n - 1) * p + 1
    print(f"p = {p:.2f}:  h = (n-1)p+1 = {h:.1f}  ->  x_({int(h)}) = {s[int(h) - 1]}")

# --- 보간법을 바꾸면 사분위수가 달라진다 ---
# numpy 는 Hyndman-Fan 의 아홉 가지를 모두 제공한다. 어느 것도 틀리지 않았다.
methods = ["inverted_cdf", "averaged_inverted_cdf", "closest_observation",
           "interpolated_inverted_cdf", "hazen", "weibull", "linear",
           "median_unbiased", "normal_unbiased"]
print(f"\n{'method':26s} {'Q1':>6s} {'Med':>6s} {'Q3':>6s} {'IQR':>6s}   위쪽 울타리")
for m in methods:
    q1, q2, q3 = (np.percentile(s, p, method=m) for p in (25, 50, 75))
    iqr = q3 - q1
    print(f"{m:26s} {q1:6.3f} {q2:6.3f} {q3:6.3f} {iqr:6.3f}   {q3 + 1.5 * iqr:7.3f}")

df.x.quantile(0.75) = 12.0
np.percentile(df.x.values, 75) = 12.0
stats.scoreatpercentile(df.x.values, 75) = 12.0
p = 0.25:  h = (n-1)p+1 = 3.0  ->  x_(3) = 6
p = 0.50:  h = (n-1)p+1 = 5.0  ->  x_(5) = 10
p = 0.75:  h = (n-1)p+1 = 7.0  ->  x_(7) = 12

method                         Q1    Med     Q3    IQR   위쪽 울타리
inverted_cdf                6.000 10.000 12.000  6.000    21.000
averaged_inverted_cdf       6.000 10.000 12.000  6.000    21.000
closest_observation         4.000 10.000 12.000  8.000    24.000
interpolated_inverted_cdf   4.500  8.500 11.750  7.250    22.625
hazen                       5.500 10.000 12.500  7.000    23.000
weibull                     5.000 10.000 13.000  8.000    25.000
linear                      6.000 10.000 12.000  6.000    21.000
median_unbiased             5.333 10.000 12.667  7.333    23.667
normal_unbiased             5.375 10.000 12.625  7.250    23.500

import matplotlib
matplotlib.use("Agg")
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

# 당 함량을 5g 간격으로 끊고, 각 지점까지 누적된 비율을 기록한 자료다.
# y가 단조 증가하고 마지막이 1.0으로 끝나는 것이 누적상대도수의 성질이다.
x = np.arange(0, 55, 5)
y = [0, 0.1, 0.1, 0.2, 0.3, 0.5, 0.6, 0.6, 0.8, 0.9, 1.0]

fig, ax = plt.subplots(figsize=(12, 3))
ax.plot(x, y, '-o', c="#1565C0")             # 점을 찍고 이어 그린다
ax.set_xlabel("당 함량 (g)")
ax.set_ylabel("누적 상대도수")
ax.set_title("스타벅스 음료 32종의 당 함량 누적상대도수 곡선")
# y 눈금을 0.1 간격으로 촘촘히 두어야 백분위수를 눈으로 읽을 수 있다
ax.set_yticks(np.arange(0, 1.1, 0.1))
ax.grid()                                    # 격자가 있어야 가로세로로 읽어 나가기 쉽다
fig.savefig("ecdf_138.png", dpi=170, facecolor="white", bbox_inches="tight")

# 그림에서 눈으로 읽는 값을 코드로도 구해 본다.
# 누적비율 y에서 가로로 이동해 곡선을 만나는 x가 그 백분위수다.
for p_ in (0.25, 0.50, 0.75):
    print(f"P{int(p_*100)} = {np.interp(p_, y, x):.1f} g")

# 손으로 푼 식과 맞는지 확인한다. p 가 든 구간의 양 끝 (x_k, y_k), (x_{k+1}, y_{k+1}) 에서
#   P_p = x_k + (x_{k+1} - x_k) * (p - y_k) / (y_{k+1} - y_k)
y = np.asarray(y, float)
for p_ in (0.25, 0.50, 0.75):
    k = int(np.searchsorted(y, p_, side="left")) - 1
    hand = x[k] + (x[k + 1] - x[k]) * (p_ - y[k]) / (y[k + 1] - y[k])
    print(f"  손으로: x_{k} = {x[k]}, y_{k} = {y[k]:.1f} -> P{int(p_*100)} = {hand:.2f} g")
print(f"IQR = {np.interp(0.75, y, x) - np.interp(0.25, y, x):.2f} g")
print(f"당 15 g 인 음료의 누적비율 = {np.interp(15, x, y):.2f}")

# 평평한 구간에서는 역이 하나로 정해지지 않는다.
# y 가 0.1 인 구간이 [5, 10], 0.6 인 구간이 [30, 35] 다.
for p_ in (0.1, 0.6):
    lo = x[np.flatnonzero(y == p_)[0]]
    hi = x[np.flatnonzero(y == p_)[-1]]
    print(f"F(x) = {p_} 인 x 의 범위 = [{lo}, {hi}],  "
          f"np.interp 의 답 = {np.interp(p_, y, x):.1f},  "
          f"inf 규약의 답 = {lo}")

P25 = 17.5 g
P50 = 25.0 g
P75 = 38.8 g
  손으로: x_3 = 15, y_3 = 0.2 -> P25 = 17.50 g
  손으로: x_4 = 20, y_4 = 0.3 -> P50 = 25.00 g
  손으로: x_7 = 35, y_7 = 0.6 -> P75 = 38.75 g
IQR = 21.25 g
당 15 g 인 음료의 누적비율 = 0.20
F(x) = 0.1 인 x 의 범위 = [5, 10],  np.interp 의 답 = 10.0,  inf 규약의 답 = 5
F(x) = 0.6 인 x 의 범위 = [30, 35],  np.interp 의 답 = 35.0,  inf 규약의 답 = 30

import matplotlib
matplotlib.use("Agg")
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

data = np.array([1, 2, 0, 0, 0, 1, 3, 1, 2, 1, 2, 4, 5, -1, -2, 0, 8])

# 다섯 수치 요약은 최소·Q1·중앙값·Q3·최대다. 모두 분위수이므로 q 만 바꿔 부른다.
quantiles = {"Min": 0, "Q1": 0.25, "Median": 0.5, "Q3": 0.75, "Max": 1}

for label, q in quantiles.items():
    print(f"{label:6} : {np.quantile(data, q)}")

# 상자그림은 이 다섯 수를 그림으로 옮긴 것이다. 상자의 위아래가 Q3와 Q1,
# 가운데 선이 중앙값이며, 수염 밖에 찍히는 점이 이상치 후보다.
fig, ax = plt.subplots(figsize=(3, 4))
ax.boxplot(data)
ax.set_xticks([1])
ax.set_xticklabels(["자료"])
ax.set_ylabel("관측값")
ax.set_title("다섯 수치 요약의 상자그림")
fig.savefig("ecdf_216.png", dpi=170, facecolor="white", bbox_inches="tight")

# --- 세 사분위수가 보간 없이 순서통계량에 그대로 떨어지는지 확인한다 ---
s = np.sort(data)
n = len(s)
print(f"\n정렬: {list(s)}")
for p in (0.25, 0.50, 0.75):
    h = (n - 1) * p + 1
    print(f"  p = {p:.2f}: h = (n-1)p+1 = {h:.1f} -> x_({int(h)}) = {s[int(h) - 1]}")

# --- 울타리와 수염 끝 ---
q1, q3 = np.quantile(data, 0.25), np.quantile(data, 0.75)
iqr = q3 - q1
lo, hi = q1 - 1.5 * iqr, q3 + 1.5 * iqr
print(f"\nIQR = {iqr}, 울타리 = [{lo}, {hi}]")
print(f"  아래 수염 끝 = {s[s >= lo].min()},  위 수염 끝 = {s[s <= hi].max()}")
print(f"  이상치 후보 = {list(s[(s < lo) | (s > hi)])}")

# --- 상자그림이 가리는 것 ---
vals, cnt = np.unique(s, return_counts=True)
print(f"\n값과 도수: {dict(zip(vals.tolist(), cnt.tolist()))}")
print(f"  평균 = {data.mean():.4f},  중앙값 = {np.median(data)}")

Min    : -2
Q1     : 0.0
Median : 1.0
Q3     : 2.0
Max    : 8

정렬: [-2, -1, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 3, 4, 5, 8]
  p = 0.25: h = (n-1)p+1 = 5.0 -> x_(5) = 0
  p = 0.50: h = (n-1)p+1 = 9.0 -> x_(9) = 1
  p = 0.75: h = (n-1)p+1 = 13.0 -> x_(13) = 2

IQR = 2.0, 울타리 = [-3.0, 5.0]
  아래 수염 끝 = -2,  위 수염 끝 = 5
  이상치 후보 = [8]

값과 도수: {-2: 1, -1: 1, 0: 4, 1: 4, 2: 3, 3: 1, 4: 1, 5: 1, 8: 1}
  평균 = 1.5882,  중앙값 = 1.0

"""표본의 분위수를 이론 분포의 분위수에 맞대어 그리는 Q-Q 그림을 만든다."""
import matplotlib
matplotlib.use("Agg")
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

plt.rcParams["font.family"] = "Apple SD Gothic Neo"
plt.rcParams["axes.unicode_minus"] = False

def plot_qq(data, dist="norm", sparams=(), title="", fname="qq.png",
            figsize=(12, 3)):
    """자료를 dist 의 분위수에 대해 그린다. 점이 직선 위에 놓이면 그 분포에 맞는다."""
    fig, ax = plt.subplots(figsize=figsize)
    # sparams 는 분포의 모양모수다. 정규처럼 위치·척도만 있는 분포는 비워 두면
    # probplot 이 자료에서 추정한다. 카이제곱처럼 모양모수가 있으면 넘겨야 한다.
    stats.probplot(data, dist=dist, sparams=sparams, plot=ax)
    ax.get_lines()[0].set_color("#1565C0")     # 자료점
    ax.get_lines()[1].set_color("#D32F2F")     # 기준선
    ax.spines[["top", "right"]].set_visible(False)
    ax.set_title(title)
    ax.set_xlabel("이론 분포의 분위수")
    ax.set_ylabel("정렬된 관측값")
    fig.savefig(fname, dpi=170, facecolor="white", bbox_inches="tight")

np.random.seed(0)
sample_data = np.random.normal(loc=0, scale=1, size=1000)
# 정규 자료를 정규에 맞댄다 — 직선이 나온다
plot_qq(sample_data, dist="norm",
        title="정규 자료를 정규 분위수에 맞댄 Q-Q 그림",
        fname="ecdf_262.png")

# --- 가로축에 놓이는 값과 기준선의 정체를 확인한다 ---
(osm, osr), (slope, intercept, r) = stats.probplot(sample_data, dist="norm")
n = len(sample_data)
i = np.arange(1, n + 1)
pp = (i - 0.3175) / (n + 0.365)       # scipy 가 쓰는 자리잡기(plotting position)
pp[0] = 1 - 0.5 ** (1 / n)
pp[-1] = 0.5 ** (1 / n)
print(f"osm 이 Phi^(-1)(p_i) 와 같은가: {np.allclose(osm, stats.norm.ppf(pp))}")
print(f"p_i + p_(n+1-i) = 1 인가: {np.allclose(pp + pp[::-1], 1)}")
print(f"osm 의 평균 = {osm.mean():.3e}  (대칭이므로 0)")
print(f"기울기 {slope:.6f}  <-  sum(m*y)/sum(m^2) = {(osm * osr).sum() / (osm ** 2).sum():.6f}")
print(f"절편   {intercept:.6f}  <-  표본평균        = {osr.mean():.6f}")
print(f"표본표준편차 = {sample_data.std(ddof=1):.6f},  r = {r:.6f},  r^2 = {r * r:.6f}")

osm 이 Phi^(-1)(p_i) 와 같은가: True
p_i + p_(n+1-i) = 1 인가: True
osm 의 평균 = 5.684e-17  (대칭이므로 0)
기울기 0.989271  <-  sum(m*y)/sum(m^2) = 0.989271
절편   -0.045257  <-  표본평균        = -0.045257
표본표준편차 = 0.987527,  r = 0.999482,  r^2 = 0.998965

# 지수분포는 오른쪽으로 심하게 치우쳐 있다. 그래도 지수 분위수에 맞대면 직선이 된다.
# Q-Q 그림이 보는 것은 치우침 자체가 아니라 "가정한 분포와 얼마나 맞는가"이다.
np.random.seed(0)
sample_data = np.random.exponential(scale=1, size=1000)
plot_qq(sample_data, dist="expon",
        title="지수 자료를 지수 분위수에 맞댄 Q-Q 그림",
        fname="ecdf_285.png")

(osm_e, osr_e), (slope_e, inter_e, r_e) = stats.probplot(sample_data, dist="expon")
print(f"osm 이 -ln(1-p_i) 와 같은가: {np.allclose(osm_e, -np.log(1 - pp))}")
print(f"기울기 {slope_e:.6f} (척도 1 을 겨냥),  절편 {inter_e:.6f} (위치 0 을 겨냥)")
print(f"절편 = ybar - b*mbar = {osr_e.mean() - slope_e * osm_e.mean():.6f}  "
      f"(osm 평균 = {osm_e.mean():.6f} 이라 0 이 아니다)")
print(f"r^2 = {r_e * r_e:.6f}   <- 정규/정규의 {r * r:.6f} 보다 낮다")
print(f"표본 왜도 = {stats.skew(sample_data):.4f}  (지수의 이론 왜도 = 2)")

osm 이 -ln(1-p_i) 와 같은가: True
기울기 1.034924 (척도 1 을 겨냥),  절편 -0.029782 (위치 0 을 겨냥)
절편 = ybar - b*mbar = -0.029782  (osm 평균 = 0.998453 이라 0 이 아니다)
r^2 = 0.995440   <- 정규/정규의 0.998965 보다 낮다
표본 왜도 = 2.0526  (지수의 이론 왜도 = 2)

# 카이제곱은 자유도라는 모양모수가 있으므로 sparams=(10,) 으로 알려 주어야 한다.
# 이 값을 틀리게 주면 자료가 맞는 분포에서 왔더라도 직선에서 벗어난다.
np.random.seed(0)
sample_data = np.random.chisquare(df=10, size=1000)
plot_qq(sample_data, dist="chi2", sparams=(10,),
        title="카이제곱 자료를 카이제곱 분위수에 맞댄 Q-Q 그림",
        fname="ecdf_293.png")

sy = sample_data.std(ddof=1)
print(f"표본평균 {sample_data.mean():.4f} (이론 10),  표본표준편차 {sy:.4f} (이론 sqrt(20) = {np.sqrt(20):.4f})")
print(f"{'가정한 df':>10s} {'기울기':>8s} {'절편':>9s} {'ybar-b*mbar':>12s} {'sy/sd(osm)':>11s} {'r^2':>9s}")
for d in (5, 10, 20):
    (m, yv), (b, a, rv) = stats.probplot(sample_data, dist="chi2", sparams=(d,))
    print(f"{d:10d} {b:8.4f} {a:9.4f} {yv.mean() - b * m.mean():12.4f} "
          f"{sy / m.std(ddof=1):11.4f} {rv * rv:9.6f}")

표본평균 9.9224 (이론 10),  표본표준편차 4.2814 (이론 sqrt(20) = 4.4721)
    가정한 df      기울기        절편  ybar-b*mbar  sy/sd(osm)       r^2
         5   1.3537    3.1583       3.1583      1.3600  0.990753
        10   0.9598    0.3281       0.3281      0.9607  0.998152
        20   0.6772   -3.6194      -3.6194      0.6789  0.995051

# 같은 카이제곱 자료를 이번에는 정규 분위수에 맞댄다.
# 오른쪽 끝이 직선 위로 휘어 오르는 것이 "정규보다 오른쪽 꼬리가 두껍다"는 신호다.
np.random.seed(0)
sample_data = np.random.chisquare(df=10, size=1000)
plot_qq(sample_data, dist="norm",
        title="카이제곱 자료를 정규 분위수에 맞댄 Q-Q 그림 (잘못된 가정)",
        fname="ecdf_240.png")

(m9, y9), (b9, a9, r9) = stats.probplot(sample_data, dist="norm")
print(f"정규 기준선: 기울기 {b9:.4f}, 절편 {a9:.4f} (= 표본평균 {sample_data.mean():.4f})")
print(f"  r^2 = {r9 * r9:.6f}   <- 올바른 카이제곱 기준의 0.998152 보다 낮다")
print(f"  왜도 {stats.skew(sample_data):.4f} (이론 sqrt(8/10) = {np.sqrt(0.8):.4f}),  "
      f"초과첨도 {stats.kurtosis(sample_data):.4f} (이론 12/10 = 1.2)")
# 가장 큰 잔차가 어디서 나는지 본다
res = y9 - (a9 + b9 * m9)
print(f"  기준선에서 가장 멀리 벗어난 점: 잔차 {res.max():.3f} (가장 큰 관측값 쪽), "
      f"{res.min():.3f} (왼쪽 꼬리 쪽)")

# --- ECDF 는 구간을 고르지 않으므로 두 표본을 바로 견줄 수 있다 ---
rng = np.random.default_rng(7)
normal_like = rng.normal(sample_data.mean(), sy, 1000)
pool = np.sort(np.concatenate([sample_data, normal_like]))
F1 = np.searchsorted(np.sort(sample_data), pool, side="right") / len(sample_data)
F2 = np.searchsorted(np.sort(normal_like), pool, side="right") / len(normal_like)
res2 = stats.ks_2samp(sample_data, normal_like)
print(f"\n두 ECDF 의 최대 수직거리 = {np.abs(F1 - F2).max():.6f}")
print(f"  ks_2samp 가 준 D = {res2.statistic:.6f},  p = {res2.pvalue:.4f}")

정규 기준선: 기울기 4.2003, 절편 9.9224 (= 표본평균 9.9224)
  r^2 = 0.958102   <- 올바른 카이제곱 기준의 0.998152 보다 낮다
  왜도 0.8519 (이론 sqrt(8/10) = 0.8944),  초과첨도 0.8951 (이론 12/10 = 1.2)
  기준선에서 가장 멀리 벗어난 점: 잔차 5.171 (가장 큰 관측값 쪽), -0.590 (왼쪽 꼬리 쪽)

두 ECDF 의 최대 수직거리 = 0.051000
  ks_2samp 가 준 D = 0.051000,  p = 0.1484

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)
band = lambda n, alpha=0.05: np.sqrt(np.log(2 / alpha) / (2 * n))

print("DKW 신뢰띠: 곡선 전체가 띠 안에 들어갈 확률")
for n in (20, 50, 200, 1000):
    eps = band(n)
    covered = 0
    B = 20_000
    for _ in range(B):
        x = np.sort(rng.normal(0, 1, n))
        F = stats.norm.cdf(x)
        upper = np.arange(1, n + 1) / n          # ECDF 는 계단함수라
        lower = upper - 1 / n                    # 각 점의 양쪽을 모두 본다
        d = max(np.max(np.abs(upper - F)), np.max(np.abs(lower - F)))
        if d <= eps:
            covered += 1
    print(f"  n={n:>5}: 띠 반폭 {eps:.4f}   실제 포함률 {covered / B:.4f}")

DKW 신뢰띠: 곡선 전체가 띠 안에 들어갈 확률
  n=   20: 띠 반폭 0.3037   실제 포함률 0.9604
  n=   50: 띠 반폭 0.1921   실제 포함률 0.9562
  n=  200: 띠 반폭 0.0960   실제 포함률 0.9536
  n= 1000: 띠 반폭 0.0429   실제 포함률 0.9489

import numpy as np

x = np.array([3., 7., 8., 12., 15., 19., 24., 31.])
print(f"자료 {x}  (n={len(x)})\n")
print(f"{'방법':>18}{'Q1':>9}{'중앙값':>9}{'Q3':>9}{'IQR':>9}")
for m in ("linear", "lower", "higher", "nearest", "midpoint",
          "inverted_cdf", "hazen", "weibull", "median_unbiased"):
    q = [np.quantile(x, p, method=m) for p in (0.25, 0.5, 0.75)]
    print(f"{m:>18}{q[0]:>9.3f}{q[1]:>9.3f}{q[2]:>9.3f}{q[2] - q[0]:>9.3f}")

자료 [ 3.  7.  8. 12. 15. 19. 24. 31.]  (n=8)

                방법       Q1      중앙값       Q3      IQR
            linear    7.750   13.500   20.250   12.500
             lower    7.000   12.000   19.000   12.000
            higher    8.000   15.000   24.000   16.000
           nearest    8.000   15.000   19.000   11.000
          midpoint    7.500   13.500   21.500   14.000
      inverted_cdf    7.000   12.000   19.000   12.000
             hazen    7.500   13.500   21.500   14.000
           weibull    7.250   13.500   22.750   15.500
   median_unbiased    7.417   13.500   21.917   14.500

import numpy as np

rng = np.random.default_rng(0)

# 두 집단: 중앙값은 같지만 퍼짐이 다르다
a = rng.normal(100, 5, 300)
b = rng.normal(100, 15, 300)

print(f"{'':8s}{'중앙값':>9s}{'Q1':>8s}{'Q3':>8s}{'IQR':>8s}")
for lab, v in [("A", a), ("B", b)]:
    q1, q2, q3 = np.percentile(v, [25, 50, 75])
    print(f"{lab:8s}{q2:>9.2f}{q1:>8.2f}{q3:>8.2f}{q3 - q1:>8.2f}")

# 두 ECDF의 세로 간격이 가장 벌어지는 지점을 찾는다.
grid = np.linspace(min(a.min(), b.min()), max(a.max(), b.max()), 500)
Fa = np.searchsorted(np.sort(a), grid, side="right") / len(a)
Fb = np.searchsorted(np.sort(b), grid, side="right") / len(b)
gap = np.abs(Fa - Fb)
i = gap.argmax()
print(f"\n두 ECDF의 최대 세로 간격 {gap[i]:.4f}  (x = {grid[i]:.2f})")
print(f"  그 지점에서  A {Fa[i]:.4f}   B {Fb[i]:.4f}")

# 꼬리에서의 차이
for t in [85, 90, 110, 115]:
    print(f"  {t}보다 작은 비율:  A {(a < t).mean():.4f}   B {(b < t).mean():.4f}")

              중앙값      Q1      Q3     IQR
A           99.62   96.47  103.40    6.93
B           99.71   89.81  110.01   20.21

두 ECDF의 최대 세로 간격 0.2500  (x = 92.52)
  그 지점에서  A 0.0433   B 0.2933
  85보다 작은 비율:  A 0.0033   B 0.1533
  90보다 작은 비율:  A 0.0267   B 0.2533
  110보다 작은 비율:  A 0.9800   B 0.7500
  115보다 작은 비율:  A 0.9967   B 0.8333

import numpy as np

rng = np.random.default_rng(3)
A = rng.normal(50, 5, 400)         # 평균 낮고 고르다
B = rng.normal(53, 12, 400)        # 평균 높고 흩어진다

print(f"A: 평균 {A.mean():.2f}  표준편차 {A.std(ddof=1):.2f}")
print(f"B: 평균 {B.mean():.2f}  표준편차 {B.std(ddof=1):.2f}\n")

print(f"{'기준값':>8}{'F_A(x)':>10}{'F_B(x)':>10}{'기준 미만이 많은 쪽':>22}")
for x in (35, 45, 50, 55, 65):
    fa, fb = (A <= x).mean(), (B <= x).mean()
    print(f"{x:>8}{fa:>10.3f}{fb:>10.3f}{('B' if fb > fa else 'A'):>16}")

A: 평균 50.17  표준편차 5.05
B: 평균 53.44  표준편차 11.92

     기준값    F_A(x)    F_B(x)           기준 미만이 많은 쪽
      35     0.003     0.062               B
      45     0.135     0.230               B
      50     0.495     0.362               A
      55     0.840     0.552               A
      65     0.998     0.830               A