히스토그램과 밀도 그림¶

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import scipy.stats as stats
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)      # scipy의 rvs 도 numpy의 전역 난수를 쓴다.
                       # 시드를 고정해야 아래 출력이 재현된다.

samples = 10_000
x = stats.norm(loc=5, scale=10).rvs(samples)     # 평균 5, 표준편차 10의 정규분포

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

# density=True 로 넓이의 합이 1이 되게 정규화한다.
# 이렇게 해야 확률밀도함수와 같은 눈금 위에 놓여 겹쳐 그릴 수 있다.
# hist는 (도수, 구간경계, 막대객체)를 돌려주므로 가운데만 받아 둔다.
_, bins, _ = ax.hist(x, bins=100, density=True, color="#1565C0",
                     label="히스토그램 (구간 100개)")

# 표본에서 추정한 모수로 정규 밀도함수를 만든다.
# 참값(5, 10)이 아니라 표본에서 잰 값을 쓴다는 점이 중요하다.
# 실제 분석에서는 참값을 모르기 때문이다.
x_mean = x.mean()
x_std = x.std(ddof=1)
pdf = stats.norm(loc=x_mean, scale=x_std).pdf(bins)

# 적합된 밀도곡선을 겹쳐 그린다.
# 한글은 $...$ 바깥에 둔다. 수식 글꼴에는 한글 글리프가 없다.
ax.plot(bins, pdf, "-", color="#D32F2F", linewidth=2,
        label="적합된 정규 밀도")
ax.set_xlabel("관측값")
ax.set_ylabel("밀도")
ax.set_title("정규 표본 10,000개의 히스토그램과 적합된 밀도곡선")
ax.legend()
ax.spines[["top", "right"]].set_visible(False)
fig.savefig("histograms_17.png", dpi=170, facecolor="white",
            bbox_inches="tight")

print(f"표본평균   {x_mean:.3f}  (참값 5)")
print(f"표본표준편차 {x_std:.3f}  (참값 10)")

# --- 두 추정값이 참값과 어긋나는가. 표준오차로 잰다 ---
se_mean = 10 / np.sqrt(samples)
se_sd = 10 / np.sqrt(2 * samples)
print(f"  평균:   SE = {se_mean:.4f},  z = {(x_mean - 5) / se_mean:+.3f}")
print(f"  표준편차: SE = {se_sd:.4f},  z = {(x_std - 10) / se_sd:+.3f}")

# --- 막대 넓이의 합이 1 인가 (높이의 합이 아니다) ---
dens, edges = np.histogram(x, bins=100, density=True)
cnt, _ = np.histogram(x, bins=100)
w = edges[1] - edges[0]
print(f"\n구간 너비 w = {w:.6f},  넓이의 합 = {(dens * w).sum():.12f},  높이의 합 = {dens.sum():.4f}")

# --- 평균이 든 구간에서 높이의 기댓값과 흔들림 ---
F = stats.norm(5, 10)
j = int(np.searchsorted(edges, 5)) - 1
p = F.cdf(edges[j + 1]) - F.cdf(edges[j])
c = (edges[j] + edges[j + 1]) / 2
print(f"평균이 든 구간 [{edges[j]:.4f}, {edges[j + 1]:.4f}),  가운데 c = {c:.4f}")
print(f"  p = {p:.6f},  E[도수] = {samples * p:.2f},  sd(도수) = {np.sqrt(samples * p * (1 - p)):.2f},"
      f"  실제 도수 = {cnt[j]}")
print(f"  E[높이] = p/w = {p / w:.8f},  실제 높이 = {dens[j]:.8f}")
print(f"  f(c)    =       {F.pdf(c):.8f}")
# 테일러 전개: (1/w)∫f = f(c) + (w^2/24) f''(c) + ...
fpp = F.pdf(c) * ((c - 5) ** 2 / 10 ** 4 - 1 / 10 ** 2)
print(f"  차이 E[높이] - f(c) = {p / w - F.pdf(c):.4e},  예측 (w^2/24)f''(c) = {w ** 2 / 24 * fpp:.4e}")
print(f"  구간당 상대 잡음 = {np.sqrt((1 - p) / (samples * p)):.4f},"
      f"  상대 치우침 = {abs(p / w - F.pdf(c)) / F.pdf(c):.6f}")
