바이올린 그림¶

import matplotlib.pyplot as plt
import numpy as np

# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 '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(0)

# --- 자료 1: 봉우리가 둘인 분포 -------------------------------------
# 0 근처 500개와 5 근처 500개를 이어 붙인다.
# 두 무리의 한가운데인 2.5 부근에는 자료가 거의 없다.
data_1 = np.concatenate([np.random.normal(0, 1, 500),
                         np.random.normal(5, 1, 500)])

# --- 자료 2: 봉우리가 하나인 분포 -----------------------------------
# 중심을 자료 1과 같은 2.5 에 맞춘다. 중심만 보아서는
# 두 자료를 구별할 수 없게 만드는 것이 목적이다.
data_2 = np.random.normal(2.5, 2, 1000)

# 요약통계량을 먼저 확인한다. 상자그림이 그리는 것이 바로 이 수들이다.
for name, d in [("data_1 (이봉)", data_1), ("data_2 (단봉)", data_2)]:
    q1, q3 = np.percentile(d, [25, 75])
    print(f"{name}: 평균 {d.mean():.2f}, 중앙값 {np.median(d):.2f}, "
          f"표준편차 {d.std():.2f}, IQR [{q1:.2f}, {q3:.2f}]")

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

# --- 왼쪽: 바이올린 그림 --------------------------------------------
# 좌우로 펼쳐진 폭이 그 높이에서의 커널밀도추정값이다.
# showmeans / showmedians 로 평균선과 중앙값선을 함께 표시한다.
ax1.violinplot([data_1, data_2], showmeans=True, showmedians=True)
ax1.set_title("바이올린 그림")
ax1.set_xticks([1, 2])
ax1.set_xticklabels(["이봉", "단봉"])

# --- 오른쪽: 같은 자료의 상자그림 -----------------------------------
# 다섯 수치 요약과 이상치만 그린다. 밀도 정보는 버려진다.
# 이름표는 boxplot 인자(labels= / tick_labels=)가 버전마다 달라지므로
# 축에 직접 달아 두면 버전에 상관없이 동작한다.
ax2.boxplot([data_1, data_2])
ax2.set_xticks([1, 2])
ax2.set_xticklabels(["이봉", "단봉"])
ax2.set_title("상자그림")

plt.tight_layout()
plt.show()

data_1 (이봉): 평균 2.45, 중앙값 2.40, 표준편차 2.67, IQR [-0.05, 4.93]
data_2 (단봉): 평균 2.53, 중앙값 2.55, 표준편차 1.94, IQR [1.19, 3.75]

from scipy.stats import gaussian_kde

n = len(data_1)
sd = data_1.std(ddof=1)
scott = n ** (-1 / 5)                 # scipy 의 'scott' 어림
h0 = sd * scott                       # 기본 띠폭
print(f"n = {n},  표본 sd = {sd:.4f},  n^(-1/5) = {scott:.4f}")
print(f"기본 띠폭  h0 = sd * n^(-1/5) = {h0:.4f}")

# 참값으로 본 임계 띠폭
mu, sig = 2.5, 1.0
h_crit = np.sqrt(mu ** 2 - sig ** 2)
print(f"참값   mu = {mu}, sigma = {sig} -> 임계 h = {h_crit:.4f}  (배율 {h_crit / h0:.4f})")

# 표본의 값으로 고쳐 본 임계 띠폭
a, b = data_1[:500], data_1[500:]
mu_h = (b.mean() - a.mean()) / 2
sg_h = np.sqrt((a.var(ddof=1) + b.var(ddof=1)) / 2)
h_crit2 = np.sqrt(mu_h ** 2 - sg_h ** 2)
print(f"표본   mu = {mu_h:.4f}, sigma = {sg_h:.4f} -> 임계 h = {h_crit2:.4f}  "
      f"(배율 {h_crit2 / h0:.4f})")


def peaks(factor):
    """띠폭을 기본의 factor 배로 준 KDE 의 봉우리 개수."""
    k = gaussian_kde(data_1, bw_method=factor * scott)
    xs = np.linspace(data_1.min() - 2, data_1.max() + 2, 4000)
    y = k(xs)
    m = (y[1:-1] > y[:-2]) & (y[1:-1] > y[2:])
    return int(m.sum()), xs, y


print("\n배율   띠폭 h   봉우리")
for f in [0.25, 0.5, 1.0, 2.0, 3.0, 3.4, 4.0]:
    print(f"{f:5.2f} {f * h0:8.4f} {peaks(f)[0]:7d}")

lo, hi = 1.0, 8.0                      # 이분법으로 실제 임계 배율을 찾는다
for _ in range(40):
    mid = (lo + hi) / 2
    lo, hi = (mid, hi) if peaks(mid)[0] >= 2 else (lo, mid)
