육각형 구간 그림과 2차원 히스토그램¶

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(0)
n = 20000
x = rng.normal(0, 1, n)
y = 0.7 * x + rng.normal(0, 0.7, n)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# hexbin: 평면을 정육각형으로 덮는다.
#   gridsize=30  가로 방향으로 육각형을 30개 놓는다 (많을수록 잘게 나뉨)
#   mincnt=1     점이 하나도 없는 칸은 그리지 않는다 (배경이 흰색으로 남음)
hb = axes[0].hexbin(x, y, gridsize=30, cmap='viridis', mincnt=1)
axes[0].set_title("hexbin — 육각형 칸")
plt.colorbar(hb, ax=axes[0], label="개수")

# hist2d: 평면을 정사각형 격자로 나눈다.
#   bins=30   가로·세로 각각 30등분 -> 900개의 칸
#   cmin=1    개수가 1 미만인 칸은 그리지 않는다 (mincnt와 같은 역할)
# hist2d는 (도수, x경계, y경계, 이미지)를 돌려주므로 [3]번이 색막대용 객체다
h = axes[1].hist2d(x, y, bins=30, cmap='viridis', cmin=1)
axes[1].set_title("hist2d — 사각형 칸")
plt.colorbar(h[3], ax=axes[1], label="개수")

for ax in axes:
    ax.set_xlabel('x')
    ax.set_ylabel('y')
plt.tight_layout()
plt.show()

# (1) 평면을 덮는 정다각형은 셋뿐이다. 내각 (k-2)*180/k 가 360 을 나누어야 한다.
print(f"{'k':>4}{'내각':>10}{'360/내각':>11}{'평면을 덮는가':>15}")
for k in range(3, 13):
    ang = (k - 2) * 180 / k
    q = 360 / ang
    print(f"{k:>4}{ang:>10.3f}{q:>11.4f}{'예' if abs(q - round(q)) < 1e-12 else '아니오':>15}")

# (2) 같은 면적에서 둘레가 짧은 타일이 원에 가깝다.
print(f"\n{'타일':>8}{'둘레/sqrt(면적)':>18}")
for name, per, area in [("정삼각형", 3.0, np.sqrt(3) / 4),
                        ("정사각형", 4.0, 1.0),
                        ("정육각형", 6.0, 3 * np.sqrt(3) / 2)]:
    print(f"{name:>8}{per / np.sqrt(area):>18.4f}")
print(f"{'원':>8}{2 * np.sqrt(np.pi):>18.4f}   <- 이론적 하한")

# 이웃까지의 거리. 변 길이를 1 로 둔 격자의 중심 간 거리를 센다.
print("\n이웃까지의 중심 간 거리 (변 길이 1)")
sq = np.array([1.0, 1.0, 1.0, 1.0, np.sqrt(2), np.sqrt(2), np.sqrt(2), np.sqrt(2)])
hx = np.full(6, np.sqrt(3))
print(f"  정사각형 격자: 이웃 {len(sq)}개  거리 {np.unique(np.round(sq, 4))}"
      f"  최대/최소 {sq.max()/sq.min():.4f}")
print(f"  정육각형 격자: 이웃 {len(hx)}개  거리 {np.unique(np.round(hx, 4))}"
      f"  최대/최소 {hx.max()/hx.min():.4f}")

# 가장 진한 칸의 도수
H, _, _ = np.histogram2d(x, y, bins=30)
print(f"\n가장 진한 칸: 육각형 {int(hb.get_array().max())}개, 사각형 {int(H.max())}개")

   k        내각     360/내각        평면을 덮는가
   3    60.000     6.0000              예
   4    90.000     4.0000              예
   5   108.000     3.3333            아니오
   6   120.000     3.0000              예
   7   128.571     2.8000            아니오
   8   135.000     2.6667            아니오
   9   140.000     2.5714            아니오
  10   144.000     2.5000            아니오
  11   147.273     2.4444            아니오
  12   150.000     2.4000            아니오

      타일       둘레/sqrt(면적)
    정삼각형            4.5590
    정사각형            4.0000
    정육각형            3.7224
       원            3.5449   <- 이론적 하한

이웃까지의 중심 간 거리 (변 길이 1)
  정사각형 격자: 이웃 8개  거리 [1.     1.4142]  최대/최소 1.4142
  정육각형 격자: 이웃 6개  거리 [1.7321]  최대/최소 1.0000

가장 진한 칸: 육각형 270개, 사각형 320개

import warnings
warnings.filterwarnings("ignore")

import os
import numpy as np
import pandas as pd
import yfinance as yf

