표본조사¶

"""단순무작위추출과 층화추출을 인공 모집단에서 비교한다."""

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)

# === 두 층으로 이루어진 모집단 ===
# 젊은 층 70%, 나이 든 층 30%. 두 층의 소득 분포가 뚜렷이 다르다.
#   Young: 평균 4만, 표준편차 1만
#   Old:   평균 7만, 표준편차 1.5만
# 층 안에서는 비교적 고르고 층 사이에서 크게 갈리는 이 구조가
# 층화추출이 이득을 보는 전형적인 상황이다.
n_pop = 10_000
stratum = rng.choice(["Young", "Old"], size=n_pop, p=[0.7, 0.3])
income = np.where(
    stratum == "Young",
    rng.normal(40_000, 10_000, n_pop),
    rng.normal(70_000, 15_000, n_pop),
)
pop = pd.DataFrame({"stratum": stratum, "income": income})
true_mean = pop["income"].mean()      # 추정하려는 참값

# === 방법 1: 단순무작위추출(SRS) 200명 ===
# 층을 무시하고 1만 명 중 200명을 그냥 뽑는다.
# 표본에 우연히 노년층이 많이(또는 적게) 들어올 수 있고, 그 우연이 추정값을 흔든다.
srs = pop.sample(200, random_state=1)

# === 방법 2: 비례 층화추출 200명 ===
# 각 층에서 그 층의 모집단 비율만큼 뽑는다: Young 140명, Old 60명.
# 층의 구성비를 표본에서 고정해 버리므로 "우연히 치우칠" 여지가 사라진다.
strat = pop.groupby("stratum", group_keys=False).apply(
    lambda x: x.sample(int(round(200 * len(x) / n_pop)), random_state=1)
)

print(f"True mean:        ${true_mean:,.0f}")
print(f"SRS estimate:     ${srs['income'].mean():,.0f}")
print(f"Stratified est.:  ${strat['income'].mean():,.0f}")

# --- 되풀이해 두 방법의 흔들림을 잰다 ---
w = np.array([0.3, 0.7])                 # Old, Young 의 모집단 비율
sd = np.array([15_000.0, 10_000.0])
mu = np.array([70_000.0, 40_000.0])
mu_bar = (w * mu).sum()
within = (w * sd ** 2).sum()
between = (w * (mu - mu_bar) ** 2).sum()
print(f"이론: 층내분산 {within:,.0f} + 층간분산 {between:,.0f} = 전체분산 {within + between:,.0f}")
print(f"이론 SE: SRS {np.sqrt((within + between) / 200):,.1f}, "
      f"층화 {np.sqrt(within / 200):,.1f},  설계효과 {within / (within + between):.4f}")

idx_y = np.flatnonzero(pop["stratum"].values == "Young")
idx_o = np.flatnonzero(pop["stratum"].values == "Old")
inc = pop["income"].values
reps = 2000
srs_m, str_m = np.empty(reps), np.empty(reps)
for r in range(reps):
    srs_m[r] = inc[rng.choice(n_pop, 200, replace=False)].mean()
    str_m[r] = 0.7 * inc[rng.choice(idx_y, 140, replace=False)].mean() \
             + 0.3 * inc[rng.choice(idx_o, 60, replace=False)].mean()
print(f"모의 {reps}회 SE: SRS {srs_m.std(ddof=1):,.1f}, 층화 {str_m.std(ddof=1):,.1f}, "
      f"비 {str_m.var(ddof=1) / srs_m.var(ddof=1):.4f}")

True mean:        $49,092
SRS estimate:     $49,455
Stratified est.:  $48,375
이론: 층내분산 137,500,000 + 층간분산 189,000,000 = 전체분산 326,500,000
이론 SE: SRS 1,277.7, 층화 829.2,  설계효과 0.4211
모의 2000회 SE: SRS 1,283.4, 층화 820.1, 비 0.4084

import numpy as np

rng = np.random.default_rng(0)
N_h = np.array([5000, 3000, 2000])
mu_h = np.array([40., 55., 90.])
sd_h = np.array([5., 10., 40.])          # 분산이 층마다 크게 다르다

pop = [rng.normal(m, s, k) for m, s, k in zip(mu_h, sd_h, N_h)]
N, n = N_h.sum(), 300
W = N_h / N

prop = np.array([150, 90, 60])                          # 비례:  W_h * n
ney = np.round(n * W * sd_h / (W * sd_h).sum()).astype(int)
ney[-1] = n - ney[:-1].sum()                            # 네이만: W_h sigma_h
print(f"비례 배분:   {prop}")
print(f"네이만 배분: {ney}")

var_of = lambda a: sum(W[i] ** 2 * sd_h[i] ** 2 / a[i] for i in range(3))
print(f"\n이론 분산  비례 {var_of(prop):.5f}   네이만 {var_of(ney):.5f}"
      f"   비 {var_of(prop) / var_of(ney):.3f}")

B = 20_000
sim = lambda a: np.array([sum(W[i] * rng.choice(pop[i], a[i], replace=False).mean()
                              for i in range(3)) for _ in range(B)])
a, b = sim(prop), sim(ney)
truth = sum(W[i] * pop[i].mean() for i in range(3))
print(f"\n참값 {truth:.4f}")
print(f"비례:   평균 {a.mean():.4f}   표준오차 {a.std(ddof=1):.4f}")
print(f"네이만: 평균 {b.mean():.4f}   표준오차 {b.std(ddof=1):.4f}")
print(f"모의 분산비 {a.var() / b.var():.3f}")