print(f"  가장 높은 막대 = {dens.max():.6f}  (참 봉우리 1/(10*sqrt(2pi)) = {1 / (10 * np.sqrt(2 * np.pi)):.6f})")

# --- 구간 100 개는 적절한가 ---
iqr = np.subtract(*np.percentile(x, [75, 25]))
rng_ = x.max() - x.min()
print(f"\n스콧 h = {3.49 * x_std * samples ** (-1 / 3):.4f} -> 구간 {rng_ / (3.49 * x_std * samples ** (-1 / 3)):.1f}개")
print(f"FD   h = {2 * iqr * samples ** (-1 / 3):.4f} -> 구간 {rng_ / (2 * iqr * samples ** (-1 / 3)):.1f}개")
print(f"스터지스 k = 1 + log2(n) = {1 + np.log2(samples):.1f}개")

표본평균   4.816  (참값 5)
표본표준편차 9.876  (참값 10)
  평균:   SE = 0.1000,  z = -1.843
  표준편차: SE = 0.0707,  z = -1.753

구간 너비 w = 0.754176,  넓이의 합 = 1.000000000000,  높이의 합 = 1.3260
평균이 든 구간 [4.5536, 5.3078),  가운데 c = 4.9307
  p = 0.030079,  E[도수] = 300.79,  sd(도수) = 17.08,  실제 도수 = 330
  E[높이] = p/w = 0.03988382,  실제 높이 = 0.04375636
  f(c)    =       0.03989327
  차이 E[높이] - f(c) = -9.4519e-06,  예측 (w^2/24)f''(c) = -9.4539e-06
  구간당 상대 잡음 = 0.0568,  상대 치우침 = 0.000237
  가장 높은 막대 = 0.045347  (참 봉우리 1/(10*sqrt(2pi)) = 0.039894)

스콧 h = 1.5998 -> 구간 47.1개
FD   h = 1.2400 -> 구간 60.8개
스터지스 k = 1 + log2(n) = 14.3개

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

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

def plot_loan_income_distribution():
    """대출 신청자 소득의 히스토그램에 정규분포를 겹쳐 그린다.

    앞 보기와 코드 구조는 같지만 결론이 정반대다.
    앞에서는 곡선이 히스토그램에 잘 맞았고, 여기서는 맞지 않는다.
    """
    url = 'https://raw.githubusercontent.com/gedeck/practical-statistics-for-data-scientists/8a6d3bb6468e979c861d4b37215e1413702dfdfa/data/loans_income.csv'
    df = pd.read_csv(url)

    # 소득의 평균과 표준편차. 이 둘만으로 정규분포가 결정된다.
    mean_income = df['x'].mean()
    std_dev_income = df['x'].std()

    fig, ax = plt.subplots(figsize=(15, 4))
    _, bins, _ = ax.hist(df['x'], bins=30, density=True,
                         color='#DCEBFB', edgecolor='#1565C0',
                         label='소득 히스토그램')

    # 같은 평균·표준편차를 갖는 정규분포를 겹쳐 그린다.
    # 두 곡선이 어긋나는 방식이 곧 "자료가 정규분포와 어떻게 다른가"를 말해 준다.
    norm_pdf = stats.norm(loc=mean_income, scale=std_dev_income).pdf(bins)
    ax.plot(bins, norm_pdf, "--", color="#D32F2F", lw=2,
            label='같은 평균·표준편차의 정규분포')

    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.set_title('대출 신청자 소득의 분포와 정규 적합')
    ax.set_xlabel('소득 (달러)')
    ax.set_ylabel('밀도')
    ax.legend()
    fig.savefig("histograms_46.png", dpi=170, facecolor="white",
                bbox_inches="tight")

    # 치우침을 숫자로 확인한다.
    # 오른쪽으로 치우치면 평균이 중앙값보다 크고 왜도가 양수다.
    print(f"평균   {mean_income:,.0f}")
    print(f"중앙값 {df['x'].median():,.0f}")
    print(f"왜도   {stats.skew(df['x']):.3f}  (0이면 대칭)")

    # --- 적합된 정규분포가 음수 소득에 주는 확률 ---
    n = len(df)
    z0 = (0 - mean_income) / std_dev_income
    p_neg = stats.norm(loc=mean_income, scale=std_dev_income).cdf(0)
    print(f"\nn = {n:,},  표준편차 {std_dev_income:,.0f}")
    print(f"정규 적합이 소득 < 0 에 주는 확률 = Phi({z0:.4f}) = {p_neg:.6f}"
          f"  -> {n * p_neg:,.0f} 명")
    print(f"  자료의 실제 최솟값 = {df['x'].min():,}  (음수인 사람 {int((df['x'] < 0).sum())} 명)")

    # --- 피어슨의 둘째 왜도계수도 같은 방향을 가리킨다 ---
    pearson2 = 3 * (mean_income - df['x'].median()) / std_dev_income
    print(f"  피어슨 둘째 왜도계수 3(mean-median)/s = {pearson2:.4f}")

    # --- 구간 개수: 스터지스는 이 자료에서 얼마나 모자라는가 ---
    iqr = np.subtract(*np.percentile(df['x'], [75, 25]))
    span = df['x'].max() - df['x'].min()
    sturges = 1 + np.log2(n)
    fd_h = 2 * iqr * n ** (-1 / 3)
    scott_h = 3.49 * std_dev_income * n ** (-1 / 3)
    print(f"\nIQR = {iqr:,.0f},  범위 = {span:,.0f},  n^(1/3) = {n ** (1 / 3):.1f}")
    print(f"  스터지스 k = 1 + log2(n) = {sturges:.1f} 개  (너비 {span / sturges:,.0f})")
    print(f"  FD       h = {fd_h:,.0f}  ->  {span / fd_h:.1f} 개")
    print(f"  스콧     h = {scott_h:,.0f}  ->  {span / scott_h:.1f} 개")
    print(f"  numpy bins='auto' 가 고른 구간 수 = {len(np.histogram_bin_edges(df['x'], bins='auto')) - 1}")