# ── 내려받을 대상과 기간 ────────────────────────────────────────────
# ^GSPC 는 S&P 500 지수, AAPL 은 애플 보통주다.
# 기간을 날짜로 고정해 두었다. datetime.now() 를 쓰면 실행할 때마다
# 자료가 달라져 아래 출력이 재현되지 않는다.
TICKERS = {"^GSPC": "GSPC", "AAPL": "AAPL"}
START, END = "2025-09-01", "2026-09-01"
CACHE = "datasets/stocks"          # 한 번 받은 자료를 저장해 둘 곳

os.makedirs(CACHE, exist_ok=True)

def load(ticker: str, name: str) -> pd.DataFrame:
    """자료를 캐시에서 읽고, 없으면 내려받아 캐시에 저장한다.

    야후 파이낸스는 짧은 시간에 여러 번 요청하면 거절한다(rate limit).
    한 번 받아 두고 다시 쓰는 편이 안전하고 빠르다.
    """
    path = f"{CACHE}/{name}.csv"
    if os.path.exists(path):
        return pd.read_csv(path, index_col=0, parse_dates=True)

    try:
        # curl_cffi 가 있으면 브라우저처럼 행세해 차단을 덜 받는다.
        from curl_cffi import requests as cr
        session = cr.Session(impersonate="chrome")
    except ImportError:
        session = None

    df = yf.Ticker(ticker, session=session).history(
        start=START, end=END, auto_adjust=False)
    df.to_csv(path)
    return df

frames = {name: load(tk, name) for tk, name in TICKERS.items()}
for name, df in frames.items():
    print(f"{name:5s}  {df.shape[0]}일  "
          f"{df.index.min().date()} ~ {df.index.max().date()}")

# ── 일별 수익률 ─────────────────────────────────────────────────────
# 배당·액면분할을 반영한 'Adj Close' 가 있으면 그것을 쓴다.
# pct_change() 는 (오늘값 - 어제값) / 어제값 이고, 첫날은 NaN 이라 버린다.
# 100 을 곱해 퍼센트 단위로 바꾼다.
def returns(df: pd.DataFrame) -> pd.Series:
    col = "Adj Close" if "Adj Close" in df.columns else "Close"
    return df[col].pct_change().dropna() * 100

r = {name: returns(df) for name, df in frames.items()}

# 두 종목의 거래일이 완전히 같다는 보장이 없으므로 교집합만 남긴다.
common = r["GSPC"].index.intersection(r["AAPL"].index)
x = r["GSPC"].loc[common].to_numpy()      # S&P 500 일별 수익률 (%)
y = r["AAPL"].loc[common].to_numpy()      # AAPL    일별 수익률 (%)

print(f"\n공통 거래일 {len(x)}일\n")
print(f"{'':6s}{'평균':>9s}{'표준편차':>10s}{'최소':>9s}{'최대':>9s}")
for lab, v in [("S&P", x), ("AAPL", y)]:
    # ddof=1 은 표본표준편차다. numpy 의 기본값은 ddof=0 이므로 명시한다.
    print(f"{lab:6s}{v.mean():>+9.4f}{v.std(ddof=1):>10.4f}"
          f"{v.min():>+9.4f}{v.max():>+9.4f}")
print(f"\n피어슨 상관계수 {np.corrcoef(x, y)[0, 1]:.4f}")

# ── (1) 의 세 값 ────────────────────────────────────────────────────
rho = np.corrcoef(x, y)[0, 1]
sx, sy = x.std(ddof=1), y.std(ddof=1)
beta = rho * sy / sx
s_eps = sy * np.sqrt(1 - rho ** 2)

print(f"\n베타  = rho * sy / sx = {rho:.4f} * {sy:.4f} / {sx:.4f} = {beta:.4f}")
# 최소제곱으로 직접 적합해 공식을 확인한다.
b_ols, a_ols = np.polyfit(x, y, 1)
print(f"최소제곱 기울기 {b_ols:.4f}   절편 {a_ols:+.4f}")
print(f"R^2 = rho^2 = {rho**2:.4f}   →  지수가 설명하는 분산 {rho**2:.1%}, "
      f"고유 변동 {1 - rho**2:.1%}")
print(f"고유 변동의 표준편차 s_eps = sy * sqrt(1 - rho^2) = {s_eps:.4f}")
resid = y - (a_ols + b_ols * x)
print(f"  실제 잔차의 표준편차 (ddof=1)             = {resid.std(ddof=1):.4f}")
print(f"분산 분해: beta^2 * sx^2 + s_eps^2 = {beta**2 * sx**2 + s_eps**2:.4f}"
      f"   sy^2 = {sy**2:.4f}")

GSPC   251일  2025-09-02 ~ 2026-08-31
AAPL   251일  2025-09-02 ~ 2026-08-31

공통 거래일 250일

             평균      표준편차       최소       최대
