기술통계를 처음부터 구현하기: 대마 가격¶

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

# 캘리포니아와 뉴욕의 상품 대마 가격을 날짜순으로 적은 것이다.
# 두 열의 i번째 값이 같은 날에 대응하므로 나중에 상관을 구할 수 있다.
CA_PRICES = np.array([
    248.75, 248.59, 248.63, 248.37, 248.02, 247.68, 247.36,
    246.85, 246.44, 246.06, 245.81, 245.48, 245.18, 244.87,
    244.55, 244.23, 243.89, 243.60, 243.34, 243.08, 242.85,
    242.64, 242.36, 242.15, 241.88, 241.64, 241.40, 241.14,
    240.91, 240.65, 240.42, 240.20, 239.96, 239.74, 239.52,
    239.28, 239.07, 238.81, 238.55, 238.34, 238.12, 237.90,
    237.66, 237.43, 237.19, 236.98, 236.76, 236.56,
])

NY_PRICES = np.array([
    350.50, 350.31, 350.02, 349.82, 349.55, 349.30, 349.04,
    348.78, 348.54, 348.27, 348.01, 347.78, 347.51, 347.26,
    346.98, 346.72, 346.48, 346.19, 345.93, 345.68, 345.44,
    345.17, 344.93, 344.67, 344.42, 344.18, 343.91, 343.68,
    343.43, 343.17, 342.93, 342.68, 342.44, 342.18, 341.93,
    341.67, 341.43, 341.16, 340.90, 340.66, 340.41, 340.16,
    339.92, 339.67, 339.41, 339.18, 338.93, 338.70,
])

# i.i.d. 표본이면 순서가 무작위이므로 내림 횟수 D 의 평균은 (n-1)/2,
# 분산은 (n+1)/12 여야 한다. 작은 n 에서는 순열을 모두 세어 확인할 수 있다.
# 모의실험이 아니라 완전열거이므로 몬테카를로 오차가 없다.
from itertools import permutations

print(f"{'n':>3}{'E[D] 열거':>12}{'(n-1)/2':>10}"
      f"{'Var(D) 열거':>14}{'(n+1)/12':>11}{'오름 1회':>10}{'2^n-n-1':>10}")
for m in (4, 5, 6, 7):
    Ds = np.array([sum(p[i] > p[i + 1] for i in range(m - 1))
                   for p in permutations(range(m))])
    one_asc = int(np.sum(Ds == m - 2))      # 내림 m-2 회 = 오름 1 회
    print(f"{m:>3}{Ds.mean():>12.4f}{(m - 1) / 2:>10.4f}"
          f"{Ds.var():>14.4f}{(m + 1) / 12:>11.4f}{one_asc:>10d}{2 ** m - m - 1:>10d}")

# 이제 실제 자료의 D 를 재고, 시간에 대한 회귀의 r^2 도 함께 본다.
n_obs = len(CA_PRICES)
mu_D, sd_D = (n_obs - 1) / 2, np.sqrt((n_obs + 1) / 12)
print(f"\nn = 48 이면 E[D] = {mu_D}, sd(D) = {sd_D:.4f}")
for name, x in (("CA", CA_PRICES), ("NY", NY_PRICES)):
    d = np.diff(x)
    D = int(np.sum(d < 0))
    t = np.arange(1, len(x) + 1)
    r2 = stats.linregress(t, x).rvalue ** 2
    print(f"  {name}: D = {D}/{len(d)},  z = {(D - mu_D) / sd_D:5.2f}"
          f"  | 시간 추세 r^2 = {r2:.5f}")

d_ca = np.diff(CA_PRICES)
print(f"CA 에서 값이 오른 곳은 한 군데뿐: 둘째 달 -> 셋째 달 {d_ca[d_ca > 0]}")
print(f"i.i.d. 라면 완전한 내림차순일 확률 = 1/48! = {1 / math.factorial(48):.2e}")
print(f"오름이 한 번 이하일 확률 = (2^48 - 48)/48! = "
      f"{(2 ** 48 - 48) / math.factorial(48):.2e}")

