표준오차¶

\(X_i\)를 \(i\)번째 사람의 물 소비량이라 하자. 독립을 가정하면 중심극한정리에 의해 표본평균 \(\bar{X}\)는 근사적으로 평균 2, 표준편차 \(0.7/\sqrt{50} \approx 0.0990\)인 정규분포를 따른다.

(1) 이론값. \(\text{Uniform}(0,1)\)은 \(\mu = 1/2\), \(\sigma = 1/\sqrt{12}\)이므로

이다. 1만 번 되풀이할 때 이 두 추정값이 흔들리는 폭은 각각 \(\operatorname{SE}/\sqrt{B} = 0.001291\)과 \(\operatorname{SE}/\sqrt{2B} = 0.000913\)이므로, 소수 셋째 자리까지만 맞기를 기대해야 한다.

화살표가 가두는 확률은 \(n = 5\)에서 정확히 \(0.6724\)다(표본평균 쪽의 보기 7에서 어윈-홀 분포로 계산했다). 정규분포의 \(0.6827\)보다 \(0.0103\) 작다.

(2) 실행.

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)

def main():
    # 크기 5짜리 균등표본을 1만 번 뽑아 표본평균을 모은다.
    X_bar = []
    for _ in range(10_000):
        x = np.random.uniform(size=(5,))
        x_bar = x.mean()
        X_bar.append(x_bar)

    # 표준오차의 정의를 그대로 실행한 것이 아래 두 줄이다.
    #   average        = 표집분포의 중심 (참 mu = 0.5 의 좋은 추정)
    #   standard_error = 표집분포의 표준편차 <- 이것이 표준오차다
    # 즉 표준오차는 새로운 개념이 아니라, 통계량의 분포에 대한 표준편차다.
    # 현실에서는 표본이 하나뿐이라 이렇게 구할 수 없어 공식 s/sqrt(n) 을 쓴다.
    X_bar = np.array(X_bar)
    average = X_bar.mean()
    standard_error = X_bar.std()

    print(f'표본평균의 표집분포 중심 : {average:.4}')
    print(f'표본평균의 표준오차      : {standard_error:.4}')

    fig, ax = plt.subplots(figsize=(12, 3))
    # 한글은 $...$ 바깥에 둔다. 수식 글꼴에는 한글 글리프가 없다.
    ax.set_title(r'표본평균 $\bar X$ 의 표집분포', fontsize=20)

    ax.hist(X_bar, bins=100, density=True, alpha=0.3)
    ax.vlines(average, ymin=0, ymax=5, alpha=1.0, color='k', ls='-', lw=5)
    ax.vlines(average + standard_error, ymin=0, ymax=5, alpha=0.7, color='k', ls='--')
    ax.vlines(average - standard_error, ymin=0, ymax=5, alpha=0.7, color='k', ls='--')

    # 평균에서 +1 표준오차까지를 양방향 화살표로 표시해
    # "표준오차 = 이만큼의 폭"임을 그림에서 직접 보여 준다.
    ax.annotate('', xy=(average, 5), xytext=(average + standard_error, 5),
                arrowprops=dict(arrowstyle='<->', color='k', linewidth=3,
                                mutation_scale=20))
    ax.text(average, 5.5, '표준오차', fontsize=15)

    ax.set_xlim(0.0, 1.0)
    ax.set_ylim(-0.1, 6)

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    plt.show()

if __name__ == "__main__":
    main()

출력:

표본평균의 표집분포 중심 : 0.4981
표본평균의 표준오차      : 0.1287

두 숫자가 띠 안에 든다. 중심 \(0.4981\)이 \(0.5\)에서 \(0.0019\) 떨어져 있어 몬테카를로 오차 \(0.001291\)의 \(1.47\)배, 표준오차 \(0.1287\)이 \(0.129099\)에서 \(0.0004\) 떨어져 있어 \(0.000913\)의 \(0.44\)배다. 표준오차는 새로 배울 개념이 아니라 "통계량 값을 잔뜩 모아 잰 표준편차"라는 것이 이 두 줄의 내용이다.

화살표의 폭. 평균에서 한 표준오차까지가 화살표이고, 양쪽을 합한 구간 \([\bar X - \operatorname{SE},\ \bar X + \operatorname{SE}]\)이 가두는 확률은 \(0.6724\)다. 정규분포라면 \(0.6827\)이니 \(1\) 퍼센트포인트 작다. \(n = 5\)의 균등 표본평균은 꼬리가 정규보다 짧아(\(\bar X\)의 초과첨도 \(-0.24\)) 한 표준오차 안에 든 질량이 오히려 적다.