S&P     +0.0755    0.8061  -2.7112  +2.9131
AAPL    +0.1427    1.5832  -7.3539  +4.8407

피어슨 상관계수 0.3704

베타  = rho * sy / sx = 0.3704 * 1.5832 / 0.8061 = 0.7274
최소제곱 기울기 0.7274   절편 +0.0878
R^2 = rho^2 = 0.1372   →  지수가 설명하는 분산 13.7%, 고유 변동 86.3%
고유 변동의 표준편차 s_eps = sy * sqrt(1 - rho^2) = 1.4706
  실제 잔차의 표준편차 (ddof=1)             = 1.4706
분산 분해: beta^2 * sx^2 + s_eps^2 = 2.5066   sy^2 = 2.5066

import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

# ── 2차원 정규분포를 자료에 맞춘다 ──────────────────────────────────
# 적합이라야 할 일이 두 줄뿐이다. 2차원 정규분포의 모수는
# 평균벡터(2개)와 공분산행렬(3개: 분산 2개 + 공분산 1개)이고,
# 최대가능도 추정값이 곧 표본평균과 표본공분산이기 때문이다.
X = np.column_stack([x, y])
mu = X.mean(axis=0)
S = np.cov(X.T)                    # np.cov 는 기본이 ddof=1
model = multivariate_normal(mean=mu, cov=S)

print(f"평균벡터   [{mu[0]:+.4f}, {mu[1]:+.4f}]")
print("공분산행렬")
print(f"   [{S[0, 0]:8.4f} {S[0, 1]:8.4f}]")
print(f"   [{S[1, 0]:8.4f} {S[1, 1]:8.4f}]")
print(f"상관계수   {S[0, 1] / np.sqrt(S[0, 0] * S[1, 1]):.4f}")
print(f"평균 로그가능도 {model.logpdf(X).mean():.4f}")

# (1)(2) 두 분모를 나란히 견준다.
S_mle = np.cov(X.T, ddof=0)                    # 최대가능도추정값: 분모 n
nobs, dim = X.shape
print(f"\n{'분모':>8}{'det':>12}{'D^2 평균':>12}{'평균 로그가능도':>18}{'상관계수':>11}")
for lab, Sg in [("n-1", S), ("n", S_mle)]:
    dev = X - mu
    D2g = np.einsum("ij,jk,ik->i", dev, np.linalg.inv(Sg), dev)
    llg = multivariate_normal(mean=mu, cov=Sg).logpdf(X).mean()
    # 닫힌 꼴:  -d/2 log(2pi) - 1/2 log det - (D^2 의 평균)/2
    closed = (-dim / 2 * np.log(2 * np.pi) - 0.5 * np.log(np.linalg.det(Sg))
              - D2g.mean() / 2)
    print(f"{lab:>8}{np.linalg.det(Sg):>12.6f}{D2g.mean():>12.6f}"
          f"{llg:>18.8f}{Sg[0,1]/np.sqrt(Sg[0,0]*Sg[1,1]):>11.6f}")
    print(f"{'':8}{'':12}{'':12}{'닫힌 꼴 ' + f'{closed:.8f}':>18}")
c = nobs / (nobs - 1)
print(f"\n최대가능도쪽이 큰 양  "
      f"{multivariate_normal(mean=mu, cov=S_mle).logpdf(X).mean() - model.logpdf(X).mean():.3e}")
print(f"예측값 log c + 1/c - 1  (c = {c:.6f})  {np.log(c) + 1/c - 1:.3e}")

# ── PDF 를 그릴 격자 ────────────────────────────────────────────────
# 자료 범위보다 1%p 씩 넓게 잡아 등고선이 잘리지 않게 한다.
gx = np.linspace(x.min() - 1, x.max() + 1, 120)
gy = np.linspace(y.min() - 1, y.max() + 1, 120)
GX, GY = np.meshgrid(gx, gy)
# dstack 으로 (120, 120, 2) 모양을 만들면 pdf 가 격자 전체를 한 번에 받는다.
Z = model.pdf(np.dstack([GX, GY]))

fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.8))

# ── 왼쪽: 관측 자료의 2차원 히스토그램 ──────────────────────────────
# cmin=1 은 비어 있는 칸을 그리지 않는다(흰색으로 남긴다).
h = axes[0].hist2d(x, y, bins=35, cmap="Blues", cmin=1)
axes[0].set_title("2차원 히스토그램 (관측)")
fig.colorbar(h[3], ax=axes[0], label="개수")

# ── 가운데: 적합된 모형의 확률밀도 ──────────────────────────────────
cf = axes[1].contourf(GX, GY, Z, levels=20, cmap="RdYlBu_r", alpha=0.85)
cl = axes[1].contour(GX, GY, Z, levels=10, colors="black",
                     alpha=0.35, linewidths=0.7)