# 수준의 퍼짐과 전월 대비 변화의 퍼짐을 견준다.
print(f"CA: 수준의 표준편차 {CA_PRICES.std(ddof=1):.4f}, "
      f"전월 대비 변화의 표준편차 {d_ca.std(ddof=1):.4f}")

  n     E[D] 열거   (n-1)/2     Var(D) 열거   (n+1)/12     오름 1회   2^n-n-1
  4      1.5000    1.5000        0.4167     0.4167        11        11
  5      2.0000    2.0000        0.5000     0.5000        26        26
  6      2.5000    2.5000        0.5833     0.5833        57        57
  7      3.0000    3.0000        0.6667     0.6667       120       120

n = 48 이면 E[D] = 23.5, sd(D) = 2.0207
  CA: D = 46/47,  z = 11.13  | 시간 추세 r^2 = 0.99314
  NY: D = 47/47,  z = 11.63  | 시간 추세 r^2 = 0.99994
CA 에서 값이 오른 곳은 한 군데뿐: 둘째 달 -> 셋째 달 [0.04]
i.i.d. 라면 완전한 내림차순일 확률 = 1/48! = 8.06e-62
오름이 한 번 이하일 확률 = (2^48 - 48)/48! = 2.27e-47
CA: 수준의 표준편차 3.7089, 전월 대비 변화의 표준편차 0.0771

def mean_from_scratch(data):
    """모든 값을 더해 개수로 나눈다. 정의 그대로다."""
    return np.sum(data) / len(data)


# 직접 구현한 값과 pandas의 값을 대조한다.
# 이 대조를 절마다 되풀이하는 것이 이 페이지의 요점이다.
print(f"직접 구현: {mean_from_scratch(CA_PRICES):.4f}")
print(f"pandas   : {pd.Series(CA_PRICES).mean():.4f}")

# (1) 을 확인한다. 먼저 Q 를 격자에서 최소화해 평균이 나오는지 본다.
n = len(CA_PRICES)
xbar = mean_from_scratch(CA_PRICES)
Q = lambda c: np.sum((CA_PRICES - c) ** 2)

grid = np.linspace(xbar - 2, xbar + 2, 40001)          # 간격 1e-4
Qs = ((CA_PRICES[None, :] - grid[:, None]) ** 2).sum(axis=1)
print(f"\n격자 최소점 = {grid[np.argmin(Qs)]:.6f},  평균 = {xbar:.6f}")

# 최솟값이 (n-1)s^2 인가.
print(f"Q(평균) = {Q(xbar):.4f},  (n-1)s^2 = {(n - 1) * np.var(CA_PRICES, ddof=1):.4f}")

# 분해식 Q(c) = (n-1)s^2 + n(c - xbar)^2 를 평균이 아닌 점에서 확인한다.
rhs = (n - 1) * np.var(CA_PRICES, ddof=1) + n * (242.0 - xbar) ** 2
print(f"Q(242)  = {Q(242.0):.4f},  분해식 = {rhs:.4f}"
      f"  (벌점 증가분 {n * (242.0 - xbar) ** 2:.4f})")

직접 구현: 242.3102
pandas   : 242.3102

격자 최소점 = 242.310208,  평균 = 242.310208
Q(평균) = 646.5273,  (n-1)s^2 = 646.5273
Q(242)  = 651.1463,  분해식 = 651.1463  (벌점 증가분 4.6190)

def median_from_scratch(data):
    """정렬한 뒤 가운데 값을 고른다."""
    sorted_data = np.sort(data)          # 중앙값은 순서에만 의존하므로 정렬이 먼저다
    n = len(sorted_data)
    mid = n // 2
    if n % 2 == 1:
        return sorted_data[mid]          # 홀수: 가운데 하나
    # 짝수: 가운데 두 값의 평균. 48개이므로 이 가지를 탄다.
    return (sorted_data[mid - 1] + sorted_data[mid]) / 2

print(f"직접 구현: {median_from_scratch(CA_PRICES):.4f}")
print(f"pandas   : {pd.Series(CA_PRICES).median():.4f}")

