고급 기술 측도¶

import numpy as np
from scipy import stats

# 6년치 연간 수익률. 큰 이익과 큰 손실이 섞여 있다.
returns = np.array([0.36, 0.23, -0.48, -0.30, 0.15, 0.31])

# 산술평균: 수익률을 그냥 더해서 나눈다.
arith_mean = np.mean(returns)

# 기하평균: 수익률이 아니라 "성장배수" (1+r)의 기하평균을 낸 뒤 1을 뺀다.
# 복리는 곱셈으로 쌓이므로 곱셈의 평균인 기하평균이 맞는 요약이다.
geo_mean = stats.mstats.gmean(1 + returns) - 1

print(f"Arithmetic mean: {arith_mean:.4f}  ({arith_mean*100:.2f}%)")
print(f"Geometric  mean: {geo_mean:.4f}  ({geo_mean*100:.2f}%)")
# 6년을 실제로 곱해 나갔을 때 1달러가 얼마가 되는지
print(f"Compound value of 1 USD: {np.prod(1 + returns):.4f}")
# 기하평균을 6번 복리로 굴려도 같은 값이 나온다. 이것이 기하평균의 정의다.
print(f"Using geo mean:          {(1 + geo_mean)**len(returns):.4f}")

# 로그를 거친 정의도 같은 값을 준다. 곱을 합으로 바꾼 것이다.
print(f"로그 정의 exp(mean(log(1+r)))-1: {np.exp(np.mean(np.log(1 + returns))) - 1:.4f}")

# (1) 손실 x 를 되돌리는 데 필요한 이익은 x/(1-x) 다.
print()
for x in (0.30, 0.48, 0.50):
    print(f"  손실 {x:.0%} 회복에 필요한 이익 = x/(1-x) = {x / (1 - x):.4f}  ({100 * x / (1 - x):.1f}%)")
print(f"  -48% 와 -30% 를 연달아 맞은 뒤 본전까지: "
      f"{1 / ((1 - 0.48) * (1 - 0.30)) - 1:.4f}  (+174.7%)")

# (2) 산술평균으로 6년 복리를 예측하면 얼마나 틀리는가.
T = len(returns)
print(f"\n  산술평균의 예측 (1+rbar)^T = {(1 + arith_mean) ** T:.4f}")
print(f"  실제            prod(1+r)  = {np.prod(1 + returns):.4f}")
print(f"  비 = {(1 + arith_mean) ** T / np.prod(1 + returns):.4f}  배 과대예측")
print(f"  두 평균의 차 rbar - r_g = {arith_mean - geo_mean:.4f}")

Arithmetic mean: 0.0450  (4.50%)
Geometric  mean: -0.0143  (-1.43%)
Compound value of 1 USD: 0.9173
Using geo mean:          0.9173
로그 정의 exp(mean(log(1+r)))-1: -0.0143

  손실 30% 회복에 필요한 이익 = x/(1-x) = 0.4286  (42.9%)
  손실 48% 회복에 필요한 이익 = x/(1-x) = 0.9231  (92.3%)
  손실 50% 회복에 필요한 이익 = x/(1-x) = 1.0000  (100.0%)
  -48% 와 -30% 를 연달아 맞은 뒤 본전까지: 1.7473  (+174.7%)

  산술평균의 예측 (1+rbar)^T = 1.3023
  실제            prod(1+r)  = 0.9173
  비 = 1.4197  배 과대예측
  두 평균의 차 rbar - r_g = 0.0593

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

# 그림이 가리는 것: 하한과 "실제"의 거리. 세 분포에서 직접 재 본다.
print(f"{'k':>6}{'체비쇼프 하한':>14}{'정규 실제':>11}{'지수 실제':>11}{'균등 실제':>11}")
for k in (1.1, 1.5, 2.0, 3.0, 4.0, 5.0):
    cheb = 1 - 1 / k ** 2
    norm = 2 * stats.norm.cdf(k) - 1
    # 지수(1): mu = 1, sigma = 1
    expo = stats.expon.cdf(1 + k) - stats.expon.cdf(max(0.0, 1 - k))
    # 균등(0,1): mu = 0.5, sigma = 1/sqrt(12)
    s = 1 / np.sqrt(12)
    unif = min(1.0, 0.5 + k * s) - max(0.0, 0.5 - k * s)
    print(f"{k:>6.2f}{cheb:>14.4f}{norm:>11.4f}{expo:>11.4f}{unif:>11.4f}")