axes[1].clabel(cl, inline=True, fontsize=7, fmt="%.3f")
axes[1].scatter(*mu, color="darkred", s=140, marker="X", edgecolors="black",
                linewidth=1.5, zorder=5, label="평균")
axes[1].set_title("적합된 2차원 정규 확률밀도 (모형)")
fig.colorbar(cf, ax=axes[1], label="확률밀도")
axes[1].legend(loc="upper left")

# ── 오른쪽: 둘을 겹쳐서 ─────────────────────────────────────────────
# 히스토그램을 흐리게(alpha=0.65) 깔고 모형의 등고선을 빨간 선으로 올린다.
h2 = axes[2].hist2d(x, y, bins=35, cmap="Blues", cmin=1, alpha=0.65)
co = axes[2].contour(GX, GY, Z, levels=12, colors="red",
                     alpha=0.8, linewidths=1.2)
axes[2].clabel(co, inline=True, fontsize=7, fmt="%.3f")
axes[2].scatter(*mu, color="red", s=140, marker="X", edgecolors="darkred",
                linewidth=2, zorder=5, label="평균")
axes[2].set_title("겹쳐 그리기 (빨강 = 모형)")
fig.colorbar(h2[3], ax=axes[2], label="개수")
axes[2].legend(loc="upper left")

for ax in axes:
    ax.set_xlabel("S&P 500 일별 수익률 (%)")
    ax.set_ylabel("AAPL 일별 수익률 (%)")
    ax.grid(alpha=0.2)

fig.suptitle("S&P 500 과 AAPL 의 일별 수익률, 2025-09 ~ 2026-08", y=1.02)
fig.tight_layout()
plt.show()

평균벡터   [+0.0755, +0.1427]
공분산행렬
   [  0.6499   0.4727]
   [  0.4727   2.5066]
상관계수   0.3704
평균 로그가능도 -3.0041

      분모         det      D^2 평균          평균 로그가능도       상관계수
     n-1    1.405476    1.992000       -3.00406518   0.370399
                                  닫힌 꼴 -3.00406518
       n    1.394255    2.000000       -3.00405715   0.370399
                                  닫힌 꼴 -3.00405715

최대가능도쪽이 큰 양  8.021e-06
예측값 log c + 1/c - 1  (c = 1.004016)  8.021e-06

# 가장 진한 칸에 몇 개가 들어가는지부터 센다.
print(f"{'bins':>6s}{'칸 수':>8s}{'평균 도수':>10s}{'최대 도수':>10s}"
      f"{'빈 칸 비율':>11s}{'균등가정':>10s}{'정확':>8s}")
for b in [35, 20, 12, 8]:
    H, _, _ = np.histogram2d(x, y, bins=b)
    lam = len(x) / b**2
    print(f"{b:>6d}{b * b:>8d}{lam:>10.3f}"
          f"{int(H.max()):>10d}{(H == 0).mean():>11.3f}"
          f"{np.exp(-lam):>10.3f}{(1 - 1/b**2)**len(x):>8.3f}")
print(f"\n평균 도수 5 를 보장하는 상한 sqrt(n/5) = {np.sqrt(len(x)/5):.3f}")

fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.5))

h1 = axes[0].hist2d(x, y, bins=35, cmap="Blues", cmin=1)
axes[0].set_title("bins=35 — 칸이 너무 잘다")
fig.colorbar(h1[3], ax=axes[0], label="개수")

h2 = axes[1].hist2d(x, y, bins=12, cmap="Blues", cmin=1)
axes[1].set_title("bins=12 — 읽을 만하다")
fig.colorbar(h2[3], ax=axes[1], label="개수")

# hexbin 은 같은 자료를 육각형 칸으로 덮는다.
hb = axes[2].hexbin(x, y, gridsize=12, cmap="Blues", mincnt=1)
axes[2].set_title("hexbin, gridsize=12")
fig.colorbar(hb, ax=axes[2], label="개수")

for ax in axes:
    ax.set_xlabel("S&P 500 일별 수익률 (%)")
    ax.set_ylabel("AAPL 일별 수익률 (%)")
    ax.grid(alpha=0.2)
fig.suptitle("같은 250일, 세 가지 칸 크기", y=1.02)
fig.tight_layout()
plt.show()

  bins     칸 수     평균 도수     최대 도수     빈 칸 비율      균등가정      정확
    35    1225     0.204         6      0.869     0.815   0.815
    20     400     0.625        11      0.752     0.535   0.535
    12     144     1.736        26      0.639     0.176   0.175
     8      64     3.906        43      0.500     0.020   0.020