# (1) 의 주장: S 는 한 점이 아니라 구간 [x_(24), x_(25)] 전체에서 최소다.
srt = np.sort(CA_PRICES)
lo, hi = srt[23], srt[24]            # 0-기반 색인이므로 24번째와 25번째 순서통계량
S = lambda c: np.sum(np.abs(CA_PRICES - c))

print(f"\n평평한 구간 = [{lo}, {hi}],  폭 {hi - lo:.2f}")
for c in (lo, 242.0, median_from_scratch(CA_PRICES), 242.1, hi):
    print(f"  S({c:8.4f}) = {S(c):.6f}")
print("구간 밖:")
for c in (lo - 0.01, hi + 0.01):
    print(f"  S({c:8.4f}) = {S(c):.6f}")

# 닫힌 꼴 최솟값: 위쪽 24개 합 - 아래쪽 24개 합
print(f"\n위 24개 합 - 아래 24개 합 = {srt[24:].sum() - srt[:24].sum():.6f}")

직접 구현: 242.0150
pandas   : 242.0150

평평한 구간 = [241.88, 242.15],  폭 0.27
  S(241.8800) = 150.670000
  S(242.0000) = 150.670000
  S(242.0150) = 150.670000
  S(242.1000) = 150.670000
  S(242.1500) = 150.670000
구간 밖:
  S(241.8700) = 150.690000
  S(242.1600) = 150.690000

위 24개 합 - 아래 24개 합 = 150.670000

def mode_from_scratch(data, decimals=1):
    """가장 자주 나오는 값. 연속 자료이므로 먼저 반올림해 묶는다."""
    rounded = np.round(data, decimals)
    # unique가 값과 그 개수를 함께 돌려준다
    values, counts = np.unique(rounded, return_counts=True)
    return values[np.argmax(counts)]     # 개수가 가장 큰 값

# 반올림 정밀도를 바꾸면 답이 달라진다는 점을 직접 확인해 보자
for d in (0, 1, 2):
    print(f"decimals={d}: 최빈값 = {mode_from_scratch(CA_PRICES, d)}")

# 왜 그 값이 나왔는지 도수를 직접 세어 본다.
# np.unique 는 값을 오름차순으로, np.argmax 는 최대가 여럿일 때 첫째를
# 돌려주므로 동점이면 "가장 작은 값"이 최빈값이 된다.
print()
for d in (0, 1, 2):
    values, counts = np.unique(np.round(CA_PRICES, d), return_counts=True)
    top = counts.max()
    tied = values[counts == top]
    shown = ", ".join(f"{v:g}" for v in tied[:3]) + (" ..." if len(tied) > 3 else "")
    print(f"decimals={d}: 서로 다른 값 {len(values):2d}개, 최대도수 {top}, "
          f"그 도수를 갖는 값 {len(tied):2d}개 -> {shown}")

# decimals=0 의 1달러 칸마다 몇 개가 들어가는지는 하강 속도가 정한다.
step = (CA_PRICES.max() - CA_PRICES.min()) / (len(CA_PRICES) - 1)
print(f"\n한 달 평균 하강폭 {step:.4f} 달러 -> 1달러 칸마다 약 {1 / step:.2f}개")
values, counts = np.unique(np.round(CA_PRICES, 0), return_counts=True)
print(f"decimals=0 의 칸별 도수: {dict(zip(values.astype(int), counts))}")

# decimals=1 에서 248.6 으로 묶인 원자료는 무엇이고 그 값은 어디에 있는가.
print(f"248.6 으로 반올림되는 원자료: {CA_PRICES[np.round(CA_PRICES, 1) == 248.6]}")
m4, s4 = CA_PRICES.mean(), CA_PRICES.std(ddof=1)
print(f"자료 범위 [{CA_PRICES.min()}, {CA_PRICES.max()}], 평균 {m4:.2f}, s = {s4:.4f}")
print(f"248.6 은 평균보다 {(248.6 - m4) / s4:.2f} 표준편차 위, "
      f"236.56 은 {(m4 - 236.56) / s4:.2f} 표준편차 아래")