현실에서는 이 그림을 그릴 수 없다는 점을 잊지 말아야 한다. 표본은 한 번만 뽑히므로 \(\bar X\) 값 하나밖에 손에 없고, 1만 개를 모아 표준편차를 재는 일은 불가능하다. \(s/\sqrt n\)이라는 공식이 이 그림을 대신하는 장치다.

유도할 식이 없는 보기다. 여기서 읽을 것은 수가 아니라 구조다.

(1) import 가 곧 실행이다. 파이썬에서 모듈을 처음 가져오면 그 파일의 최상위 문장이 전부 한 번 실행된다. 이 파일의 최상위에는 세 가지가 있다.

1. ARGS = parser.parse_args() — 명령줄을 읽어 --seed 값을 정한다.
2. np.random.seed(ARGS.seed) — 전역 난수 상태를 그 값으로 고정한다.
3. plt.rcParams[...] = ... — 글꼴과 마이너스 기호 설정을 바꾼다.

그러므로 이 파일을 가져오는 행위 자체가 "씨앗 고정 + 그림 설정"을 수행한다. 뒤의 두 파일은 np.random.seed를 한 줄도 쓰지 않는데도 재현 가능한 결과를 내는데, 그 일을 여기서 대신 해 주기 때문이다.

씨앗을 한 곳에서 한 번만 정한다는 점이 핵심이다. 파일마다 np.random.seed(1)을 적어 두면 어느 한 곳을 고치는 순간 조용히 어긋나기 시작한다. 모듈 하나로 모아 두면 그럴 수 없다. 설정이 여러 군데 흩어지지 않게 한다는 뜻에서 단일 진실 공급원을 둔 것이다.

(2) 쓰이지 않는 이름을 가져오는 까닭. from global_name_space import ARGS라고 적으면 린터는 "쓰지 않는 import"라고 경고할 것이다. 그러나 이 줄을 지우면 모듈이 아예 로드되지 않고, 따라서 씨앗도 글꼴도 적용되지 않는다. 가져오는 이름이 아니라 가져오는 행위가 목적인 import다. 이런 것을 부작용을 노린 import 라 부르며, 드물지만 설정 모듈에서는 흔한 관용이다.

같은 이유로 import global_name_space라고만 적어도 효과는 같다. ARGS를 명시하는 것은 "이 파일이 설정을 들고 있다"는 신호를 코드에 남기려는 선택이다.

깨지는 경우. parse_args()를 모듈 최상위에서 부르는 것이 값이다. 이 파서는 --seed만 알고 있으므로, 이 모듈을 가져오는 어떤 프로그램이든 자기 명령줄 인자를 그 파서에 내주게 된다.

$ python standard_error_of_s_square.py --bins 50
usage: standard_error_of_s_square.py [-h] [--seed S]
standard_error_of_s_square.py: error: unrecognized arguments: --bins 50

standard_error_of_s_square.py는 --bins라는 옵션을 알지 못하지만, 그 오류를 내는 주체도 그 파일이 아니다. global_name_space가 가진 파서다. 같은 이유로 pytest나 주피터에서 이 모듈을 가져오면 그쪽의 명령줄 인자를 읽으려 들어 SystemExit로 죽는다.

고치는 방법은 파싱을 함수 안으로 내리고(def get_args(): ...) 호출 시점을 if __name__ == "__main__": 안으로 옮기는 것이다. 다만 그러면 씨앗 고정도 명시적으로 불러야 하므로 "가져오기만 하면 설정이 끝난다"는 편의는 사라진다. 편의와 안전을 맞바꾸는 자리이고, 이 쪽에서는 보기 셋이 같은 씨앗을 쓰는 것을 보이려고 편의 쪽을 택했다.

import argparse

import matplotlib.pyplot as plt
import numpy as np

# 이 파일은 여러 스크립트가 공유하는 설정을 한곳에 모아 두는 용도다.
# 다른 모듈에서 `from global_name_space import ARGS` 로 가져다 쓴다.
parser = argparse.ArgumentParser(description='Standard error simulation')
parser.add_argument('--seed', type=int, default=1, metavar='S',
                    help='random seed (default: 1)')
ARGS = parser.parse_args()

# 시드를 여기서 한 번만 고정하면 이 설정을 가져다 쓰는 모든 스크립트가
# 같은 난수열을 쓰게 되어 결과가 재현된다.
np.random.seed(ARGS.seed)