if __name__ == "__main__":
    plot_loan_income_distribution()

평균   68,761
중앙값 62,000
왜도   1.049  (0이면 대칭)

n = 50,000,  표준편차 32,872
정규 적합이 소득 < 0 에 주는 확률 = Phi(-2.0918) = 0.018230  -> 911 명
  자료의 실제 최솟값 = 4,000  (음수인 사람 0 명)
  피어슨 둘째 왜도계수 3(mean-median)/s = 0.6170

IQR = 40,000,  범위 = 195,000,  n^(1/3) = 36.8
  스터지스 k = 1 + log2(n) = 16.6 개  (너비 11,740)
  FD       h = 2,172  ->  89.8 개
  스콧     h = 3,114  ->  62.6 개
  numpy bins='auto' 가 고른 구간 수 = 90

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import os
import pandas as pd
import tarfile
import urllib.request

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

# 캘리포니아 주택 자료. 구역마다 소득·집값·방 수 등 아홉 개 변수가 들어 있다.
DOWNLOAD_ROOT = "https://raw.githubusercontent.com/ageron/handson-ml2/7b7e23e7267356f8355580877eff98c43cda1bd0/"
HOUSING_PATH = os.path.join("datasets", "housing")
HOUSING_URL = DOWNLOAD_ROOT + "datasets/housing/housing.tgz"

def fetch_housing_data(housing_url=HOUSING_URL, housing_path=HOUSING_PATH):
    """압축 파일을 내려받아 풀어 둔다. 이미 받아 두었으면 다시 받지 않는다."""
    if not os.path.isdir(housing_path):
        os.makedirs(housing_path)
    tgz_path = os.path.join(housing_path, "housing.tgz")
    urllib.request.urlretrieve(housing_url, tgz_path)
    with tarfile.open(tgz_path) as housing_tgz:
        housing_tgz.extractall(path=housing_path)

def load_housing_data(housing_path=HOUSING_PATH):
    """풀어 둔 csv 를 자료틀로 읽는다."""
    csv_path = os.path.join(housing_path, "housing.csv")
    return pd.read_csv(csv_path)

fetch_housing_data()
df = load_housing_data()

# 3×3 격자에 아홉 변수를 한꺼번에 그린다. 자료를 처음 만났을 때
# 어느 변수가 치우쳤는지, 어디가 잘렸는지 한눈에 훑는 방법이다.
fig, axes = plt.subplots(3, 3, figsize=(12, 9))
df.hist(bins=50, ax=axes, color="#1565C0")