decimals=0: 최빈값 = 237.0
decimals=1: 최빈값 = 248.6
decimals=2: 최빈값 = 236.56

decimals=0: 서로 다른 값 13개, 최대도수 5, 그 도수를 갖는 값  2개 -> 237, 240
decimals=1: 서로 다른 값 47개, 최대도수 2, 그 도수를 갖는 값  1개 -> 248.6
decimals=2: 서로 다른 값 48개, 최대도수 1, 그 도수를 갖는 값 48개 -> 236.56, 236.76, 236.98 ...

한 달 평균 하강폭 0.2594 달러 -> 1달러 칸마다 약 3.86개
decimals=0 의 칸별 도수: {237: 5, 238: 4, 239: 4, 240: 5, 241: 4, 242: 4, 243: 4, 244: 3, 245: 4, 246: 3, 247: 2, 248: 3, 249: 3}
248.6 으로 반올림되는 원자료: [248.59 248.63]
자료 범위 [236.56, 248.75], 평균 242.31, s = 3.7089
248.6 은 평균보다 1.70 표준편차 위, 236.56 은 1.55 표준편차 아래

def variance_from_scratch(data):
    """평균에서의 편차를 제곱해 더하고 n-1로 나눈다."""
    m = mean_from_scratch(data)
    # n이 아니라 n-1로 나누는 것이 베셀 보정이다.
    # 편차를 "참 평균"이 아니라 "표본평균"에서 재기 때문에
    # 제곱합이 체계적으로 작아지는데, 그것을 되돌리는 보정이다.
    return np.sum((data - m) ** 2) / (len(data) - 1)

def std_from_scratch(data):
    """표준편차는 분산의 제곱근. 단위가 원자료와 같아진다."""
    return np.sqrt(variance_from_scratch(data))

# pandas의 var()/std()도 기본값이 ddof=1, 즉 n-1로 나눈다.
print(f"분산   직접 {variance_from_scratch(CA_PRICES):.4f} | "
      f"pandas {pd.Series(CA_PRICES).var():.4f}")
print(f"표준편차 직접 {std_from_scratch(CA_PRICES):.4f} | "
      f"pandas {pd.Series(CA_PRICES).std():.4f}")

# numpy는 기본이 ddof=0(n으로 나눔)이라 값이 다르다. 흔한 함정이다.
print(f"numpy 기본(ddof=0): {np.var(CA_PRICES):.4f}  <- 다르다")
print(f"numpy ddof=1      : {np.var(CA_PRICES, ddof=1):.4f}")

# (1) 의 주장: 비가 자료와 무관하게 n/(n-1) 이고 표준편차에서는 그 제곱근이다.
v0, v1 = np.var(CA_PRICES), np.var(CA_PRICES, ddof=1)
print(f"\n비 v1/v0 = {v1 / v0:.10f},  n/(n-1) = {n}/{n - 1} = {n / (n - 1):.10f}")
print(f"v0 * n/(n-1) = {v0 * n / (n - 1):.10f},  v1 = {v1:.10f}")

s0, s1 = np.std(CA_PRICES), np.std(CA_PRICES, ddof=1)
print(f"비 s1/s0 = {s1 / s0:.10f},  sqrt(n/(n-1)) = {np.sqrt(n / (n - 1)):.10f}")
print(f"어림식 1 + 1/(2(n-1)) = {1 + 1 / (2 * (n - 1)):.10f}")
print(f"ddof=0 표준편차 {s0:.6f} -> ddof=1 {s1:.6f}  (차이 {s1 - s0:.6f}, "
      f"{100 * (s1 / s0 - 1):.2f}%)")
print(f"분산은 {100 * (v1 / v0 - 1):.2f}% 차이")

분산   직접 13.7559 | pandas 13.7559
표준편차 직접 3.7089 | pandas 3.7089
numpy 기본(ddof=0): 13.4693  <- 다르다
numpy ddof=1      : 13.7559