print()
for target in (0.75, 0.9, 0.99, 0.999):
    print(f"  하한 {target:.3f} 에 이르는 k = 1/sqrt(1-p) = {1 / np.sqrt(1 - target):.4f}")
print(f"  k=1   에서 하한 = {1 - 1 / 1 ** 2:.4f}  (참이지만 내용이 없다)")
print(f"  k=0.5 에서 하한 = {1 - 1 / 0.5 ** 2:.4f}  (음수라 더욱 공허하다)")

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

def chebyshev(k):
    """평균에서 k 표준편차 안에 들어가는 자료의 최소 비율."""
    return 1 - 1 / k**2

# k=1 에서는 1 - 1 = 0 이라 아무 말도 하지 못한다. 그래서 1.1 부터 그린다.
z_vals = np.arange(1.1, 10, 0.1)
cheb_vals = [chebyshev(z) for z in z_vals]

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(z_vals, cheb_vals, lw=2, color="#33691E")
ax.set_xlabel("$k$ (표준편차의 배수)")
ax.set_ylabel(r"$k\sigma$ 안에 있는 최소 비율")
ax.set_title("체비쇼프 부등식의 하한 $1 - 1/k^2$")
ax.axhline(0.75, color="#90A4AE", linestyle=":", alpha=0.8)
ax.annotate("$k = 2$ 에서 75% 이상", (2, 0.75), fontsize=9,
            xytext=(4, 0.6), arrowprops=dict(arrowstyle="->"))
ax.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
fig.savefig("descriptive_measures_advanced_97.png", dpi=170,
            facecolor="white", bbox_inches="tight")

     k       체비쇼프 하한      정규 실제      지수 실제      균등 실제
  1.10        0.1736     0.7287     0.8775     0.6351
  1.50        0.5556     0.8664     0.9179     0.8660
  2.00        0.7500     0.9545     0.9502     1.0000
  3.00        0.8889     0.9973     0.9817     1.0000
  4.00        0.9375     0.9999     0.9933     1.0000
  5.00        0.9600     1.0000     0.9975     1.0000

  하한 0.750 에 이르는 k = 1/sqrt(1-p) = 2.0000
  하한 0.900 에 이르는 k = 1/sqrt(1-p) = 3.1623
  하한 0.990 에 이르는 k = 1/sqrt(1-p) = 10.0000
  하한 0.999 에 이르는 k = 1/sqrt(1-p) = 31.6228
  k=1   에서 하한 = 0.0000  (참이지만 내용이 없다)
  k=0.5 에서 하한 = -3.0000  (음수라 더욱 공허하다)

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

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

np.random.seed(42)

# 모집단을 직접 만들어 두었으므로 참 분산을 우리는 알고 있다.
population = np.random.normal(170, 10, 1000)
pop_var = np.var(population)

n_samples = 10000
sample_size = 100

biased_vars = np.empty(n_samples)
unbiased_vars = np.empty(n_samples)

# 같은 표본에서 n 으로 나눈 값과 n-1 로 나눈 값을 나란히 기록한다.
for i in range(n_samples):
    sample = np.random.choice(population, size=sample_size, replace=False)
    biased_vars[i] = np.var(sample, ddof=0)
    unbiased_vars[i] = np.var(sample, ddof=1)

# 1만 번의 평균을 참값과 견준다. ddof=0 쪽만 아래로 치우쳐 있다.
print(f"모분산            : {pop_var:.4f}")
print(f"편향 추정 평균 (ddof=0): {biased_vars.mean():.4f}")
print(f"불편 추정 평균 (ddof=1): {unbiased_vars.mean():.4f}")