비례 배분:   [150  90  60]
네이만 배분: [ 56  67 177]

이론 분산  비례 1.20833   네이만 0.60752   비 1.989

참값 54.7497
비례:   평균 54.7432   표준오차 1.0633
네이만: 평균 54.7464   표준오차 0.7492
모의 분산비 2.014

import numpy as np

rng = np.random.default_rng(1)
N = 7000
day = np.arange(N) % 7
y = 50 + np.where(day == 5, 40., 0.) + rng.normal(0, 5, N)   # 토요일에 매출 급증

print(f"모평균 {y.mean():.4f}")

sys7 = [y[np.arange(s, N, 7)].mean() for s in range(7)]
print(f"\nk=7 계통표집, 시작점별 평균:")
print(f"  {np.round(sys7, 2)}")
print(f"  표준편차 {np.std(sys7):.4f}")

sys8 = [y[np.arange(s, N, 8)].mean() for s in range(8)]
print(f"\nk=8 계통표집, 시작점별 평균:")
print(f"  {np.round(sys8, 2)}")
print(f"  표준편차 {np.std(sys8):.4f}")

srs = np.array([y[rng.choice(N, 1000, replace=False)].mean() for _ in range(5000)])
print(f"\nSRS n=1000: 평균 {srs.mean():.4f}   표준오차 {srs.std(ddof=1):.4f}")

모평균 55.6948

k=7 계통표집, 시작점별 평균:
  [49.98 50.31 49.98 49.82 50.11 89.89 49.77]
  표준편차 13.9603

k=8 계통표집, 시작점별 평균:
  [55.76 55.67 55.78 55.78 55.81 55.39 55.58 55.79]
  표준편차 0.1361

SRS n=1000: 평균 55.7078   표준오차 0.4347

import numpy as np

rng = np.random.default_rng(2)
N = 100_000
W_true = np.array([0.40, 0.35, 0.25])              # 알려진 모집단 연령 분포
age = rng.choice([0, 1, 2], N, p=W_true)
y = np.array([30., 50., 70.])[age] + rng.normal(0, 10, N)

p_resp = np.array([0.10, 0.25, 0.50])[age]         # 나이 많을수록 응답률↑

B = 5000
naive = np.empty(B)
post = np.empty(B)
for b in range(B):
    pool = rng.choice(N, 4000, replace=False)
    r = pool[rng.random(4000) < p_resp[pool]]      # 실제 응답자
    naive[b] = y[r].mean()
    post[b] = sum(W_true[h] * y[r][age[r] == h].mean() for h in range(3))

print(f"참 모평균 {y.mean():.4f}")
print(f"단순평균   {naive.mean():.4f}   편향 {naive.mean() - y.mean():+.4f}"
      f"   표준오차 {naive.std(ddof=1):.4f}")
print(f"사후층화   {post.mean():.4f}   편향 {post.mean() - y.mean():+.4f}"
      f"   표준오차 {post.std(ddof=1):.4f}")

참 모평균 47.0006
단순평균   56.7327   편향 +9.7320   표준오차 0.5630
사후층화   46.9737   편향 -0.0270   표준오차 0.3868

import numpy as np

rng = np.random.default_rng(3)
N, n, B = 5000, 200, 20_000
size = rng.lognormal(3, 0.8, N)          # 사업체 규모
pi = n * size / size.sum()               # 규모비례확률
cum = np.cumsum(pi)
madow = lambda u: np.searchsorted(cum, u + np.arange(n))   # 고정크기 PPS

cases = [("y 가 규모에 거의 비례", 2.5 * size + rng.normal(0, 3, N)),
         ("y 가 규모와 무관",      np.full(N, 50.) + rng.normal(0, 20, N))]

for label, y in cases:
    T = y.sum()
    ht = np.empty(B); naive = np.empty(B); srs = np.empty(B)
    for b in range(B):
        s = madow(rng.random())
        ht[b] = (y[s] / pi[s]).sum()          # 호비츠–톰프슨
        naive[b] = y[s].mean() * N            # 가중을 잊은 경우
        srs[b] = y[rng.choice(N, n, replace=False)].mean() * N
    print(f"{label}   (참 총계 {T:,.0f})")
    for nm, v in [("단순 확대", naive), ("PPS + HT", ht), ("SRS", srs)]:
        print(f"   {nm:>9}: 편향 {v.mean() - T:>+10,.0f}   표준오차 {v.std(ddof=1):>10,.0f}")
    print(f"   → HT 의 표준오차는 SRS 의 {ht.std(ddof=1) / srs.std(ddof=1):.2f}배\n")

y 가 규모에 거의 비례   (참 총계 347,672)
       단순 확대: 편향   +309,513   표준오차     23,963
    PPS + HT: 편향        +19   표준오차      1,438
         SRS: 편향       +108   표준오차     22,740
   → HT 의 표준오차는 SRS 의 0.06배

y 가 규모와 무관   (참 총계 251,386)
       단순 확대: 편향       -786   표준오차      6,175
    PPS + HT: 편향       +175   표준오차     19,796
         SRS: 편향        -11   표준오차      6,796
   → HT 의 표준오차는 SRS 의 2.91배