평균 도수 5 를 보장하는 상한 sqrt(n/5) = 7.071

from scipy import stats

# ── 마할라노비스 거리 ───────────────────────────────────────────────
d = X - mu
D2 = np.einsum("ij,jk,ik->i", d, np.linalg.inv(S), d)

n = len(D2)
# 관측된 D² 를 정렬하고, 같은 개수의 이론 분위수와 짝짓는다.
q_theory = stats.chi2.ppf((np.arange(1, n + 1) - 0.5) / n, df=2)
D2_sorted = np.sort(D2)

# (1) D² 의 평균은 자료와 무관하게 d(n-1)/n 이다. 다른 자료로도 확인한다.
print(f"d(n-1)/n = 2 * {n-1}/{n} = {2*(n-1)/n:.6f}")
rng_chk = np.random.default_rng(3)
for trial in range(3):
    Z = rng_chk.standard_t(3, size=(n, 2))     # 꼬리가 무거운 t(3) 자료
    mz, Sz = Z.mean(0), np.cov(Z.T)
    dz = Z - mz
    print(f"  t(3) 자료 {trial + 1}: D² 의 평균 "
          f"{np.einsum('ij,jk,ik->i', dz, np.linalg.inv(Sz), dz).mean():.6f}")

print(f"\nD² 의 평균   {D2.mean():.4f}   (이론값 2)")
print(f"D² 의 95분위 {np.quantile(D2, 0.95):.4f}   "
      f"(이론값 {stats.chi2.ppf(0.95, 2):.4f})")
print(f"D² 의 99분위 {np.quantile(D2, 0.99):.4f}   "
      f"(이론값 {stats.chi2.ppf(0.99, 2):.4f})")
print(f"D² 최댓값    {D2.max():.4f}")

# 95% 타원 밖으로 나간 날이 몇 %인가 (이론상 5%)
print(f"\n95% 타원 밖 관측 비율 {np.mean(D2 > stats.chi2.ppf(0.95, 2)):.4f} (이론 0.05)")
print(f"99% 타원 밖 관측 비율 {np.mean(D2 > stats.chi2.ppf(0.99, 2)):.4f} (이론 0.01)")

# ── 꼬리 확률: 모형과 실제가 얼마나 다른가 ──────────────────────────
# 두 자산이 "함께 크게 빠지는 날"이 투자자에게 가장 중요한데,
# 2차원 정규분포가 그 확률을 제대로 맞히는지 직접 세어 본다.
# 모형 확률은 결합 CDF 로 바로 구한다. 모의실험으로 근사하면
# 배율이 표본오차 때문에 몇 %씩 흔들리므로 정확한 값을 쓰는 편이 낫다.
print(f"\n{'사건':>26s}{'실제':>9s}{'모형':>9s}{'배율':>8s}")
for lab, thr in [("둘 다 -1% 아래", -1.0), ("둘 다 -2% 아래", -2.0)]:
    emp = np.mean((x < thr) & (y < thr))
    mod = model.cdf([thr, thr])
    print(f"{lab:>26s}{emp:>9.4f}{mod:>9.4f}{emp / mod:>8.2f}")

# ── 주변분포의 첨도 ─────────────────────────────────────────────────
print(f"\n초과첨도 (정규분포는 0)")
for lab, v in [("S&P 500", x), ("AAPL", y)]:
    print(f"  {lab:9s} {stats.kurtosis(v):+.4f}"
          f"   왜도 {stats.skew(v):+.4f}")

# ── 그림 ────────────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 3, figsize=(16.5, 4.5))

axes[0].scatter(q_theory, D2_sorted, s=14, alpha=0.7, color="steelblue")
lim = max(q_theory.max(), D2_sorted.max()) * 1.05
axes[0].plot([0, lim], [0, lim], "r--", lw=1.5, label="y = x")
axes[0].set_xlabel("자유도 2인 카이제곱 분위수")
axes[0].set_ylabel("관측된 $D^2$")
axes[0].set_title("마할라노비스 QQ 그림")
axes[0].legend()

# 95% 타원을 자료 위에 직접 그려 본다.
# 단위원을 촐레스키 인자로 늘려 주면 곧 등고선 타원이 된다.
theta = np.linspace(0, 2 * np.pi, 200)
circle = np.column_stack([np.cos(theta), np.sin(theta)])
L = np.linalg.cholesky(S)
for lvl, c in [(0.50, "tab:green"), (0.95, "tab:orange"), (0.99, "tab:red")]:
    e = mu + (circle * np.sqrt(stats.chi2.ppf(lvl, 2))) @ L.T
    axes[1].plot(e[:, 0], e[:, 1], color=c, lw=1.8, label=f"{lvl:.0%} 타원")