# --- 과녁이 둘이다. 비복원추출에서 ddof=1 이 맞추는 것은 어느 쪽인가 ---
N = len(population)
S2 = np.var(population, ddof=1)          # 분모 N-1
print(f"\n분모 N   인 모분산 = {pop_var:.4f}")
print(f"분모 N-1 인 S^2    = {S2:.4f}   비 N/(N-1) = {S2 / pop_var:.9f}")
print(f"이론: E[ddof=1] = S^2         = {S2:.4f}")
print(f"이론: E[ddof=0] = (n-1)/n S^2 = {(sample_size - 1) / sample_size * S2:.4f}")

se = unbiased_vars.std(ddof=1) / np.sqrt(n_samples)
print(f"\n모의 평균 - S^2    = {unbiased_vars.mean() - S2:+.4f}  ({(unbiased_vars.mean() - S2) / se:+.2f} SE)")
print(f"모의 평균 - 모분산  = {unbiased_vars.mean() - pop_var:+.4f}  ({(unbiased_vars.mean() - pop_var) / se:+.2f} SE)")
print(f"몬테카를로 표준오차 SE = {se:.4f},  두 과녁의 간격 = {S2 - pop_var:.4f}")
print(f"-> 간격을 2 SE 로 분간하려면 반복 "
      f"{(unbiased_vars.std(ddof=1) / ((S2 - pop_var) / 2)) ** 2:,.0f}회 필요")

# 작은 모집단에서는 모든 표본을 열거해 정확히 확인할 수 있다.
print(f"\n완전열거로 E[s^2] = S^2 확인 (몬테카를로 오차 없음)")
print(f"{'N':>4}{'n':>4}{'표본수':>8}{'E[s^2]':>14}{'S^2':>14}{'E[ddof=0]':>14}{'(n-1)/n S^2':>14}")
for Ns, ns in ((5, 3), (6, 4), (7, 3)):
    rng = np.random.default_rng(7)
    y = rng.normal(0, 1, Ns)
    v1 = [np.var(np.array(c), ddof=1) for c in combinations(y, ns)]
    v0 = [np.var(np.array(c), ddof=0) for c in combinations(y, ns)]
    print(f"{Ns:>4}{ns:>4}{len(v1):>8}{np.mean(v1):>14.10f}{np.var(y, ddof=1):>14.10f}"
          f"{np.mean(v0):>14.10f}{(ns - 1) / ns * np.var(y, ddof=1):>14.10f}")

모분산            : 95.7905
편향 추정 평균 (ddof=0): 94.9057
불편 추정 평균 (ddof=1): 95.8643

분모 N   인 모분산 = 95.7905
분모 N-1 인 S^2    = 95.8864   비 N/(N-1) = 1.001001001
이론: E[ddof=1] = S^2         = 95.8864
이론: E[ddof=0] = (n-1)/n S^2 = 94.9275

모의 평균 - S^2    = -0.0220  (-0.17 SE)
모의 평균 - 모분산  = +0.0738  (+0.56 SE)
몬테카를로 표준오차 SE = 0.1321,  두 과녁의 간격 = 0.0959
-> 간격을 2 SE 로 분간하려면 반복 75,935회 필요

완전열거로 E[s^2] = S^2 확인 (몬테카를로 오차 없음)
   N   n     표본수        E[s^2]           S^2     E[ddof=0]   (n-1)/n S^2
   5   3      10  0.2040262412  0.2040262412  0.1360174941  0.1360174941
   6   4      15  0.2514938167  0.2514938167  0.1886203625  0.1886203625
   7   3      35  0.2379084415  0.2379084415  0.1586056276  0.1586056276

# 그림이 가리는 것을 먼저 수치로 재 둔다.
S2 = np.var(population, ddof=1)
print(f"표본마다 ddof=1 / ddof=0 의 비: 최소 {np.min(unbiased_vars / biased_vars):.10f}, "
      f"최대 {np.max(unbiased_vars / biased_vars):.10f}")
print(f"  n/(n-1) = {sample_size / (sample_size - 1):.10f}")
print(f"\n{'분위수':>8}{'ddof=0':>12}{'ddof=1':>12}{'비':>14}")
for q in (5, 25, 50, 75, 95):
    b, u = np.percentile(biased_vars, q), np.percentile(unbiased_vars, q)
    print(f"{q:>7}%{b:>12.4f}{u:>12.4f}{u / b:>14.10f}")