# 그림 설정도 여기 모아 둔다. 이 파일을 가져다 쓰는 스크립트는 글꼴을
# 따로 손대지 않아도 된다.
# 그림에 한글을 쓰므로 한글 글꼴을 지정한다. 맥이면 'Apple SD Gothic Neo',
# 윈도우면 'Malgun Gothic', 리눅스면 'NanumGothic' 정도가 무난하다.
# 글꼴을 바꾸면 마이너스 기호가 깨지므로 unicode_minus 도 함께 꺼 준다.
plt.rcParams['font.family'] = 'Apple SD Gothic Neo'
plt.rcParams['axes.unicode_minus'] = False

(1) 판정 기준. 보기 1은 np.random.seed(0), 이 파일은 global_name_space.py의 기본값 --seed 1을 쓴다. 뽑히는 난수열이 다르니 결과가 달라지는 것이 정상이다.

문제는 "얼마나 달라도 정상인가"다. 두 실행은 같은 양 \(\mu = 0.5\)와 \(\operatorname{SE} = 1/\sqrt{60} = 0.129099\)를 각각 1만 개의 표본평균으로 추정한다. 추정값의 요동은

이다. 그러므로 두 실행의 중심이 \(0.0013\)의 두세 배 안에서 다르고 표준오차가 \(0.0009\)의 두세 배 안에서 다르면 요동이고, 그보다 크게 벌어지면 코드가 다른 것이다. 각 실행이 참값에서 얼마나 떨어졌는지를 재면 바로 판정된다.

(2) 실행.

import matplotlib.pyplot as plt
import numpy as np

# 이 import 는 ARGS 라는 이름이 필요해서가 아니다. global_name_space 를 불러오는
# 순간 그 파일의 np.random.seed(ARGS.seed) 와 글꼴 설정이 실행된다.
# 즉 "가져오는 것" 자체가 설정을 적용하는 일이다.
from global_name_space import ARGS

def draw(values, name, symbol, ymax, ytop, xlim):
    """표집분포의 히스토그램에 평균선과 ±1 표준오차를 그린다.

    표본평균이든 표본분산이든 재는 방법도 그리는 방법도 같으므로
    한곳에 모아 둔다. values 에 무엇을 담아 넘기느냐만 바뀐다.
    """
    values = np.array(values)
    average = values.mean()
    standard_error = values.std()

    print(f'{name}의 표집분포 중심 : {average:.4}')
    print(f'{name}의 표준오차      : {standard_error:.4}')

    fig, ax = plt.subplots(figsize=(12, 3))
    # 한글은 $...$ 바깥에 둔다. 수식 글꼴에는 한글 글리프가 없다.
    ax.set_title(f'{name} ${symbol}$ 의 표집분포', fontsize=20)

    ax.hist(values, bins=100, density=True, alpha=0.3)
    ax.vlines(average, ymin=0, ymax=ymax, alpha=1.0, color='k', ls='-', lw=5)
    ax.vlines(average + standard_error, ymin=0, ymax=ymax, alpha=0.7, color='k', ls='--')
    ax.vlines(average - standard_error, ymin=0, ymax=ymax, alpha=0.7, color='k', ls='--')

    # 평균에서 +1 표준오차까지를 양방향 화살표로 표시해
    # "표준오차 = 이만큼의 폭"임을 그림에서 직접 보여 준다.
    ax.annotate('', xy=(average, ymax), xytext=(average + standard_error, ymax),
                arrowprops=dict(arrowstyle='<->', color='k', linewidth=3,
                                mutation_scale=20))
    ax.text(average, ymax * 1.1, '표준오차', fontsize=15)

    ax.set_xlim(*xlim)
    ax.set_ylim(-0.1, ytop)

    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    plt.show()

def main():
    X_bar = []
    for _ in range(10_000):
        x = np.random.uniform(size=(5,))
        x_bar = x.mean()
        X_bar.append(x_bar)

    draw(X_bar, '표본평균', r'\bar X', ymax=5, ytop=6, xlim=(0.0, 1.0))

if __name__ == "__main__":
    main()

출력:

표본평균의 표집분포 중심 : 0.4993
표본평균의 표준오차      : 0.1281