print(f"\n실측 임계 배율 = {lo:.4f}   (h = {lo * h0:.4f})")

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 3.8),
                               gridspec_kw={"width_ratios": [1.4, 1]})
cols = {0.25: "#E65100", 1.0: "#1565C0", 3.4: "#6A1B9A"}
for f, c in cols.items():
    npk, xs, y = peaks(f)
    ax1.plot(xs, y, color=c, lw=1.8,
             label=f"배율 {f} (h = {f * h0:.2f}), 봉우리 {npk}개")
ax1.hist(data_1, bins=40, density=True, color="#CFD8DC", zorder=0)
ax1.set_xlabel("값")
ax1.set_ylabel("밀도")
ax1.set_title("같은 자료, 띠폭만 바꾼 KDE", fontsize=11)
ax1.legend(fontsize=9)
ax1.spines[["top", "right"]].set_visible(False)

for i, (f, c) in enumerate(cols.items(), start=1):
    pp = ax2.violinplot([data_1], positions=[i], widths=0.8,
                        showextrema=False, bw_method=f * scott)
    pp["bodies"][0].set_facecolor(c)
    pp["bodies"][0].set_alpha(0.55)
ax2.set_xticks([1, 2, 3])
ax2.set_xticklabels([f"배율 {f}" for f in cols])
ax2.set_title("같은 자료, 세 가지 바이올린", fontsize=11)
ax2.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
plt.show()

n = 1000,  표본 sd = 2.6706,  n^(-1/5) = 0.2512
기본 띠폭  h0 = sd * n^(-1/5) = 0.6708
참값   mu = 2.5, sigma = 1.0 -> 임계 h = 2.2913  (배율 3.4157)
표본   mu = 2.4801, sigma = 0.9878 -> 임계 h = 2.2749  (배율 3.3912)

배율   띠폭 h   봉우리
 0.25   0.1677       3
 0.50   0.3354       2
 1.00   0.6708       2
 2.00   1.3416       2
 3.00   2.0124       2
 3.40   2.2808       1
 4.00   2.6833       1

실측 임계 배율 = 3.3120   (h = 2.2217)

import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt

# 그림에 한글이 들어가므로 한글 글꼴을 지정한다(지정하지 않으면 네모로 깨진다).
plt.rcParams['font.family'] = 'Apple SD Gothic Neo'
plt.rcParams['axes.unicode_minus'] = False

# 자료를 인터넷에서 내려받으므로 실행에 연결이 필요하다.
url = "https://raw.githubusercontent.com/datasciencedojo/datasets/f0ccab6a7ceafdff780052166fb6fab3311398eb/titanic.csv"
df = pd.read_csv(url)

# 그림에 앞서 숫자로 먼저 확인한다. Age에 결측이 있어 count가 891보다 작다.
print(df.groupby(["Pclass", "Sex"])["Age"]
        .agg(["count", "median", "mean"]).round(1))

fig, ax = plt.subplots(figsize=(10, 4))

# x   : 바이올린을 나눌 기준 (객실 등급 1, 2, 3)
# y   : 분포를 볼 값 (나이)
# hue : 색으로 구분할 두 번째 범주 (성별)
# split=True : 두 색을 하나의 바이올린 좌우에 붙여 그린다.
#              범주가 정확히 둘일 때만 쓸 수 있고, 같은 등급 안에서
#              남녀를 곧바로 견주어 볼 수 있게 해 준다.
# 범례에 한글이 나오도록 값 자체를 한글로 바꾼 열을 따로 만든다.
# ax.legend(labels=[...]) 로 이름만 갈아 끼우면 색과 이름이 어긋날 위험이 있다.
df["성별"] = df["Sex"].map({"male": "남성", "female": "여성"})

sns.violinplot(data=df, x="Pclass", y="Age", hue="성별",
               hue_order=["여성", "남성"], split=True, ax=ax)
ax.set_title("객실 등급과 성별에 따른 나이 분포 (타이타닉)")
ax.set_xlabel("객실 등급")
ax.set_ylabel("나이 (세)")
plt.show()

               count  median  mean
Pclass Sex
1      female     85    35.0  34.6
       male      101    40.0  41.3
2      female     74    28.0  28.7
       male       99    30.0  30.7
3      female    102    21.5  21.8
       male      253    25.0  26.5

import numpy as np
from scipy.stats import gaussian_kde, norm