sd = unbiased_vars.std(ddof=1)
bias = (sample_size - 1) / sample_size * S2 - S2
print(f"\n추정량의 표준편차 {sd:.4f},  ddof=0 의 편향 {bias:.4f}")
print(f"  |편향| / 표준편차 = {abs(bias) / sd:.4f}   <- 그림에서 겹쳐 보이는 까닭")

print(f"\n과녁을 S^2 = {S2:.4f} 로 두고 MSE 를 비교하면")
for name, v in (("ddof=0", biased_vars), ("ddof=1", unbiased_vars)):
    print(f"  {name}: 편향 {v.mean() - S2:+.4f}, 분산 {v.var(ddof=1):.4f}, "
          f"MSE {np.mean((v - S2) ** 2):.4f}")
print(f"  MSE 비 ddof=0 / ddof=1 = {np.mean((biased_vars - S2) ** 2) / np.mean((unbiased_vars - S2) ** 2):.4f}")

SS = biased_vars * sample_size            # 제곱합을 되살린다
print(f"\n제곱합을 c 로 나눌 때의 MSE")
for c in (98, 99, 100, 101, 102, 103):
    tag = {99: "  <- ddof=1", 100: "  <- ddof=0", 101: "  <- 최소"}.get(c, "")
    print(f"  c = {c:>3}: MSE {np.mean((SS / c - S2) ** 2):.4f}{tag}")

# 두 추정값의 분포를 겹쳐 그리고 참값 자리에 세로선을 긋는다.
# 두 분포의 모양은 같고 위치만 조금 다르다. 그 차이가 편향이다.
fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(biased_vars, bins=40, alpha=0.55, color="#1565C0",
        label="편향 (ddof=0)", density=True)
ax.hist(unbiased_vars, bins=40, alpha=0.55, color="#E65100",
        label="불편 (ddof=1)", density=True)
ax.axvline(pop_var, color="#D32F2F", linestyle="--", lw=2,
           label=f"참값 $\\sigma^2$ = {pop_var:.1f}")
ax.set_xlabel("분산 추정값")
ax.set_ylabel("밀도")
ax.set_title("모분산과 표본분산의 추정값 분포")
ax.legend(fontsize=8)
ax.spines[["top", "right"]].set_visible(False)
fig.tight_layout()
fig.savefig("descriptive_measures_advanced_179.png", dpi=170,
            facecolor="white", bbox_inches="tight")

표본마다 ddof=1 / ddof=0 의 비: 최소 1.0101010101, 최대 1.0101010101
  n/(n-1) = 1.0101010101

     분위수      ddof=0      ddof=1             비
      5%     74.3806     75.1319  1.0101010101
     25%     85.8831     86.7506  1.0101010101
     50%     94.3173     95.2700  1.0101010101
     75%    103.2273    104.2700  1.0101010101
     95%    117.4685    118.6551  1.0101010101

추정량의 표준편차 13.2114,  ddof=0 의 편향 -0.9589
  |편향| / 표준편차 = 0.0726   <- 그림에서 겹쳐 보이는 까닭

과녁을 S^2 = 95.8864 로 두고 MSE 를 비교하면
  ddof=0: 편향 -0.9807, 분산 171.0674, MSE 172.0121
  ddof=1: 편향 -0.0220, 분산 174.5408, MSE 174.5238
  MSE 비 ddof=0 / ddof=1 = 0.9856

제곱합을 c 로 나눌 때의 MSE
  c =  98: MSE 179.0175
  c =  99: MSE 174.5238  <- ddof=1
  c = 100: MSE 172.0121  <- ddof=0
  c = 101: MSE 171.3677  <- 최소
  c = 102: MSE 172.4829
  c = 103: MSE 175.2559

import numpy as np

rng = np.random.default_rng(0)
SIGMA2 = 100.0          # 참 모분산
B = 20_000

print(f"참 모분산 {SIGMA2:.1f},  모의실험 {B:,}회")
print(f"{'n':>5s}{'ddof=0 평균':>13s}{'ddof=1 평균':>13s}"
      f"{'과소추정':>10s}{'n/(n-1)':>10s}")