비 v1/v0 = 1.0212765957,  n/(n-1) = 48/47 = 1.0212765957
v0 * n/(n-1) = 13.7558999557,  v1 = 13.7558999557
비 s1/s0 = 1.0105823053,  sqrt(n/(n-1)) = 1.0105823053
어림식 1 + 1/(2(n-1)) = 1.0106382979
ddof=0 표준편차 3.670057 -> ddof=1 3.708895  (차이 0.038838, 1.06%)
분산은 2.13% 차이

def covariance_from_scratch(x, y):
    """두 변수의 편차를 곱해 더하고 n-1로 나눈다.

    분산이 (x-mx)를 제곱한 것이라면, 공분산은 그 제곱을 (x-mx)(y-my)로
    바꾼 것이다. 즉 Cov(X, X) = Var(X)다.
    """
    n = len(x)
    mx, my = mean_from_scratch(x), mean_from_scratch(y)
    return np.sum((x - mx) * (y - my)) / (n - 1)

def correlation_from_scratch(x, y):
    """공분산을 각자의 표준편차로 나누어 [-1, 1]로 표준화한다."""
    return covariance_from_scratch(x, y) / (std_from_scratch(x) * std_from_scratch(y))

print(f"공분산 직접 {covariance_from_scratch(CA_PRICES, NY_PRICES):.4f} | "
      f"pandas {pd.Series(CA_PRICES).cov(pd.Series(NY_PRICES)):.4f}")
print(f"상관   직접 {correlation_from_scratch(CA_PRICES, NY_PRICES):.4f} | "
      f"pandas {pd.Series(CA_PRICES).corr(pd.Series(NY_PRICES)):.4f}")

# Cov(X, X) = Var(X) 임을 확인해 본다
print(f"Cov(CA, CA) = {covariance_from_scratch(CA_PRICES, CA_PRICES):.4f}"
      f"  =  Var(CA) = {variance_from_scratch(CA_PRICES):.4f}")

# (1) 달러를 센트로 바꾸면 공분산은 10^4 배, 상관은 그대로여야 한다.
cov_d = covariance_from_scratch(CA_PRICES, NY_PRICES)
cov_c = covariance_from_scratch(100 * CA_PRICES, 100 * NY_PRICES)
print(f"\n센트 단위(x100)")
print(f"  공분산 {cov_c:.4f}  (달러의 {cov_c / cov_d:.0f}배)")
print(f"  상관   {correlation_from_scratch(100 * CA_PRICES, 100 * NY_PRICES):.10f}"
      f"  (달러: {correlation_from_scratch(CA_PRICES, NY_PRICES):.10f})")

# (2) 같은 두 열에서 "무엇을 변수로 부르는가"만 바꾸어 r 을 네 번 구한다.
t = np.arange(1, len(CA_PRICES) + 1)
resid = lambda x, deg: x - np.polyval(np.polyfit(t, x, deg), t)
print(f"\nr(가격 수준)           = {correlation_from_scratch(CA_PRICES, NY_PRICES):.4f}")
print(f"r(1차 추세 제거 잔차)   = "
      f"{correlation_from_scratch(resid(CA_PRICES, 1), resid(NY_PRICES, 1)):.4f}")
print(f"r(2차 추세 제거 잔차)   = "
      f"{correlation_from_scratch(resid(CA_PRICES, 2), resid(NY_PRICES, 2)):.4f}")

dCA, dNY = np.diff(CA_PRICES), np.diff(NY_PRICES)
r_d = correlation_from_scratch(dCA, dNY)
m_d = len(dCA)
z_d = np.arctanh(r_d)                    # 피셔 변환으로 신뢰구간을 만든다
se_d = 1 / np.sqrt(m_d - 3)
print(f"r(전월 대비 변화)       = {r_d:.4f}")
print(f"  n = {m_d},  p = {stats.pearsonr(dCA, dNY).pvalue:.4f},  "
      f"95% 신뢰구간 [{np.tanh(z_d - 1.96 * se_d):.4f}, "
      f"{np.tanh(z_d + 1.96 * se_d):.4f}]")