# 열 이름은 영어이므로 그림에 쓸 한글 이름을 따로 준비한다.
KOREAN = {
    "longitude": "경도", "latitude": "위도",
    "housing_median_age": "주택 연식 중앙값 (년)",
    "total_rooms": "구역의 방 수", "total_bedrooms": "구역의 침실 수",
    "population": "구역 인구", "households": "구역 가구 수",
    "median_income": "소득 중앙값 (만 달러)",
    "median_house_value": "집값 중앙값 (달러)",
}

# 격자를 1차원으로 펴서 아홉 축을 차례로 다듬는다.
for ax in axes.reshape((-1,)):
    ax.set_title(KOREAN.get(ax.get_title(), ax.get_title()), fontsize=10)
    ax.set_ylabel("구역 수", fontsize=9)
    ax.grid(False)
    ax.spines[["top", "right"]].set_visible(False)

fig.suptitle("캘리포니아 주택 자료 아홉 변수의 히스토그램")
fig.tight_layout()
fig.savefig("histograms_83.png", dpi=170, facecolor="white",
            bbox_inches="tight")

# --- 눈으로 읽은 것을 수로 확인한다 ---
# 그림에서 오른쪽 끝에 솟은 막대는 "잘린 자료"의 흔적이다.
# 몇 개가 거기 쌓여 있는지는 그림이 말해 주지 않으므로 직접 센다.
from scipy import stats

num = df.select_dtypes("number")
print(f"구역 {len(df):,} 개,  수치형 변수 {num.shape[1]} 개")
for c in num.columns:
    s = num[c]
    at_max = int((s == s.max()).sum())
    print(f"  {KOREAN[c]}")
    print(f"      결측 {int(s.isna().sum()):>4},  왜도 {stats.skew(s.dropna()):>7.3f},"
          f"  최댓값 {s.max():>12,.4f} 에 쌓인 구역 {at_max:>5,d} ({100 * at_max / len(s):5.2f}%)")

구역 20,640 개,  수치형 변수 9 개
  경도
      결측    0,  왜도  -0.298,  최댓값    -114.3100 에 쌓인 구역     1 ( 0.00%)
  위도
      결측    0,  왜도   0.466,  최댓값      41.9500 에 쌓인 구역     2 ( 0.01%)
  주택 연식 중앙값 (년)
      결측    0,  왜도   0.060,  최댓값      52.0000 에 쌓인 구역 1,273 ( 6.17%)
  구역의 방 수
      결측    0,  왜도   4.147,  최댓값  39,320.0000 에 쌓인 구역     1 ( 0.00%)
  구역의 침실 수
      결측  207,  왜도   3.459,  최댓값   6,445.0000 에 쌓인 구역     1 ( 0.00%)
  구역 인구
      결측    0,  왜도   4.935,  최댓값  35,682.0000 에 쌓인 구역     1 ( 0.00%)
  구역 가구 수
      결측    0,  왜도   3.410,  최댓값   6,082.0000 에 쌓인 구역     1 ( 0.00%)
  소득 중앙값 (만 달러)
      결측    0,  왜도   1.647,  최댓값      15.0001 에 쌓인 구역    49 ( 0.24%)
  집값 중앙값 (달러)
      결측    0,  왜도   0.978,  최댓값 500,001.0000 에 쌓인 구역   965 ( 4.68%)

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

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, index_col='PassengerId')

# 다섯 변수를 한 줄에 나란히 그린다.
# 자료를 처음 받았을 때 모든 변수를 한눈에 훑는 표준적인 방법이다.
fig, axes = plt.subplots(1, 5, figsize=(13, 3))
columns = ("Sex", "Survived", "Age", "Pclass", "Age")
titles = ("성별 (명목형)", "생존 여부 (이진형)", "나이 (연속형)",
          "객실등급 (순서형)", "나이 (구간 20개)")
bins = (10, 10, 10, 10, 20)

for ax, col, title, b in zip(axes, columns, titles, bins):
    # 변수 유형이 섞여 있다는 점에 주목하라.
    #   Sex      문자열 범주형 -> 히스토그램이 사실상 막대그림이 된다
    #   Survived 0/1 이진형    -> 막대 두 개
    #   Pclass   1/2/3 순서형  -> 막대 세 개
    #   Age      연속형        -> 진짜 히스토그램
    # 범주형에 히스토그램을 쓰는 것은 원칙적으로 맞지 않지만,
    # 탐색 단계에서 빠르게 훑을 때는 흔히 이렇게 한다.
    ax.hist(df[col], bins=b, density=True, color="#DCEBFB",
            edgecolor="#1565C0")
    ax.set_title(title, fontsize=10)
    ax.set_ylabel("밀도", fontsize=9)
    ax.spines[["top", "right"]].set_visible(False)