요동이다. 이 실행의 중심 \(0.4993\)은 \(0.5\)에서 \(0.00067\) 떨어져 있어 오차 \(0.001291\)의 \(0.52\)배, 표준오차 \(0.1281\)은 \(0.129099\)에서 \(0.00100\) 떨어져 있어 \(0.000913\)의 \(1.09\)배다. 보기 1은 각각 \(1.47\)배와 \(0.44\)배였다. 네 값이 모두 \(1.5\)배 안에 들어 있으므로 두 실행은 같은 것을 재고 있다.

씨앗을 바꿔 가며 몇 번 더 돌려 보면 감이 잡힌다. --seed 7이면 중심 \(0.4996\), 표준오차 \(0.1283\)이 나온다. 세 실행의 표준오차가 \(0.1287\), \(0.1281\), \(0.1283\)으로 셋째 자리에서 흔들리고 넷째 자리는 의미가 없다. 1만 번으로 살 수 있는 정밀도가 거기까지라는 것이 (1)에서 계산한 \(0.0009\)의 뜻이다.

draw()를 뗀 이득. 이 함수가 받는 것은 "값의 목록, 이름, 기호, 축 범위"뿐이고 값이 표본평균인지 표본분산인지 묻지 않는다. 표준오차를 재는 두 줄

average = values.mean()
standard_error = values.std()

도 통계량의 종류와 무관하다. 표준오차의 정의가 대상을 가리지 않기 때문에 코드도 가릴 필요가 없다. 그래서 다음 보기의 파일이 열 줄 남짓으로 끝난다. 모듈로 나눈 것의 값은 파일 수를 늘린 데 있지 않고, "재는 방법이 하나"라는 사실을 코드 구조로 드러낸 데 있다.

(1) 이론값. \(\text{Uniform}(0,1)\)은 \(\sigma^2 = 1/12\), \(\mu_4 = 1/80\), 첨도 \(\beta_2 = \mu_4/\sigma^4 = 1.8\)이다. 불편성에서

이고, 폭은 \(n = 5\)에서 \((n-3)/(n-1) = 1/2\)이므로

다. 1만 번 되풀이할 때의 요동은 중심이 \(0.042492/100 = 0.000425\), 표준오차가 \(0.042492/\sqrt{20000} = 0.000300\)이다.

눈여겨볼 비가 하나 있다. \(\operatorname{SE}(S^2)/E[S^2] = 0.042492/0.083333 = 0.51\)이다. 분산을 재는 일이 평균을 재는 일보다 훨씬 거칠다. 같은 \(n = 5\)에서 \(\bar X\)의 상대적 흔들림은 \(0.129099/0.5 = 0.26\)으로 그 절반이다.

ddof를 빼먹으면. x.var()는 \(n\)으로 나누므로 \(\frac{n-1}{n}S^2\)을 돌려준다. 따라서 중심이

로 나온다. 참값보다 \(20\%\) 작다. \(0.0833\)이 아니라 \(0.0667\)이 찍히면 불편성이 깨진 것처럼 보이지만 깨진 것은 ddof다.

(2) 실행.

import numpy as np

# 그리는 부분은 앞 파일에 이미 있으므로 가져다 쓴다.
# 앞 파일에 if __name__ == "__main__": 가드가 있어서, import 해도
# 그쪽 main() 은 실행되지 않고 함수 정의만 넘어온다. 가드를 두는 이유가 이것이다.
# 이 import 가 global_name_space 도 함께 불러오므로 시드도 같이 고정된다.
from standard_error_of_x_bar import draw


def main():
    S_square = []
    for _ in range(10_000):
        x = np.random.uniform(size=(5,))
        # ddof=1 을 반드시 적어야 한다. numpy의 기본값은 ddof=0 이라
        # x.var() 나 x.std()**2 는 n 으로 나눈 편향추정량을 준다.
        S_square.append(x.var(ddof=1))

    draw(S_square, '표본분산', 'S^2', ymax=8, ytop=10, xlim=(0.0, 0.25))

if __name__ == "__main__":
    main()

출력:

표본분산의 표집분포 중심 : 0.08406
표본분산의 표준오차      : 0.04263

두 숫자가 맞는다. 중심 \(0.08406\)이 이론 \(0.083333\)에서 \(0.00073\) 떨어져 있어 오차 \(0.000425\)의 \(1.72\)배, 표준오차 \(0.04263\)이 이론 \(0.042492\)에서 \(0.00014\) 떨어져 있어 \(0.000300\)의 \(0.45\)배다. 중심 쪽이 \(1.7\)배까지 벌어진 것은 \(S^2\)의 표집분포가 오른쪽으로 치우쳐 있어 1만 번으로도 평균이 안정되기 더 어렵기 때문이고, \(2\)배 안쪽이므로 어긋났다고 볼 근거는 없다.