공분산 직접 13.1002 | pandas 13.1002
상관   직접 0.9970 | pandas 0.9970
Cov(CA, CA) = 13.7559  =  Var(CA) = 13.7559

센트 단위(x100)
  공분산 131002.0980  (달러의 10000배)
  상관   0.9970102134  (달러: 0.9970102134)

r(가격 수준)           = 0.9970
r(1차 추세 제거 잔차)   = 0.7417
r(2차 추세 제거 잔차)   = 0.3257
r(전월 대비 변화)       = 0.0515
  n = 47,  p = 0.7310,  95% 신뢰구간 [-0.2392, 0.3337]

data = CA_PRICES
s = pd.Series(data)

# 뒤의 시각화에서 쓰려고 평균과 중앙값을 이름에 담아 둔다
m   = mean_from_scratch(data)
med = median_from_scratch(data)

# 직접 구현한 함수와 pandas 메서드를 한 줄씩 나란히 찍는다.
checks = [
    ("mean",   m,                           s.mean()),
    ("median", med,                          s.median()),
    ("var",    variance_from_scratch(data), s.var()),
    ("std",    std_from_scratch(data),      s.std()),
]

print(f"{'stat':<8}{'from scratch':>14}{'pandas':>14}{'diff':>12}")
for name, mine, theirs in checks:
    print(f"{name:<8}{mine:>14.6f}{theirs:>14.6f}{mine - theirs:>12.2e}")

# 두 지역 사이의 관계
cov  = covariance_from_scratch(CA_PRICES, NY_PRICES)
corr = correlation_from_scratch(CA_PRICES, NY_PRICES)
print(f"\n공분산 = {cov:.4f}")
print(f"상관   = {corr:.4f}")

# 평균의 일치는 순서 덕인가? 자료를 섞어 가며 200번 비교해 본다.
rng = np.random.default_rng(0)
bad = sum(np.sum(p) / len(p) != pd.Series(p).mean()
          for p in (rng.permutation(data) for _ in range(200)))
print(f"\n자료 순서를 200번 섞어 평균을 비교: 어긋난 횟수 {bad}")

# 분산은 어느 쪽이 얼마나 틀렸는가. math.fsum 은 정확히 반올림된 합을 준다.
dev2 = (data - m) ** 2
n = len(data)
v_exact = math.fsum(dev2) / (n - 1)           # 기준값
v_np = np.sum(dev2) / (n - 1)                 # np.sum: 쌍합(pairwise summation)
v_loop = 0.0
for z in dev2:                                # 순차 덧셈
    v_loop += z
v_loop /= n - 1

ulp = np.spacing(v_exact)
eps = np.finfo(float).eps
print(f"\n1 ULP = {ulp:.3e},  eps = {eps:.3e}")
print(f"{'구현':<8}{'분산':>22}{'ULP':>8}{'상대오차/eps':>14}")
for name, v in (("fsum", v_exact), ("np.sum", v_np),
                ("loop", v_loop), ("pandas", s.var())):
    print(f"{name:<8}{v!r:>22}{(v - v_exact) / ulp:>8.1f}"
          f"{(v - v_exact) / v_exact / eps:>14.2f}")
print(f"np.sum 과 pandas 의 간격 = {(v_np - s.var()) / ulp:.1f} ULP = "
      f"{v_np - s.var():.2e}")

stat      from scratch        pandas        diff
mean        242.310208    242.310208    0.00e+00
median      242.015000    242.015000    0.00e+00
var          13.755900     13.755900   -5.33e-15
std           3.708895      3.708895   -8.88e-16

공분산 = 13.1002
상관   = 0.9970

자료 순서를 200번 섞어 평균을 비교: 어긋난 횟수 0

1 ULP = 1.776e-15,  eps = 2.220e-16
구현                          분산     ULP      상대오차/eps
fsum        13.755899955673765     0.0          0.00
np.sum      13.755899955673762    -2.0         -1.16
loop        13.755899955673767     1.0          0.58
pandas      13.755899955673767     1.0          0.58
np.sum 과 pandas 의 간격 = -3.0 ULP = -5.33e-15