fig.suptitle("타이타닉 자료: 종류가 섞인 다섯 변수를 히스토그램으로 훑기")
fig.tight_layout()
fig.savefig("histograms_115.png", dpi=170, facecolor="white",
            bbox_inches="tight")

print(df[["Sex", "Survived", "Age", "Pclass"]].dtypes)

# --- 이진 변수에 density=True 를 쓰면 높이가 무엇이 되는가 ---
# 0/1 자료를 [0,1] 위에서 k 등분하면 0 은 첫 구간, 1 은 마지막 구간에 모인다.
# 그 높이는 비율/너비 = k x 비율이므로 k 를 바꾸면 그대로 따라 변한다.
p0 = (df.Survived == 0).mean()
p1 = (df.Survived == 1).mean()
print(f"\n사망 비율 {p0:.6f},  생존 비율 {p1:.6f}")
for k in (10, 20):
    h, e = np.histogram(df.Survived, bins=k, density=True)
    nz = [round(v, 4) for v in h if v > 0]
    print(f"  구간 {k:2d}개 (너비 {e[1] - e[0]:.4f}): 0 이 아닌 높이 {nz}"
          f"   예측 {k * p0:.4f}, {k * p1:.4f}")

# --- 나이: 몇 명이 그림에서 빠졌는가 ---
age = df.Age
print(f"\n전체 {len(df)} 명 중 나이 결측 {int(age.isna().sum())} 명"
      f" ({100 * age.isna().mean():.2f}%) — 히스토그램은 말없이 빼고 그린다")
a = age.dropna()
print(f"  그려진 자료 {len(a)} 명,  범위 [{a.min()}, {a.max()}],  "
      f"평균 {a.mean():.4f},  중앙값 {a.median()}")

# --- 구간 수를 바꾸면 봉우리가 몇 개로 세어지는가 ---
def count_modes(x, k):
    """구간 k 개로 끊었을 때 지역 최대인 구간의 가운데 값을 돌려준다."""
    h, e = np.histogram(x, bins=k)
    mid = (e[:-1] + e[1:]) / 2
    out = []
    for i in range(k):
        left = h[i - 1] if i > 0 else -1
        right = h[i + 1] if i < k - 1 else -1
        if h[i] > left and h[i] > right:
            out.append(round(mid[i], 2))
    return out

for k in (10, 20):
    peaks = count_modes(a, k)
    print(f"  구간 {k:2d}개: 봉우리 {len(peaks)} 개 — 나이 {peaks}")

Sex          object
Survived      int64
Age         float64
Pclass        int64
dtype: object

사망 비율 0.616162,  생존 비율 0.383838
  구간 10개 (너비 0.1000): 0 이 아닌 높이 [6.1616, 3.8384]   예측 6.1616, 3.8384
  구간 20개 (너비 0.0500): 0 이 아닌 높이 [12.3232, 7.6768]   예측 12.3232, 7.6768

전체 891 명 중 나이 결측 177 명 (19.87%) — 히스토그램은 말없이 빼고 그린다
  그려진 자료 714 명,  범위 [0.42, 80.0],  평균 29.6991,  중앙값 28.0
  구간 10개: 봉우리 2 개 — 나이 [4.4, 20.32]
  구간 20개: 봉우리 3 개 — 나이 [2.41, 22.3, 70.05]

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

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

def compute_bins_widths_heights():
    """계급의 폭이 제각각인 도수분포표에서 막대의 높이를 구한다.

    폭이 다르면 도수를 그대로 높이로 쓸 수 없다. 넓이가 비율을 나타내야
    하므로 높이는 비율을 폭으로 나눈 값, 곧 밀도가 된다.
    """
    bins = [0, 1_000, 2_000, 3_000, 4_000, 5_000,
            6_000, 7_000, 10_000, 15_000, 25_000, 50_000]
    widths = [right - left for left, right in zip(bins[:-1], bins[1:])]
    percents = [1, 2, 3, 4, 5, 5, 5, 15, 26, 26, 8]
    heights = [p / w for w, p in zip(widths, percents)]
    return bins, widths, heights

def draw_line(start, end, ax):
    """두 점을 잇는 검은 선분 하나."""
    ax.plot([start[0], end[0]], [start[1], end[1]], '-k')