for n in [3, 5, 10, 20, 50, 100, 500]:
    x = rng.normal(0, np.sqrt(SIGMA2), size=(B, n))
    v0 = x.var(axis=1, ddof=0).mean()
    v1 = x.var(axis=1, ddof=1).mean()
    print(f"{n:>5d}{v0:>13.4f}{v1:>13.4f}"
          f"{1 - v0 / SIGMA2:>10.4f}{n / (n - 1):>10.4f}")

참 모분산 100.0,  모의실험 20,000회
    n    ddof=0 평균    ddof=1 평균      과소추정   n/(n-1)
    3      66.8261     100.2392    0.3317    1.5000
    5      80.6105     100.7631    0.1939    1.2500
   10      90.0330     100.0366    0.0997    1.1111
   20      95.0686     100.0722    0.0493    1.0526
   50      97.8564      99.8535    0.0214    1.0204
  100      98.9133      99.9125    0.0109    1.0101
  500      99.7507      99.9506    0.0025    1.0020

import numpy as np

rng = np.random.default_rng(0)
normal = rng.normal(0, 1, 3_000_000)
expo = rng.exponential(1, 3_000_000)

print(f"{'k':>5}{'체비셰프':>11}{'칸텔리':>11}{'정규 실제':>12}{'지수 실제':>12}")
for k in (1, 2, 3):
    print(f"{k:>5}{1 / k ** 2:>11.4f}{1 / (1 + k ** 2):>11.4f}"
          f"{np.mean(normal >= k):>12.5f}{np.mean(expo - 1 >= k):>12.5f}")

    k       체비셰프        칸텔리       정규 실제       지수 실제
    1     1.0000     0.5000     0.15852     0.13515
    2     0.2500     0.2000     0.02291     0.04948
    3     0.1111     0.1000     0.00133     0.01826

import numpy as np

k = 2.0
p = 1 / (2 * k * k)
values = np.array([-k, 0.0, k])
probs = np.array([p, 1 - 2 * p, p])

mu = (values * probs).sum()
var = ((values - mu) ** 2 * probs).sum()

print(f"값 {values}   확률 {probs.round(4)}")
print(f"  평균 {mu:.4f}   분산 {var:.4f}   표준편차 {np.sqrt(var):.4f}")
print(f"\n  P(|X - mu| >= {k} sigma) = {2 * p:.4f}")
print(f"  체비쇼프 한계 1/k^2      = {1 / k ** 2:.4f}   → 정확히 일치")

값 [-2.  0.  2.]   확률 [0.125 0.75  0.125]
  평균 0.0000   분산 1.0000   표준편차 1.0000

  P(|X - mu| >= 2.0 sigma) = 0.2500
  체비쇼프 한계 1/k^2      = 0.2500   → 정확히 일치

import numpy as np
from scipy import stats

rng = np.random.default_rng(2)
x = rng.lognormal(0, 0.8, 500_000)

print(f"AM {x.mean():.5f}  >=  GM {stats.gmean(x):.5f}  >=  HM {stats.hmean(x):.5f}")
print(f"\n옌센 격차 (볼록·오목 함수별)")
for name, g, convex in [("exp(x)", np.exp, True), ("log(x)", np.log, False),
                        ("1/x", lambda v: 1 / v, True), ("x^2", np.square, True)]:
    lhs = g(x.mean())
    rhs = g(x).mean()
    ok = (lhs <= rhs) if convex else (lhs >= rhs)
    print(f"  {name:>8}: g(E[X]) {lhs:>10.5f}   E[g(X)] {rhs:>10.5f}   "
          f"{'볼록' if convex else '오목'}  부등식 성립: {ok}")

AM 1.37840  >=  GM 1.00123  >=  HM 0.72678

옌센 격차 (볼록·오목 함수별)
    exp(x): g(E[X])    3.96853   E[g(X)] 41846194417799664.00000   볼록  부등식 성립: True
    log(x): g(E[X])    0.32092   E[g(X)]    0.00123   오목  부등식 성립: True
       1/x: g(E[X])    0.72548   E[g(X)]    1.37593   볼록  부등식 성립: True
       x^2: g(E[X])    1.89997   E[g(X)]    3.61346   볼록  부등식 성립: True