aged = df.dropna(subset=["Age"])
print(f"{'집단':<10}{'n':>5}{'최소':>7}{'bw':>8}{'예측 아래끝':>12}{'0세 미만 질량':>13}")
for (pc, sx), g in aged.groupby(["Pclass", "성별"]):
    x = g["Age"].values
    bw = np.sqrt(gaussian_kde(x, bw_method="scott").covariance[0, 0])
    lo = x.min() - 2 * bw                      # cut = 2 가 기본값
    print(f"{f'{pc}등급 {sx}':<10}{len(x):>5}{x.min():>7.2f}{bw:>8.4f}{lo:>12.4f}"
          f"{norm.cdf((0 - x) / bw).mean():>13.4f}")

fig, axes = plt.subplots(1, 2, figsize=(12, 4), sharey=True)
for ax, cut in zip(axes, [2, 0]):
    sns.violinplot(data=df, x="Pclass", y="Age", hue="성별",
                   hue_order=["여성", "남성"], split=True, cut=cut, ax=ax)
    ax.axhline(0, color="#D32F2F", ls="--", lw=1.4)
    ax.set_title(f"cut = {cut}", fontsize=11)
    ax.set_xlabel("객실 등급")
    ax.set_ylabel("나이 (세)" if cut == 2 else "")
    ax.spines[["top", "right"]].set_visible(False)
    if cut:
        ax.legend(fontsize=9, loc="upper right")
    else:
        ax.get_legend().remove()

# 그려 놓은 다각형에서 꼭짓점을 직접 꺼내 예측과 맞춰 본다.
fig.canvas.draw()
print("\n왼쪽(cut=2) 바이올린 꼭짓점의 실제 아래끝:")
for coll in axes[0].collections:
    if coll.get_paths():
        print(f"   {coll.get_paths()[0].vertices[:, 1].min():.4f}")
print("오른쪽(cut=0):")
for coll in axes[1].collections:
    if coll.get_paths():
        print(f"   {coll.get_paths()[0].vertices[:, 1].min():.4f}")
plt.tight_layout()
plt.show()

집단            n     최소      bw      예측 아래끝     0세 미만 질량
1등급 남성      101   0.92  6.0152    -11.1104       0.0072
1등급 여성       85   2.00  5.5981     -9.1962       0.0045
2등급 남성       99   0.67  5.9014    -11.1328       0.0332
2등급 여성       74   2.00  5.4428     -8.8856       0.0218
3등급 남성      253   0.42  4.0206     -7.6212       0.0155
3등급 여성      102   0.75  5.0479     -9.3457       0.0480

왼쪽(cut=2) 바이올린 꼭짓점의 실제 아래끝:
   -9.1962
   -11.1104
   -8.8856
   -11.1328
   -9.3457
   -7.6212
오른쪽(cut=0):
   2.0000
   0.9200
   2.0000
   0.6700
   0.7500
   0.4200

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(0)

def n_modes(x):
    g = np.linspace(x.min() - 1, x.max() + 1, 600)
    d = gaussian_kde(x)(g)
    return sum(1 for i in range(1, len(d) - 1) if d[i] > d[i - 1] and d[i] > d[i + 1])

print("완전히 단봉인 N(0,1) 자료에서 KDE 가 봉우리 2개 이상을 보일 확률")
for n in (8, 15, 30, 100):
    print(f"  n={n:>4}: {np.mean([n_modes(rng.normal(0, 1, n)) >= 2 for _ in range(2000)]):.4f}")

samples = [rng.normal(0, 1, 8) for _ in range(4)]
fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
axes[0].violinplot(samples, showmedians=True)
axes[0].set_title("바이올린 (n=8 씩) — 넷이 달라 보인다", fontsize=10)
for i, s in enumerate(samples):
    axes[1].scatter(np.full(len(s), i + 1), s, s=30, alpha=0.8)
axes[1].set_title("원자료 — 같은 모집단에서 8개씩", fontsize=10)
for ax in axes:
    ax.set_xticks([1, 2, 3, 4])
    ax.set_xticklabels([f"표본 {i+1}" for i in range(4)])
fig.tight_layout()
plt.show()

완전히 단봉인 N(0,1) 자료에서 KDE 가 봉우리 2개 이상을 보일 확률
  n=   8: 0.1385
  n=  15: 0.1640
  n=  30: 0.1650
  n= 100: 0.1790

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)
x = np.concatenate([rng.normal(-2, 0.7, 150), rng.normal(2, 0.7, 150)])
grid = np.linspace(-6, 6, 800)

# 주의: gaussian_kde 의 bw_method 에 넣는 수는 대역폭 h 가 아니라 '계수'다.
# 실제 대역폭은 h = 계수 x 자료의 표준편차 이므로, 아래에서 h 를 따로 계산한다.
factors = [0.08, 0.2, "scott", 0.8, 1.2]
sd = x.std(ddof=1)