outside = D2 > stats.chi2.ppf(0.95, 2)
axes[1].scatter(x[~outside], y[~outside], s=12, alpha=0.5, color="steelblue")
axes[1].scatter(x[outside], y[outside], s=28, color="crimson",
                edgecolors="black", linewidth=0.4, label="95% 타원 밖")
axes[1].set_xlabel("S&P 500 일별 수익률 (%)")
axes[1].set_ylabel("AAPL 일별 수익률 (%)")
axes[1].set_title("자료와 모형의 등고선 타원")
axes[1].legend(fontsize=8)

# 주변분포의 정규 QQ 그림 (AAPL)
# probplot 은 축 이름을 영어로 붙이므로 그린 뒤에 덮어쓴다.
stats.probplot(y, dist="norm", plot=axes[2])
axes[2].set_title("AAPL 주변분포의 정규 QQ 그림")
axes[2].set_xlabel("이론 분위수"); axes[2].set_ylabel("관측 분위수")
axes[2].get_lines()[0].set_markersize(4)

for ax in axes:
    ax.grid(alpha=0.2)
fig.suptitle("2차원 정규분포는 잘 맞는가?", y=1.02)
fig.tight_layout()
plt.show()

d(n-1)/n = 2 * 249/250 = 1.992000
  t(3) 자료 1: D² 의 평균 1.992000
  t(3) 자료 2: D² 의 평균 1.992000
  t(3) 자료 3: D² 의 평균 1.992000

D² 의 평균   1.9920   (이론값 2)
D² 의 95분위 6.2782   (이론값 5.9915)
D² 의 99분위 12.8549   (이론값 9.2103)
D² 최댓값    29.8324

95% 타원 밖 관측 비율 0.0560 (이론 0.05)
99% 타원 밖 관측 비율 0.0280 (이론 0.01)

                        사건       실제       모형      배율
                둘 다 -1% 아래   0.0440   0.0434    1.01
                둘 다 -2% 아래   0.0080   0.0019    4.18

초과첨도 (정규분포는 0)
  S&P 500   +1.2193   왜도 -0.2340
  AAPL      +3.0215   왜도 -0.4122

둘 다 -1% 아래로 빠지는 날   →  실제 4.40%,  모형 4.34%   (배율 1.01)
둘 다 -2% 아래로 빠지는 날   →  실제 0.80%,  모형 0.19%   (배율 4.18)

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n = 20_000
x = rng.normal(0, 1, n)
y = 0.7 * x + rng.normal(0, 0.7, n)

fig, ax = plt.subplots()
print(f"{'gridsize':>10}{'칸 수':>10}{'최대 칸':>10}{'평균(빈칸 제외)':>18}")
for g in (10, 30, 60, 120, 200):
    counts = ax.hexbin(x, y, gridsize=g).get_array()
    nz = np.asarray(counts)[np.asarray(counts) > 0]
    print(f"{g:>10}{len(counts):>10}{int(counts.max()):>10}{nz.mean():>18.2f}")
    ax.clear()
plt.close(fig)

  gridsize       칸 수      최대 칸         평균(빈칸 제외)
        10       116      2327            307.69
        30      1068       270             44.44
        60      4175        81             14.37
       120     16750        29              5.03
       200     46316        14              2.59

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(0)
n = 20_000
x = rng.normal(0, 1, n)
y = 0.7 * x + rng.normal(0, 0.7, n)

fig, ax = plt.subplots()
counts = np.asarray(ax.hexbin(x, y, gridsize=40).get_array())
plt.close(fig)
nz = counts[counts > 0]