씨앗을 바꾸면 흔들림이 보인다. --seed 7로 돌리면 중심 \(0.0842\), 표준오차 \(0.04259\)가 나온다. 표준오차 쪽은 넷째 자리까지 꽤 안정적인데 중심 쪽이 더 흔들린다는 것이 위 두 오차의 비(\(0.000425\) 대 \(0.000300\))와 치우침이 함께 만든 결과다.

모양이 다른 까닭. draw()도, 재는 두 줄도, 씨앗도, 되풀이 횟수도 같다. 바뀐 것은 values에 무엇을 담았는지뿐인데 그림이 전혀 다르다.

• 보기 3의 \(\bar X\)는 \(0.5\)를 중심으로 대칭이다. 균등분포의 왜도가 \(0\)이고 \(\bar X\)의 왜도는 \(\gamma_1/\sqrt n\)이므로 \(0\)이다.
• 보기 4의 \(S^2\)은 오른쪽으로 끌려 있다. 모집단이 대칭이어도 \(S^2 \ge 0\)이라 왼쪽이 \(0\)에서 막히고 오른쪽만 열려 있기 때문이다. 모의 왜도는 \(+0.41\)이다.

그래서 점선 두 개를 "전형적인 범위"로 읽을 때 두 그림의 성격이 다르다. 보기 3에서는 그 구간이 대칭인 분포를 대칭으로 가둔 것이라 자연스럽지만, 보기 4에서는 \(\mu - 2\operatorname{SE} = 0.0833 - 0.0850 < 0\)이어서 두 표준오차만 내려가도 가능한 값의 바깥이다. 같은 도구가 대상에 따라 다르게 읽혀야 한다는 뜻이며, 표본분산의 표준오차 쪽에서 이 치우침을 수로 갈라 센다.

(a) \(\mathrm{SE}_{25} = 50/\sqrt{25} = 10\). \(\mathrm{SE}_{100} = 50/\sqrt{100} = 5\). \(n\)이 네 배가 되면 표준오차가 절반이 된다.

(b) \(\sigma/\sqrt{16} = 5\)에서 \(\sigma = 20\)이다. \(n = 64\)에서 \(\mathrm{SE} = 20/\sqrt{64} = 2.5\)이다. \(n\)을 네 배로 하면 표준오차가 절반이 된다.

(a) (\(\sigma\) 대신 \(s\)를 사용한) 추정 표준오차: \(\mathrm{SE} = 8/\sqrt{100} = 0.8\).

(b) SD = 8: 자료 집합에서 개별 관측값이 얼마나 달라지는지를 기술한다. 전형적인 개체는 평균에서 약 8단위 떨어져 있다.

SE = 0.8: 표본이 달라질 때 표본평균이 얼마나 달라지는지를 기술한다. 참 모평균은 50에서 약 1.6단위(≈ 표준오차 2개) 이내에 있을 가능성이 높다.

자료를 더 모아도 표준편차는 달라지지 않지만 표준오차는 \(1/\sqrt n\)의 비율로 줄어든다. 표준오차를 말해야 할 자리에 표준편차를 보고하거나 그 반대로 하는 것은 과학 논문에서 흔한 오류이다.

i.i.d. 표본 \(X_1, \ldots, X_n\)이 주어졌을 때:

1. 원래 표본에서 복원추출하여 붓스트랩 표본 \(X_1^*, \ldots, X_n^*\)을 뽑는다.
2. 붓스트랩 표본에서 \(\bar X^*\)를 계산한다.
3. 1–2단계를 \(B\)번 반복하여(보통 \(B = 1000\)에서 10000) \(\bar X^*_1, \ldots, \bar X^*_B\)을 얻는다.
4. 붓스트랩 복제값들의 표본표준편차로 표준오차를 추정한다: \(\hat{\mathrm{SE}}_{\text{boot}} = \sqrt{(1/(B-1))\sum(\bar X^*_b - \bar X^*_\cdot)^2}\).

왜 작동하는가: 붓스트랩 분포가 반복추출에서의 표본분포를 근사한다. 점근적으로 \(\hat{\mathrm{SE}}_{\text{boot}} \to \sigma/\sqrt n\)이지만, 붓스트랩은 정규근사보다 분포의 모양(치우침, 두꺼운 꼬리)을 더 잘 포착한다.

닫힌 형태의 표준오차가 없는 경우(중앙값, 비, 복잡한 모형의 회귀계수)에 특히 유용하다.