fig, axes = plt.subplots(1, 5, figsize=(16, 3.6), sharey=True)
for ax, bw in zip(axes, factors):
    kde = gaussian_kde(x, bw_method=bw)
    d = kde(grid)
    ax.fill_betweenx(grid, -d, d, alpha=0.6)
    peaks = sum(1 for i in range(1, len(d) - 1) if d[i] > d[i - 1] and d[i] > d[i + 1])
    ax.set_title(f"계수 {bw}  (h={kde.factor * sd:.2f})\n봉우리 {peaks}개", fontsize=9)
    ax.set_xticks([])
axes[0].set_ylabel("값")
fig.tight_layout()
plt.show()

print(f"자료의 표준편차 {sd:.3f}")
for bw in factors:
    kde = gaussian_kde(x, bw_method=bw)
    d = kde(grid)
    peaks = sum(1 for i in range(1, len(d) - 1) if d[i] > d[i - 1] and d[i] > d[i + 1])
    print(f"계수 {str(bw):>7}  ->  대역폭 h={kde.factor * sd:.3f}: 봉우리 {peaks}개")

자료의 표준편차 2.091
계수    0.08  ->  대역폭 h=0.167: 봉우리 5개
계수     0.2  ->  대역폭 h=0.418: 봉우리 2개
계수   scott  ->  대역폭 h=0.668: 봉우리 2개
계수     0.8  ->  대역폭 h=1.672: 봉우리 2개
계수     1.2  ->  대역폭 h=2.509: 봉우리 1개

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(2)
x = rng.exponential(1.0, 800)                  # 반드시 0 이상인 자료
grid = np.linspace(-1.5, 6, 900)

naive = gaussian_kde(x)(grid)
log_kde = gaussian_kde(np.log(x))
pos = grid > 0
transformed = np.zeros_like(grid)
transformed[pos] = log_kde(np.log(grid[pos])) / grid[pos]

print(f"단순 KDE 가 x<0 에 배정한 질량 {np.trapz(naive[grid < 0], grid[grid < 0]):.4f}")
print(f"x=0.05 에서: 단순 {gaussian_kde(x)(0.05)[0]:.4f}   "
      f"로그변환 {transformed[np.argmin(np.abs(grid - 0.05))]:.4f}   "
      f"참값 {np.exp(-0.05):.4f}")

fig, axes = plt.subplots(1, 3, figsize=(13, 4), sharey=True)
for ax, (d, title) in zip(axes, [
        (naive, "단순 KDE — 0 아래로 새어 나간다"),
        (np.where(grid >= 0, naive, 0), "단순히 잘라 내기 — 편향은 남는다"),
        (transformed, "로그변환 후 되돌리기")]):
    ax.fill_betweenx(grid, -d, d, alpha=0.6)
    ax.axhline(0, color="red", ls="--", lw=1.2)
    ax.set_title(title, fontsize=9)
    ax.set_xticks([])
axes[0].set_ylabel("값")
fig.tight_layout()
plt.show()

단순 KDE 가 x<0 에 배정한 질량 0.0790
x=0.05 에서: 단순 0.5156   로그변환 0.8476   참값 0.9512

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(0)
likert = rng.choice([1, 2, 3, 4, 5], 3000, p=[.05, .15, .40, .30, .10]).astype(float)
grid = np.linspace(-1, 7, 800)
d = gaussian_kde(likert)(grid)

print(f"실제 분포: {np.round([np.mean(likert == k) for k in range(1, 6)], 3)}")
print(f"KDE 가 1 미만에 배정한 질량 {np.trapz(d[grid < 1], grid[grid < 1]):.4f}")
print(f"KDE 가 5 초과에 배정한 질량 {np.trapz(d[grid > 5], grid[grid > 5]):.4f}")
print(f"KDE 봉우리 개수 {sum(1 for i in range(1, len(d) - 1) if d[i] > d[i-1] and d[i] > d[i+1])}")

fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].fill_betweenx(grid, -d, d, alpha=0.6)
axes[0].set_title("바이올린 — 있지도 않은 0.5 나 5.5 를 그린다", fontsize=10)
axes[0].set_xticks([])
levels, counts = np.unique(likert, return_counts=True)
axes[1].bar(levels, counts / counts.sum(), width=0.6)
axes[1].set_title("막대그래프 — 자료 그대로", fontsize=10)
axes[1].set_xlabel("응답")
fig.tight_layout()
plt.show()

실제 분포: [0.053 0.156 0.398 0.295 0.098]
KDE 가 1 미만에 배정한 질량 0.0259
KDE 가 5 초과에 배정한 질량 0.0474
KDE 봉우리 개수 5