print(f"칸 개수: 최소 {nz.min():.0f}  중앙값 {np.median(nz):.0f}  최대 {nz.max():.0f}")
print(f"\n중앙값 칸의 상대 밝기")
print(f"  선형 색척도 {np.median(nz) / nz.max():.4f}   ← 거의 배경과 구별되지 않는다")
print(f"  로그 색척도 {np.log(np.median(nz)) / np.log(nz.max()):.4f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
hb1 = axes[0].hexbin(x, y, gridsize=40, cmap="viridis")
fig.colorbar(hb1, ax=axes[0], label="개수")
axes[0].set_title("선형 색척도 — 중심만 보인다", fontsize=10)
hb2 = axes[1].hexbin(x, y, gridsize=40, cmap="viridis", norm=LogNorm())
fig.colorbar(hb2, ax=axes[1], label="개수 (로그)")
axes[1].set_title("로그 색척도 — 꼬리까지 보인다", fontsize=10)
fig.tight_layout()
plt.show()

칸 개수: 최소 1  중앙값 9  최대 164

중앙값 칸의 상대 밝기
  선형 색척도 0.0549   ← 거의 배경과 구별되지 않는다
  로그 색척도 0.4308

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(1)
n = 12_000
x = np.concatenate([rng.normal(-1, 0.6, n // 2), rng.normal(1.5, 0.5, n // 2)])
y = np.concatenate([rng.normal(-1, 0.6, n // 2), rng.normal(1.0, 0.8, n // 2)])

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
axes[0].hexbin(x, y, gridsize=35, cmap="viridis")
axes[0].set_title("hexbin", fontsize=10)
axes[1].hist2d(x, y, bins=35, cmap="viridis")
axes[1].set_title("hist2d (사각 격자)", fontsize=10)

xi, yi = np.mgrid[x.min():x.max():120j, y.min():y.max():120j]
zi = gaussian_kde(np.vstack([x, y]))(np.vstack([xi.ravel(), yi.ravel()]))
axes[2].contourf(xi, yi, zi.reshape(xi.shape), levels=14, cmap="viridis")
axes[2].set_title("2차원 KDE (등고선)", fontsize=10)
fig.tight_layout()
plt.show()

# 계산량 비교: KDE 는 격자점마다 모든 관측을 훑는다
print(f"자료 {len(x):,}개, 평가 격자 {xi.size:,}점")
print(f"  hist2d/hexbin: 관측을 한 번 훑는다        → {len(x):,} 연산 규모")
print(f"  2차원 KDE:     격자점마다 모든 관측을 본다 → {len(x) * xi.size:,} 연산 규모")

자료 12,000개, 평가 격자 14,400점
  hist2d/hexbin: 관측을 한 번 훑는다        → 12,000 연산 규모
  2차원 KDE:     격자점마다 모든 관측을 본다 → 172,800,000 연산 규모

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(2)
n = 6000
# 중심에 밀집하고 외곽으로 갈수록 희박한 자료
r = rng.exponential(1.0, n)
theta = rng.uniform(0, 2 * np.pi, n)
x, y = r * np.cos(theta), r * np.sin(theta)
value = 50 + 5 * r + rng.normal(0, 15, n)      # 참 평균은 반지름에 따라 완만히 증가

fig, ax = plt.subplots()
means = np.asarray(ax.hexbin(x, y, C=value, gridsize=25,
                             reduce_C_function=np.mean, mincnt=1).get_array())
cnt = np.asarray(ax.hexbin(x, y, C=value, gridsize=25,
                           reduce_C_function=len, mincnt=1).get_array())
plt.close(fig)

valid = ~np.isnan(means)
print(f"칸 {valid.sum()}개")
for lo, hi, label in [(1, 1, "점 1개"), (2, 4, "점 2~4개"), (5, 20, "점 5~20개"), (21, 10**9, "21개 이상")]:
    m = valid & (cnt >= lo) & (cnt <= hi)
    if m.sum():
        print(f"  {label:>10} ({m.sum():>3}칸): 평균값의 범위 "
              f"{means[m].min():>7.1f} ~ {means[m].max():>7.1f}   "
              f"표준편차 {means[m].std():>6.2f}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
h1 = axes[0].hexbin(x, y, C=value, gridsize=25, reduce_C_function=np.mean, cmap="viridis")
fig.colorbar(h1, ax=axes[0]); axes[0].set_title("mincnt 없음 — 외곽이 튄다", fontsize=10)
h2 = axes[1].hexbin(x, y, C=value, gridsize=25, reduce_C_function=np.mean,
                    mincnt=8, cmap="viridis")
fig.colorbar(h2, ax=axes[1]); axes[1].set_title("mincnt=8 — 안정된 칸만", fontsize=10)
fig.tight_layout()
plt.show()

칸 251개
        점 1개 ( 90칸): 평균값의 범위    40.3 ~   110.6   표준편차  15.45
      점 2~4개 ( 56칸): 평균값의 범위    46.9 ~   106.3   표준편차  10.42
     점 5~20개 ( 61칸): 평균값의 범위    49.8 ~    79.5   표준편차   5.39
      21개 이상 ( 44칸): 평균값의 범위    50.9 ~    64.9   표준편차   3.48

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(3)
nA, nB = 20_000, 5_000                       # 표본 크기가 다르다
xa, ya = rng.normal(0, 1, nA), rng.normal(0, 1, nA)
xb = rng.normal(0.4, 1, nB)
yb = rng.normal(0.2, 1, nB)

edges = [np.linspace(-4, 4, 31), np.linspace(-4, 4, 31)]
ha, _, _ = np.histogram2d(xa, ya, bins=edges)
hb, _, _ = np.histogram2d(xb, yb, bins=edges)

print(f"A: {nA}개, 최대 칸 {ha.max():.0f}")
print(f"B: {nB}개,  최대 칸 {hb.max():.0f}")
print(f"→ 같은 색척도로 나란히 그리면 B 가 통째로 옅어 보인다\n")

pa, pb = ha / ha.sum(), hb / hb.sum()        # 밀도로 정규화
diff = pa - pb
print(f"밀도 차이의 범위 {diff.min():+.5f} ~ {diff.max():+.5f}")
print(f"차이가 가장 큰 칸의 위치: A 우세 {np.unravel_index(diff.argmax(), diff.shape)}, "
      f"B 우세 {np.unravel_index(diff.argmin(), diff.shape)}")

fig, axes = plt.subplots(1, 3, figsize=(14, 4))
axes[0].hexbin(xa, ya, gridsize=25, extent=(-4, 4, -4, 4), cmap="viridis")
axes[0].set_title(f"A (n={nA})", fontsize=10)
axes[1].hexbin(xb, yb, gridsize=25, extent=(-4, 4, -4, 4), cmap="viridis")
axes[1].set_title(f"B (n={nB}) — 같은 색척도라면 옅다", fontsize=10)
m = np.abs(diff).max()
im = axes[2].imshow(diff.T, origin="lower", extent=(-4, 4, -4, 4),
                    cmap="RdBu_r", vmin=-m, vmax=m)
fig.colorbar(im, ax=axes[2], label="A 밀도 - B 밀도")
axes[2].set_title("밀도 차이 (발산형 색지도)", fontsize=10)
fig.tight_layout()
plt.show()

A: 20000개, 최대 칸 246
B: 5000개,  최대 칸 67
→ 같은 색척도로 나란히 그리면 B 가 통째로 옅어 보인다

밀도 차이의 범위 -0.00495 ~ +0.00455
차이가 가장 큰 칸의 위치: A 우세 (11, 13), B 우세 (18, 14)

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams["font.family"] = "Apple SD Gothic Neo"   # 한글 글꼴
plt.rcParams["axes.unicode_minus"] = False            # 음수 기호가 네모가 되지 않게

rng = np.random.default_rng(4)

small = (rng.normal(0, 1, 60), rng.normal(0, 1, 60))                    # n 이 작다
grid_x = rng.choice([1, 2, 3, 4, 5], 3000)                              # 이산 격자
grid_y = rng.choice([1, 2, 3, 4, 5], 3000)

fig, ax = plt.subplots()
c_small = np.asarray(ax.hexbin(*small, gridsize=20).get_array())
ax.clear()
c_grid = np.asarray(ax.hexbin(grid_x, grid_y, gridsize=20).get_array())
plt.close(fig)

print(f"n=60 자료: 비어 있지 않은 칸 {np.sum(c_small > 0)}개, "
      f"그중 1개짜리 {np.sum(c_small == 1)}개 ({np.mean(c_small[c_small>0] == 1):.1%})")
print(f"이산 격자 자료: 서로 다른 (x,y) 조합 {len(set(zip(grid_x, grid_y)))}개, "
      f"비어 있지 않은 칸 {np.sum(c_grid > 0)}개")
print("  → 25개 격자점이 여러 칸에 쪼개지거나 합쳐진다")

fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
axes[0].hexbin(*small, gridsize=20, cmap="viridis")
axes[0].set_title("n=60 에 hexbin — 대부분 1개짜리 칸", fontsize=10)
axes[1].scatter(*small, s=25, alpha=0.7)
axes[1].set_title("같은 자료의 산점도", fontsize=10)
fig.tight_layout()
plt.show()

n=60 자료: 비어 있지 않은 칸 52개, 그중 1개짜리 45개 (86.5%)
이산 격자 자료: 서로 다른 (x,y) 조합 25개, 비어 있지 않은 칸 25개
  → 25개 격자점이 여러 칸에 쪼개지거나 합쳐진다

import numpy as np

print(f"{'n':>8}{'권장 도구':>28}{'이유':>30}")
rules = [
    (50, "산점도", "모든 점이 보인다"),
    (500, "산점도 + 투명도", "겹침이 시작된다"),
    (5_000, "hexbin 또는 alpha", "겹침이 심하다"),
    (50_000, "hexbin", "개별 점은 무의미"),
    (5_000_000, "hexbin (구간화 필수)", "산점도는 그려지지도 않는다"),
]
for n, tool, why in rules:
    print(f"{n:>8,}{tool:>28}{why:>30}")

       n                       권장 도구                            이유
      50                         산점도                     모든 점이 보인다
     500                   산점도 + 투명도                      겹침이 시작된다
   5,000             hexbin 또는 alpha                       겹침이 심하다
  50,000                      hexbin                     개별 점은 무의미
5,000,000             hexbin (구간화 필수)                산점도는 그려지지도 않는다