# 그림에서 읽을 수치를 먼저 찍어 둔다.
counts, edges = np.histogram(data, bins=15)
print(f"히스토그램 15칸: 칸 너비 {edges[1] - edges[0]:.4f}, "
      f"칸마다 기대 도수 {len(data) / 15:.2f}")
print(f"  실제 도수 {counts}")
print(f"평균 {m:.4f}  중앙값 {med:.4f}  평균-중앙값 = {m - med:.4f}")
print(f"왜도 g1 = {stats.skew(data):.4f}  초과첨도 = {stats.kurtosis(data):.4f}")
q1, q2, q3 = np.percentile(data, [25, 50, 75])
iqr = q3 - q1
print(f"Q1 {q1:.4f}  Q2 {q2:.4f}  Q3 {q3:.4f}  IQR {iqr:.4f}")
print(f"울타리 [{q1 - 1.5 * iqr:.4f}, {q3 + 1.5 * iqr:.4f}], "
      f"자료 범위 [{data.min()}, {data.max()}]")
print(f"같은 범위의 균등분포라면 s = (max-min)/sqrt(12) = "
      f"{(data.max() - data.min()) / np.sqrt(12):.4f}  (실제 s = {std_from_scratch(data):.4f})")
print(f"월별 변화의 범위 [{np.diff(data).min():.2f}, {np.diff(data).max():.2f}], "
      f"평균 {np.diff(data).mean():.4f}")
print(f"맨 위 칸에 든 값 {np.sort(data[data >= edges[14]])} "
      f"(그 안의 간격 {np.round(np.diff(np.sort(data[data >= edges[14]])), 2)})")

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

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 왼쪽: 히스토그램에 평균과 중앙값을 세로선으로 얹는다.
axes[0].hist(data, bins=15, color="#DCEBFB", edgecolor="#1565C0", alpha=0.9)
axes[0].axvline(m, color="#E65100", linestyle="--", lw=2, label=f"평균 {m:.1f}")
axes[0].axvline(med, color="#33691E", linestyle=":", lw=2, label=f"중앙값 {med:.1f}")
axes[0].set_title("캘리포니아 가격 분포")
axes[0].set_xlabel("가격 (달러)")
axes[0].set_ylabel("도수")
axes[0].legend(fontsize=8)

# 가운데: 같은 자료의 상자그림. 다섯 수치 요약을 한 눈에 보여 준다.
axes[1].boxplot(data, vert=True, labels=["캘리포니아"])
axes[1].set_title("상자그림 — 캘리포니아")
axes[1].set_ylabel("가격 (달러)")

# 오른쪽: 두 주의 가격을 짝지어 찍는다. 점이 직선에 가까울수록 상관이 크다.
axes[2].scatter(CA_PRICES, NY_PRICES, alpha=0.7, color="#1565C0")
axes[2].set_xlabel("캘리포니아 가격 (달러)")
axes[2].set_ylabel("뉴욕 가격 (달러)")
axes[2].set_title(f"캘리포니아 대 뉴욕  ($r$ = {corr:.3f})")

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

히스토그램 15칸: 칸 너비 0.8127, 칸마다 기대 도수 3.20
  실제 도수 [4 4 3 4 3 4 3 3 3 3 3 2 2 2 5]
평균 242.3102  중앙값 242.0150  평균-중앙값 = 0.2952
왜도 g1 = 0.2091  초과첨도 = -1.1259
Q1 239.2275  Q2 242.0150  Q3 245.2550  IQR 6.0275
울타리 [230.1862, 254.2962], 자료 범위 [236.56, 248.75]
같은 범위의 균등분포라면 s = (max-min)/sqrt(12) = 3.5189  (실제 s = 3.7089)
월별 변화의 범위 [-0.51, 0.04], 평균 -0.2594
맨 위 칸에 든 값 [248.02 248.37 248.59 248.63 248.75] (그 안의 간격 [0.35 0.22 0.04 0.12])

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
walk = lambda n: np.cumsum(rng.normal(0, 1, n))     # 서로 완전히 독립인 확률보행