델타 방법: \(\sqrt n(\hat\theta - \theta) \xrightarrow{d} N(0, \sigma^2)\)이면, \(g'(\theta) \ne 0\)인 미분가능한 \(g\)에 대해:

표준오차 형태로는 \(\mathrm{SE}(g(\hat\theta)) \approx |g'(\hat\theta)| \cdot \mathrm{SE}(\hat\theta)\)이다.

예: \(g(\hat p) = \log(\hat p/(1 - \hat p))\)(로짓)에 대해 \(g'(\hat p) = 1/(\hat p(1 - \hat p))\)이므로 \(\mathrm{SE}(\hat\eta) = \mathrm{SE}(\hat p)/(\hat p(1 - \hat p))\)이다.

추정량이 표준오차를 아는 더 단순한 추정량의 변환일 때, 델타 방법은 표준오차 계산의 주된 도구가 된다.

독립성에 의해 \(\mathrm{Var}(\bar X_1 - \bar X_2) = \sigma_1^2/n_1 + \sigma_2^2/n_2\)이다.

(a) 합동 (\(\sigma_1 = \sigma_2 = \sigma\) 가정):

합동 분산: \(s_p^2 = ((n_1 - 1)s_1^2 + (n_2 - 1)s_2^2)/(n_1 + n_2 - 2)\).

\(\mathrm{SE}_{\text{pool}} = s_p \sqrt{1/n_1 + 1/n_2}\).

분산이 같다고 볼 만할 때 표준적인 두 표본 \(t\) 검정에서 사용한다. 이 가정이 성립하면 약간 더 효율적이다.

(b) Welch (이분산):

\(\mathrm{SE}_{\text{Welch}} = \sqrt{s_1^2/n_1 + s_2^2/n_2}\).

합동하지 않고 각 표본이 자신의 분산을 기여한다. \(t\) 임계값에는 (정수가 아닌) Welch 자유도를 사용한다.

현대적 권고: 기본적으로 Welch를 선호하라. 강한 등분산 가정을 요구하지 않으면서 분산이 같을 때에도 거의 같은 성능을 낸다.

유한모집단에서 비복원추출을 하면:

따라서 \(\mathrm{SE}(\bar X) = (\sigma/\sqrt n) \sqrt{1 - n/N}\)이다. 인수 \(\sqrt{1 - n/N}\)이 유한모집단 수정(FPC)이다.

극한:

• \(n/N \to 0\) (추출 비율이 매우 작을 때): FPC \(\to 1\)로 표준적인 \(\sigma/\sqrt n\)이 복원된다. 전국 조사(\(n = 1000, N \approx 10^8\))에 해당한다.
• \(n/N \to 1\) (전수조사): FPC \(\to 0\)으로 표본추출 변동성이 없다. 전수조사는 결정론적 추정값을 준다.

FPC가 중요한 경우: 감사(송장 500건 중 100건 추출, \(n/N = 0.2\), FPC \(\approx 0.89\)). FPC가 신뢰구간을 약 11% 좁히므로 무시할 수 없다.

입문 통계학의 공식 대부분이 FPC를 무시하는 이유는 일반적인 과학 표본에서 \(n/N\)이 작기 때문이다.

\(\hat\beta_1 = \sum_i w_iy_i\)(\(w_i = (x_i-\bar x)/S_{xx}\), \(S_{xx}=\sum(x_i-\bar x)^2\))로 쓸 수 있으므로 \(y_i\)가 독립이고 분산이 \(\sigma^2\)이면

이다. \(S_{xx} = (n-1)s_x^2\)이므로 위 식이 나온다. \(\square\)

설계에 주는 조언. 표준오차를 줄이는 길이 세 가지다.

1. \(n\)을 늘린다. \(\sqrt{n-1}\)에 반비례하므로 익숙한 \(1/\sqrt n\) 규칙이다.
2. 오차분산 \(\sigma\)를 줄인다. 측정 정밀도를 높이거나, 다른 설명변수를 넣어 잔차를 줄이거나, 실험 조건을 통제한다.
3. \(x\)를 넓게 퍼뜨린다. \(s_x\)에 반비례한다. 이것이 관측연구에는 없고 실험 설계에만 있는 지렛대다.

세 번째가 특히 중요하다. 관측이 \(n\)개로 고정되어 있다면, \(x\) 값을 가능한 범위의 양 끝에 몰아 놓는 것이 \(S_{xx}\)를 최대로 하고 따라서 기울기를 가장 정밀하게 추정한다. 약물 용량 반응 실험에서 중간 용량을 생략하고 최저·최고 용량에 표본을 몰아 배치하는 설계가 그 예다.

다만 이 설계에는 대가가 있다. 양 끝만 보면 관계가 직선인지 확인할 수 없다. 굽은 관계를 직선으로 잘못 읽어도 알아챌 방법이 없다. 그래서 실무에서는 중간에도 일부 배치해 모형 적합성을 점검할 여지를 남긴다. 최적 설계와 모형 진단 사이의 맞바꿈이다.

\(n=25\)에서 세 막대의 반길이는 다음과 같다.

SD 막대가 SE 막대의 5배이고, 신뢰구간은 SE의 2.06배다. 같은 자료인데 그림의 인상이 완전히 달라진다.

왜 반드시 밝혀야 하는가. 세 막대가 서로 다른 질문에 답하기 때문이다.

• SD 막대: "개별 관측값이 얼마나 흩어져 있는가." \(n\)이 커져도 줄지 않는다. 자료의 변동성 자체를 보여 주고 싶을 때 쓴다.
• SE 막대: "평균을 얼마나 정확히 알고 있는가." \(n\)과 함께 줄어든다.
• 신뢰구간: "참 평균이 어디에 있을 법한가." 해석이 가장 직접적이다.

독자가 SE 막대를 SD로 오해하면 자료가 실제보다 훨씬 균일하다고 믿게 되고, 반대로 오해하면 추정이 실제보다 부정확하다고 믿게 된다.

더 흔한 함정. "두 집단의 SE 막대가 겹치니 차이가 유의하지 않다"는 판단은 틀렸다. 두 집단의 표본크기와 분산이 같을 때 차이의 표준오차는 \(\sqrt2\operatorname{SE}\)이므로, 두 SE 막대가 살짝 겹쳐도 \(t\)가 2를 넘을 수 있다. 대략 두 SE 막대의 끝이 서로 상대의 평균에 닿지 않을 정도로 겹치면 아직 유의할 수 있다.

반대로 95% 신뢰구간 막대가 겹치지 않으면 차이는 거의 확실히 유의하다(보수적인 판단이다). 겹치더라도 유의할 수 있으므로, 결국 그림으로 유의성을 판정하지 말고 차이에 대한 구간을 따로 보고하는 것이 옳다.

군집 크기가 \(m=20\)이므로 설계효과는

이고 유효 표본크기는

이다. 1000명을 조사했지만 독립 표본 513명만큼의 정보밖에 없다.

과소평가의 크기. 분산이 1.95배이므로 표준오차는 \(\sqrt{1.95} = 1.40\)배다. 독립을 가정하면 표준오차를 29% 과소평가한다(\(1 - 1/1.40\)). 그 결과

• 신뢰구간이 실제보다 29% 좁아지고,
• 검정통계량이 1.40배 부풀려지며,
• 명목 5% 검정의 실제 오류율이 15% 안팎으로 뛴다.

ICC가 작아 보여도 방심할 수 없다. 0.05는 낮은 값인데도 설계효과가 2에 가깝다. 괄호 안이 \((m-1)\times\text{ICC}\)이므로 군집이 클수록 영향이 커지기 때문이다. 같은 ICC라도 학급당 5명이면 DEFF가 1.2에 그치고, 100명이면 5.95가 된다.

설계에 주는 함의. 같은 예산이라면 군집을 크게 하기보다 군집 수를 늘리는 편이 낫다. 학교 50곳에서 20명씩 뽑는 것보다 학교 100곳에서 10명씩 뽑는 쪽이 유효 표본크기가 크다(\(\text{DEFF} = 1.45\), \(n_{\text{eff}} = 690\)). 다만 학교를 추가하는 비용이 학생을 추가하는 비용보다 훨씬 크므로, 실제로는 두 비용의 비와 ICC를 함께 넣어 최적 군집 크기를 계산한다.

일반식.

이고, 지연 \(k\)인 쌍이 \(2(n-k)\)개(\(k\ge1\))이므로

이다.

\(n\)이 클 때. \(\sum_k|\rho_k| < \infty\)이면 \(1-k/n \to 1\)이므로

이다. 괄호 안이 설계효과에 해당하며 적분시간척도라 부른다.

AR(1). \(\rho_k = \phi^k\)이므로 등비급수에서

이다. 따라서

\(\phi=0.8\)이면

로 분산이 아홉 배다. 표준오차는 3배이고, 유효 표본크기가 \(n/9\)로 줄어든다. 관측 900개가 독립 관측 100개만큼의 정보밖에 없다.

3장에서 본 \(\phi = 0.9\)가 가운데 있다. 관측 만 개가 독립 관측 526개 값어치로 줄어드는 그 경우다.

실무적 함의.

• 시계열을 독립 표본처럼 다루면 표준오차를 심하게 과소평가한다. \(\phi=0.95\)라면 세 배가 아니라 여섯 배 과소평가한다.
• \(\phi < 0\)(음의 자기상관)이면 반대로 \(\operatorname{Var}(\bar X)\)가 줄어든다. 계통표집이 단순무작위표집보다 나을 수 있는 이유가 이것이다.
• MCMC 표본이 정확히 이 구조를 갖는다. 그래서 유효표본크기(ESS)를 보고하고, 사슬을 솎아 내거나(thinning) 더 길게 돌리는 것이다.
• 실무에서는 \(\rho_k\)를 모르므로 뉴이-웨스트 추정량처럼 표본 자기공분산을 가중합한 HAC 표준오차를 쓴다.

표본 \(p\) 분위수 \(\hat{x}_p\)는 점근적으로

를 따른다. 즉 표준오차가

이다. 분자는 \(p = 0.5\)에서 최대이지만, 분모의 밀도 \(f(x_p)\)가 꼬리에서 급격히 작아지므로 극단 분위수의 표준오차가 훨씬 크다.

import numpy as np
from scipy import stats

rng = np.random.default_rng(1)
B, n = 40_000, 500
d = rng.normal(0, 1, (B, n))
se_mean = 1 / np.sqrt(n)

print(f"N(0,1), n={n}")
print(f"{'p':>6}{'모의 SE':>11}{'이론 SE':>11}{'평균의 SE 대비':>15}")
for p in (0.5, 0.75, 0.90, 0.95, 0.99):
    q = np.quantile(d, p, axis=1)
    theory = np.sqrt(p * (1 - p) / n) / stats.norm.pdf(stats.norm.ppf(p))
    print(f"{p:>6.2f}{q.std():>11.5f}{theory:>11.5f}{q.std() / se_mean:>15.3f}")

출력:

N(0,1), n=500
     p      모의 SE      이론 SE      평균의 SE 대비
  0.50    0.05590    0.05605          1.250
  0.75    0.06114    0.06094          1.367
  0.90    0.07609    0.07645          1.701
  0.95    0.09410    0.09450          2.104
  0.99    0.15955    0.16696          3.568

\(p \le 0.95\)에서는 이론과 모의가 소수점 셋째 자리까지 맞는다. \(p = 0.99\)만 \(0.160\) 대 \(0.167\)로 벌어지는데, \(n = 500\)에서 상위 \(1\%\)에 해당하는 관측이 다섯 개뿐이라 점근 공식이 아직 자리를 잡지 못한 탓이다. 극단 분위수가 어렵다는 이 문제의 요지가 이론식 자체에도 그대로 나타난 셈이다.

\(99\)번째 백분위수를 중앙값만큼 정밀하게 추정하려면 표본이 약 \(8\)배 필요하다(\(3.57^2 \approx 12.7\), 중앙값 대비로는 \((0.1596/0.0559)^2 \approx 8.2\)).

실무적 함의.

• 위험관리의 VaR은 보통 \(99\%\)나 \(99.9\%\) 분위수다. 이 계산이 말하는 바는 그 추정값이 본질적으로 불안정하다는 것이다. \(99.9\%\) 분위수를 안정적으로 추정하려면 관측이 수만 개 필요하다.
• 관측 범위를 넘는 분위수는 추정할 수 없다. \(n = 500\)이면 경험적으로 \(99.8\)번째 백분위수가 최댓값이다. 그보다 극단적인 분위수를 말하려면 반드시 모형 가정(극단값 이론 등) 을 끌어들여야 하며, 그 순간 결과는 가정에 의존하게 된다.
• 중앙값이 특별히 정밀한 것도 아니다. 표본평균보다 \(1.25\)배 나쁘다(정규분포에서 효율 \(2/\pi\)의 제곱근인 \(1/\sqrt{0.637} = 1.25\)).

밀도가 작은 곳은 추정이 어렵다는 것이 일반 원리다. 자료가 드문 영역에 대해 정밀한 진술을 하려면 자료를 훨씬 많이 모으거나, 구조에 대한 가정을 빌려 와야 한다. \(\square\)