def draw_box(x_left, x_right, height, ax):
    """막대 하나를 네 선분으로 직접 그린다."""
    draw_line([x_left, 0], [x_right, 0], ax)
    draw_line([x_right, 0], [x_right, height], ax)
    draw_line([x_right, height], [x_left, height], ax)
    draw_line([x_left, height], [x_left, 0], ax)

def main():
    """계급마다 폭과 높이가 다른 막대를 이어 붙여 히스토그램을 만든다."""
    bins, widths, heights = compute_bins_widths_heights()
    fig, ax = plt.subplots(figsize=(12, 3))
    for x_left, x_right, height in zip(bins[:-1], bins[1:], heights):
        draw_box(x_left, x_right, height, ax)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['bottom'].set_position("zero")
    ax.set_xlabel("소득 수준 (달러)")
    ax.set_ylabel("밀도 (달러당 백분율)")
    ax.set_title("폭이 다른 계급의 밀도 히스토그램 — 넓이가 백분율이다")
    fig.savefig("histograms_164.png", dpi=170, facecolor="white",
                bbox_inches="tight")

    # 폭·백분율·높이·넓이를 한 표로 적어 본다.
    # 넓이 열의 합이 100 이 되는지가 이 보기의 전부다.
    percents = [1, 2, 3, 4, 5, 5, 5, 15, 26, 26, 8]
    print(f"{'계급':>18}{'폭':>8}{'백분율':>7}{'높이(1000달러당 %)':>20}{'넓이':>7}")
    for lo, hi, w, p, h in zip(bins[:-1], bins[1:], widths, percents, heights):
        print(f"{lo:>7,}-{hi:<9,}{w:>8,}{p:>6}{1000 * h:>14.4f}{w * h:>12.1f}")
    print(f"{'합계':>18}{'':>8}{sum(percents):>6}{'':>14}"
          f"{sum(w * h for w, h in zip(widths, heights)):>12.1f}")
    print(f"\n가장 높은 막대 = {1000 * max(heights):.2f}  "
          f"({bins[heights.index(max(heights))]:,}"
          f"-{bins[heights.index(max(heights)) + 1]:,} 달러)")
    print(f"백분율이 가장 큰 계급 = {max(percents)}%  (두 곳)")

if __name__ == "__main__":
    main()

                계급       폭    백분율       높이(1000달러당 %)     넓이
      0-1,000       1,000     1        1.0000         1.0
  1,000-2,000       1,000     2        2.0000         2.0
  2,000-3,000       1,000     3        3.0000         3.0
  3,000-4,000       1,000     4        4.0000         4.0
  4,000-5,000       1,000     5        5.0000         5.0
  5,000-6,000       1,000     5        5.0000         5.0
  6,000-7,000       1,000     5        5.0000         5.0
  7,000-10,000      3,000    15        5.0000        15.0
 10,000-15,000      5,000    26        5.2000        26.0
 15,000-25,000     10,000    26        2.6000        26.0
 25,000-50,000     25,000     8        0.3200         8.0
                합계           100                     100.0

가장 높은 막대 = 5.20  (10,000-15,000 달러)
백분율이 가장 큰 계급 = 26%  (두 곳)

import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(0)
x = np.concatenate([rng.normal(-1.1, 0.55, 150), rng.normal(1.1, 0.55, 150)])
peaks = lambda h: sum(1 for i in range(1, len(h) - 1)
                      if h[i] > h[i - 1] and h[i] > h[i + 1])

w = 1.3
print(f"같은 자료 {len(x)}개, 같은 구간 너비 {w}, 시작점만 다르게")
for off in (0.0, 0.2, 0.4, 0.6, 0.8):
    edges = np.arange(x.min() - w + off * w, x.max() + w, w)
    h, _ = np.histogram(x, bins=edges)
    print(f"  오프셋 {off:.1f}: 도수 {h}  → 봉우리 {peaks(h)}개")

grid = np.linspace(x.min() - 1, x.max() + 1, 1000)
print("\nKDE 에는 '시작점'이라는 개념이 없다")
for bw in ("scott", 0.2, 0.35, 0.6):
    print(f"  대역폭 {str(bw):>8}: 봉우리 {peaks(gaussian_kde(x, bw_method=bw)(grid))}개")