print("독립인 두 확률보행의 표본상관 (참 상관은 0)")
print(f"{'n':>7}{'|r| > 0.5 비율':>16}{'p < 0.05 비율':>16}")
for n in (50, 200, 1000):
    rs = np.array([np.corrcoef(walk(n), walk(n))[0, 1] for _ in range(4000)])
    ps = np.array([stats.pearsonr(walk(n), walk(n))[1] for _ in range(2000)])
    print(f"{n:>7}{np.mean(np.abs(rs) > 0.5):>16.4f}{np.mean(ps < 0.05):>16.4f}")

독립인 두 확률보행의 표본상관 (참 상관은 0)
      n    |r| > 0.5 비율     p < 0.05 비율
     50          0.4115          0.6455
    200          0.3987          0.8245
   1000          0.3987          0.9260

import numpy as np

x = np.array([10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5], float)
datasets = {
    "I":   (x, np.array([8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68])),
    "II":  (x, np.array([9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74])),
    "III": (x, np.array([7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73])),
    "IV":  (np.array([8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8], float),
            np.array([6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89])),
}

print(f"{'집합':>5}{'x 평균':>9}{'y 평균':>9}{'x 분산':>9}{'y 분산':>9}{'상관':>9}{'기울기':>9}{'절편':>9}")
for name, (xx, yy) in datasets.items():
    b, a = np.polyfit(xx, yy, 1)
    print(f"{name:>5}{xx.mean():>9.3f}{yy.mean():>9.3f}{xx.var(ddof=1):>9.3f}"
          f"{yy.var(ddof=1):>9.3f}{np.corrcoef(xx, yy)[0, 1]:>9.3f}{b:>9.3f}{a:>9.3f}")

   집합     x 평균     y 평균     x 분산     y 분산       상관      기울기       절편
    I    9.000    7.501   11.000    4.127    0.816    0.500    3.000
   II    9.000    7.501   11.000    4.128    0.816    0.500    3.001
  III    9.000    7.500   11.000    4.123    0.816    0.500    3.002
   IV    9.000    7.501   11.000    4.123    0.817    0.500    3.002

import numpy as np

rng = np.random.default_rng(0)
n = 200_000
X = rng.normal(0, 1, n)
Y = 0.8 * X + rng.normal(0, 0.6, n)
rho_true = np.corrcoef(X, Y)[0, 1]
print(f"참 상관 {rho_true:.4f}\n")

print(f"{'측정오차 sd':>12}{'신뢰도':>10}{'관측 상관':>12}{'이론 감쇠':>12}")
for se in (0.0, 0.5, 1.0, 2.0):
    X_obs = X + rng.normal(0, se, n)
    reliability = 1 / (1 + se ** 2)
    print(f"{se:>12.1f}{reliability:>10.4f}"
          f"{np.corrcoef(X_obs, Y)[0, 1]:>12.4f}{rho_true * np.sqrt(reliability):>12.4f}")

참 상관 0.7993

     측정오차 sd       신뢰도       관측 상관       이론 감쇠
         0.0    1.0000      0.7993      0.7993
         0.5    0.8000      0.7143      0.7149
         1.0    0.5000      0.5654      0.5652
         2.0    0.2000      0.3590      0.3575

import numpy as np

print(f"{'더한 상수':>12}{'교과서 공식':>18}{'올바른 값':>14}")
for off in (0, 1e6, 1e8, 1e9):
    x = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + off
    n = len(x)
    naive = (np.sum(x ** 2) - np.sum(x) ** 2 / n) / (n - 1)
    good = np.var(x, ddof=1)
    print(f"{off:>12.0e}{naive:>18.8f}{good:>14.8f}")

       더한 상수            교과서 공식         올바른 값
       0e+00        2.50000000    2.50000000
       1e+06        2.50000000    2.50000000
       1e+08        2.00000000    2.50000000
       1e+09        0.00000000    2.50000000