같은 자료 300개, 같은 구간 너비 1.3, 시작점만 다르게
  오프셋 0.0: 도수 [  0  68  89 110  33]  → 봉우리 1개
  오프셋 0.2: 도수 [  5  95  75 109  16]  → 봉우리 2개
  오프셋 0.4: 도수 [  7 114  73  99   7]  → 봉우리 2개
  오프셋 0.6: 도수 [ 24 116  87  68   5]  → 봉우리 1개
  오프셋 0.8: 도수 [ 47 104 101  46   2]  → 봉우리 1개

KDE 에는 '시작점'이라는 개념이 없다
  대역폭    scott: 봉우리 2개
  대역폭      0.2: 봉우리 2개
  대역폭     0.35: 봉우리 2개
  대역폭      0.6: 봉우리 2개

import numpy as np
from scipy.stats import gaussian_kde

rng = np.random.default_rng(2)
x = rng.exponential(1.0, 5000)              # 반드시 0 이상인 자료
grid = np.linspace(-1, 5, 1200)
kde = gaussian_kde(x)

print(f"x < 0 구간에 배정된 확률질량 {np.trapz(kde(grid[grid < 0]), grid[grid < 0]):.4f}")
print("  ← 있을 수 없는 영역이다\n")

log_kde = gaussian_kde(np.log(x))           # 로그 변환 후 추정하고 되돌린다
transformed = lambda t: log_kde(np.log(t)) / t

print(f"{'x':>6}{'단순 KDE':>11}{'로그변환 KDE':>15}{'참값':>10}")
for t in (0.05, 0.5, 1.0, 2.0):
    print(f"{t:>6.2f}{kde(t)[0]:>11.4f}{transformed(np.array([t]))[0]:>15.4f}"
          f"{np.exp(-t):>10.4f}")

x < 0 구간에 배정된 확률질량 0.0595
  ← 있을 수 없는 영역이다

     x     단순 KDE       로그변환 KDE        참값
  0.05     0.5157         0.9805    0.9512
  0.50     0.6254         0.6110    0.6065
  1.00     0.3802         0.3622    0.3679
  2.00     0.1406         0.1310    0.1353

import numpy as np

rng = np.random.default_rng(0)
print(f"{'n':>8}{'FD':>10}{'스콧':>10}{'스터지스':>11}{'n^(-1/3)':>11}")
for n in (100, 1000, 10_000, 100_000):
    d = rng.normal(0, 1, n)
    iqr = np.subtract(*np.percentile(d, [75, 25]))
    fd = 2 * iqr / n ** (1 / 3)
    scott = 3.49 * d.std(ddof=1) / n ** (1 / 3)
    sturges = (d.max() - d.min()) / (1 + np.log2(n))
    print(f"{n:>8}{fd:>10.4f}{scott:>10.4f}{sturges:>11.4f}{n ** (-1 / 3):>11.4f}")

       n        FD        스콧       스터지스   n^(-1/3)
     100    0.5842    0.7271     0.5661     0.2154
    1000    0.2574    0.3447     0.6352     0.1000
   10000    0.1264    0.1620     0.5147     0.0464
  100000    0.0581    0.0752     0.5239     0.0215

import numpy as np

rng = np.random.default_rng(2)
n = 20_000
true = rng.normal(170, 8, n)
# 60% 는 5 단위로, 40% 는 1 단위로 반올림해 보고한다
reported = np.where(rng.random(n) < 0.6, np.round(true / 5) * 5, np.round(true))

last = reported.astype(int) % 10
print("보고된 키의 끝자리 분포 (균등하다면 각 0.1)")
for d in range(10):
    print(f"  끝자리 {d}: {np.mean(last == d):.4f}")

print(f"\n평균:     참 {true.mean():.4f}   보고 {reported.mean():.4f}")
print(f"표준편차: 참 {true.std():.4f}   보고 {reported.std():.4f}")

보고된 키의 끝자리 분포 (균등하다면 각 0.1)
  끝자리 0: 0.3396
  끝자리 1: 0.0400
  끝자리 2: 0.0391
  끝자리 3: 0.0390
  끝자리 4: 0.0407
  끝자리 5: 0.3390
  끝자리 6: 0.0402
  끝자리 7: 0.0429
  끝자리 8: 0.0398
  끝자리 9: 0.0398

평균:     참 170.0558   보고 170.0445
표준편차: 참 7.9885   보고 8.0645