기댓값과 선형성¶

\(3.5\)는 주사위가 결코 낼 수 없는 값이다. 기댓값은 "기대되는 값"이 아니라 평균이라는 점을 여기서 확인하고 넘어가는 것이 좋다.

import numpy as np

# 공정한 주사위의 기댓값
values = np.arange(1, 7)
probs = np.ones(6) / 6
expected = np.sum(values * probs)
print(f"E[fair die] = {expected:.4f}")

# 모의실험
np.random.seed(42)
rolls = np.random.randint(1, 7, size=100_000)
print(f"Simulated mean = {rolls.mean():.4f}")

출력:

E[fair die] = 3.5000
Simulated mean = 3.5031

이다. 도착률이 \(\lambda\)면 평균 대기시간이 \(1/\lambda\)라는 익숙한 관계다.

위 표의 세 성질이 이 하나의 상황에서 모두 읽힌다.

• \(E[X + Y] = E[X] + E[Y]\) — 두 회사 주식을 한 주씩 가지고 있으면 내년에 받을 배당금은 \(X + Y\)이고, 그 기댓값은 각각의 기댓값을 더한 것이다. 두 배당금은 강하게 종속이라는 점이 중요하다. 같은 반도체 업황을 타므로 한쪽이 좋은 해는 다른 쪽도 좋기 쉽다. 그래도 식은 그대로 성립한다. 결합분포를 몰라도, 두 회사의 관계를 몰라도 된다.
• \(E[aX] = a\,E[X]\) — 삼성전자를 100주 가지고 있으면 받을 배당금은 \(100X\)이고 기댓값은 \(100\,E[X]\)다. 단위를 원에서 달러로 바꾸는 환산도 상수를 곱하는 일이므로 마찬가지다.
• \(E[a] = a\) — 현금 100만 원처럼 확정된 금액은 확률변수가 아니므로 기댓값이 자기 자신이다. 따라서 주식 한 주와 현금 \(c\)를 함께 들고 있으면 \(E[X + c] = E[X] + c\)다.

이 셋을 합치면 임의의 포트폴리오로 확장된다. 삼성전자 \(a\)주, SK하이닉스 \(b\)주, 현금 \(c\)원을 들고 있을 때 내년에 받을 금액의 기댓값은

이다. 종목 사이의 상관관계를 전혀 따지지 않고 각 종목의 기대 배당금만 알면 끝난다.

다만 기댓값이 같다고 같은 포트폴리오는 아니다. 한 종목에 몰아넣은 것과 두 종목에 나눠 담은 것은 기대 배당금이 같아도 해마다 받는 금액의 출렁임이 다르다. 그 차이를 재려면 분산이 필요하고, 그때는 두 종목이 함께 움직이는지가 결정적으로 중요해진다. 선형성이 무시해도 좋았던 바로 그 종속성이 다음 쪽에서 공분산이라는 이름으로 돌아온다.

(1) 해석적으로. 조화수를 그대로 더하면

이다. 그런데 어림식 \(n\ln n\) 에 넣으면

로 \(29.4\) 나 낮다. 상대오차가 \(13\%\) 이니 "약 \(225\) 번" 과는 거리가 멀다.

빠진 것은 상수항이다. 조화수의 점근전개는

이고, 여기에 \(n\) 을 곱하면

이다. \(n\ln n\) 다음 항 \(n\gamma\) 가 \(n\) 에 비례해서 커지므로 결코 무시할 수 없다. \(n = 50\) 에서 \(50\gamma = 28.86\) 이고, 이것이 방금 본 차이 \(29.4\) 의 거의 전부다. 여기에 상수 \(\tfrac12\) 까지 보태면

로 정확값 \(224.9603\) 과 소수 둘째 자리까지 맞는다. \(n\ln n\) 은 비율로만 맞는 어림이고(\(nH_n / (n\ln n) \to 1\)), 차이는 오히려 커진다.

(2) 몬테카를로 오차로 설명된다. 되풀이 \(10{,}000\) 번으로 추정한 평균의 표준오차는 \(\operatorname{sd}(T)/\sqrt{10000}\) 이다. \(\operatorname{sd}(T)\) 를 모의실험에서 직접 재면 \(62.22\) 이므로

이고, 벗어남 \(225.52 - 224.96 = 0.56\) 은 \(0.90\) 표준오차에 지나지 않는다. 이 문제에서 \(T\) 는 매우 넓게 퍼져 있다는 점이 요점이다. 표준편차가 \(62\) 로 평균의 \(28\%\) 나 되므로, 평균을 소수 첫째 자리까지 맞히려면 되풀이가 \(10{,}000\) 번으로는 모자란다.

(3) 수치적으로.

import numpy as np

def coupon_collector_simulation(n_coupons, n_trials=10_000):
    """쿠폰 수집가 문제: n종을 모두 모으려면 몇 개를 사야 하는가."""
    np.random.seed(42)
    totals = []
    for _ in range(n_trials):
        collected = set()      # set이라 중복은 저절로 걸러진다
        count = 0
        # 종류를 다 모을 때까지 무작위로 하나씩 뽑는다
        while len(collected) < n_coupons:
            collected.add(np.random.randint(0, n_coupons))
            count += 1
        totals.append(count)

    simulated = np.mean(totals)

    # 이론값 유도: 이미 k종을 모았을 때 새 종이 나올 확률은 (n-k)/n 이므로
    # 새 종 하나를 더 얻기까지 기대 횟수는 n/(n-k) 다.
    # 이를 k = 0..n-1 로 모두 더하면
    #   n/n + n/(n-1) + ... + n/1 = n * (1 + 1/2 + ... + 1/n) = n * H_n
    # 기댓값의 선형성 덕분에 각 단계가 독립이 아니어도 그냥 더할 수 있다.
    H_n = sum(1/k for k in range(1, n_coupons + 1))      # 조화수 H_n
    theoretical = n_coupons * H_n

    print(f"n = {n_coupons}")
    print(f"Simulated E[T] = {simulated:.1f}")
    print(f"Theoretical E[T] = n·Hₙ = {theoretical:.1f}")

    # (1) 어림식 n ln n 과 견준다. 빠진 것이 오일러-마스케로니 상수 항이다.
    gamma = 0.5772156649015329
    approx1 = n_coupons * np.log(n_coupons)
    approx2 = n_coupons * (np.log(n_coupons) + gamma)
    approx3 = approx2 + 0.5
    print(f"  정확값      n*H_n              = {theoretical:.4f}")
    print(f"  어림식 1    n*ln n             = {approx1:.4f}"
          f"   (오차 {approx1 - theoretical:+.4f})")
    print(f"  어림식 2    n*(ln n + gamma)   = {approx2:.4f}"
          f"   (오차 {approx2 - theoretical:+.4f})")
    print(f"  어림식 3    + 1/2              = {approx3:.4f}"
          f"   (오차 {approx3 - theoretical:+.4f})")

    # (2) 모의값이 정확값에서 벗어난 것이 몬테카를로 오차로 설명되는가.
    sd = np.std(totals, ddof=1)
    se = sd / np.sqrt(n_trials)
    print(f"  모의 표준편차 = {sd:.2f},  평균의 SE = {se:.4f}")
    print(f"  (모의 - 정확)/SE = {(simulated - theoretical) / se:+.3f}")

coupon_collector_simulation(50)

출력:

n = 50
Simulated E[T] = 225.5
Theoretical E[T] = n·Hₙ = 225.0
  정확값      n*H_n              = 224.9603
  어림식 1    n*ln n             = 195.6012   (오차 -29.3591)
  어림식 2    n*(ln n + gamma)   = 224.4619   (오차 -0.4983)
  어림식 3    + 1/2              = 224.9619   (오차 +0.0017)
  모의 표준편차 = 62.22,  평균의 SE = 0.6222
  (모의 - 정확)/SE = +0.900

\(n\ln n\) 의 오차가 \(-29.36\), 여기에 \(n\gamma\) 를 더하면 \(-0.50\), 상수 \(\tfrac12\) 까지 더하면 \(+0.0017\) 로 줄어든다. 유도한 전개가 항마다 그대로 확인된다. 모의값 \(225.52\) 는 정확값에서 \(0.900\) 표준오차 떨어져 있어 어긋남이 아니다.

어림식을 쓸 때 무엇을 버렸는지 알아야 한다. \(n\ln n\) 은 \(T\) 가 커지는 차수를 말해 주지만 \(n = 50\) 같은 실제 값에서는 \(13\%\) 를 놓친다. 위 본문이 "\(n = 50\) 이면 약 \(225\) 번" 이라 한 것은 어림식이 아니라 정확식 \(nH_n\) 을 쓴 값이다.

(1) 해석적으로. 균등분포의 적률은 한 줄로 나온다. \(k \ge 0\) 에 대해

이므로 \(E[X] = \tfrac12\), \(E[Y] = E[X^2] = \tfrac13\) 이고

이다. 한쪽을 알면 다른 쪽이 완전히 정해지는 극단적인 종속인데도 기댓값은 그냥 더해진다. 선형성의 증명 어디에도 독립이 쓰이지 않았기 때문이다. 합의 기댓값은 결합분포 위에서 \(x + y\) 를 적분한 것이고, 그 적분이 \(x\) 의 적분과 \(y\) 의 적분으로 쪼개지는 것은 덧셈의 성질이지 확률의 성질이 아니다.

(2) 다른 둘은 깨진다. 곱 규칙부터 본다. \(XY = X \cdot X^2 = X^3\) 이므로

로 \(1.5\) 배 차이가 난다. 그 차이가 곧 공분산이다.

분산도 보자. \(E[X^4] = \tfrac15\) 이므로

이고, 교차항을 포함한 참값은

인데 교차항을 빠뜨리면 \(\tfrac1{12} + \tfrac4{45} = 0.1722\) 다. 거의 두 배(1.97배)를 놓친다. 상관계수를 내 보면 얼마나 세게 얽혀 있는지 드러난다.

\(Y\) 가 \(X\) 의 함수인데도 \(\rho\) 가 \(1\) 이 아닌 것은 그 함수가 직선이 아니기 때문이다. 상관계수는 직선 관계만 재므로, 완전한 종속조차 \(0.9682\) 로밖에 잡아내지 못한다. 이 어긋남을 정면으로 다루는 곳이 독립성과 무상관성의 차이다.

(3) 수치적으로.

import numpy as np
import matplotlib.pyplot as plt

def linearity_demonstration():
    """종속인 두 변수에서도 기댓값의 선형성이 성립함을 확인한다.

    E[X + Y] = E[X] + E[Y] 는 X와 Y가 **독립이 아니어도** 성립한다.
    독립이 필요한 것은 곱의 기댓값 E[XY] = E[X]E[Y] 이나
    분산의 덧셈 Var(X+Y) = Var(X) + Var(Y) 쪽이다.
    이 구분이 확률론에서 가장 자주 헷갈리는 지점 중 하나다.
    """
    np.random.seed(42)
    n_sim = 100_000

    # X ~ Uniform(0,1), Y = X^2. X를 알면 Y가 완전히 결정되므로 극단적으로 종속이다.
    X = np.random.rand(n_sim)
    Y = X ** 2

    print("X and Y = X² are dependent, but linearity still holds:")
    print(f"E[X] = {X.mean():.4f} (theoretical: 0.5)")
    print(f"E[Y] = {Y.mean():.4f} (theoretical: 0.3333)")
    print(f"E[X + Y] = {(X + Y).mean():.4f}")
    print(f"E[X] + E[Y] = {X.mean() + Y.mean():.4f}")

    # 곱 규칙은 깨진다. E[XY] = E[X^3] = 1/4 이지만 E[X]E[Y] = 1/6 이다.
    print("\n곱 규칙은 깨진다:")
    print(f"E[XY] = {(X * Y).mean():.4f} (theoretical: {1/4:.4f})")
    print(f"E[X]E[Y] = {X.mean() * Y.mean():.4f} (theoretical: {1/6:.4f})")
    print(f"Cov(X,Y) = {np.cov(X, Y, ddof=0)[0, 1]:.4f} (theoretical: {1/12:.4f})")

    # 분산의 덧셈도 깨진다. 교차항 2Cov 가 남는다.
    print("\n분산의 덧셈도 깨진다:")
    vx, vy = X.var(), Y.var()
    print(f"Var(X) = {vx:.4f} (theoretical: {1/12:.4f})")
    print(f"Var(Y) = {vy:.4f} (theoretical: {4/45:.4f})")
    print(f"Var(X+Y) = {(X + Y).var():.4f} (theoretical: {61/180:.4f})")
    print(f"Var(X)+Var(Y) = {vx + vy:.4f} (theoretical: {1/12 + 4/45:.4f})")
    print(f"rho = {np.corrcoef(X, Y)[0, 1]:.4f} (theoretical: {np.sqrt(15)/4:.4f})")

linearity_demonstration()

출력:

X and Y = X² are dependent, but linearity still holds:
E[X] = 0.4995 (theoretical: 0.5)
E[Y] = 0.3326 (theoretical: 0.3333)
E[X + Y] = 0.8321
E[X] + E[Y] = 0.8321

곱 규칙은 깨진다:
E[XY] = 0.2492 (theoretical: 0.2500)
E[X]E[Y] = 0.1661 (theoretical: 0.1667)
Cov(X,Y) = 0.0830 (theoretical: 0.0833)

분산의 덧셈도 깨진다:
Var(X) = 0.0831 (theoretical: 0.0833)
Var(Y) = 0.0885 (theoretical: 0.0889)
Var(X+Y) = 0.3377 (theoretical: 0.3389)
Var(X)+Var(Y) = 0.1716 (theoretical: 0.1722)
rho = 0.9683 (theoretical: 0.9682)

선형성 쪽은 \(E[X+Y] = 0.8321\) 과 \(E[X] + E[Y] = 0.8321\) 이 소수 넷째 자리까지 한 글자도 다르지 않다. 모의실험의 흔들림조차 양쪽에 똑같이 들어가기 때문이며, 이 등식은 표본마다 성립하는 항등식이다. 이론값 \(0.8333\) 과는 \(0.0012\) 차이인데, \(\operatorname{sd}(X+Y) = \sqrt{0.3389} = 0.582\) 이므로 표준오차가 \(0.582/\sqrt{100000} = 0.0018\) 이고 벗어남은 그 안이다.

나머지 값들도 유도와 맞는다. \(E[XY] = 0.2492\) 대 \(0.25\), 공분산 \(0.0830\) 대 \(1/12 = 0.0833\), \(\operatorname{Var}(X+Y) = 0.3377\) 대 \(61/180 = 0.3389\), \(\rho = 0.9683\) 대 \(\sqrt{15}/4 = 0.9682\) 다.

정리하면 이렇다. 같은 \((X, Y)\) 한 쌍에서 선형성은 정확히 성립하고 곱 규칙은 \(1.5\) 배, 분산의 덧셈은 \(1.97\) 배 빗나간다. 독립이 필요한 자리와 필요 없는 자리가 이렇게 갈린다.

\(g(X) = X^2\)에 대해 LOTUS를 쓰면

이다.

기댓값의 선형성에 의해

이다. 각 \(X_i\)는 \(E[X_i] = P(X_i = 1) = 0.5\)인 베르누이 확률변수이므로

이다. 100번 던지면 앞면이 50번 나올 것으로 기대한다. 중요한 점은 선형성이 던지기의 독립 여부와 무관하게 성립한다는 것이다. 던지기가 종속이더라도 같은 답이 나온다.

LOTUS에 의해 \(E[Y] = E[(X-3.5)^2] = \sum_{x=1}^{6} (x - 3.5)^2 \cdot P(X=x)\)이다. 각 값에 대해 \(P(X=x) = 1/6\)이므로

이다. 참고로 \(E[X] = 3.5\)이므로 이는 정확히 공정한 주사위의 \(\text{Var}(X)\)다.

\(X\)와 \(Y\)가 독립이므로

이다. \(E[(X+Y)^2]\)의 경우 제곱을 전개하면

이다.

연속인 경우: 푸비니 정리에 의해

이며, 순서 교환은 비음성에 의해 정당화된다.

정숫값인 경우:

용도: 꼬리합 공식을 쓰면 생존확률로부터 기댓값을 계산할 수 있다(표준적인 확률밀도함수 적분보다 쉬울 때가 있다). 예: 첫 성공까지의 시행 횟수를 세는 기하확률변수에서 \(P(X > n) = (1 - p)^n\)이므로 \(\mathbb{E}[X] = \sum_{n=0}^\infty (1 - p)^n = 1/p\)이다.

정의에 의해 \(g(y) = \mathbb{E}[X \mid Y = y] = \int x f_{X \mid Y}(x \mid y) dx\)이다.

\(\mathbb{E}[g(Y)] = \int g(y) f_Y(y) dy = \int \int x f_{X \mid Y}(x \mid y) f_Y(y) dx \, dy = \int \int x f_{X, Y}(x, y) dx \, dy = \int x f_X(x) dx = \mathbb{E}[X]\).

\(\square\)

확률, 통계, 그리고 기댓값에 대한 동적계획법 접근 전반에서 쓰인다.

• 최적 예측자로서의 조건부 기댓값: \(\mathbb{E}[X \mid Y]\)는 모든 함수 \(g\)에 대해 \(\mathbb{E}[(X - g(Y))^2]\)을 최소화한다.
• 탑 성질(반복 기댓값): \(\mathbb{E}[X] = \mathbb{E}[\mathbb{E}[X \mid Y]] = \mathbb{E}[\mathbb{E}[\mathbb{E}[X \mid Y, Z] \mid Y]]\) 등.
• MCMC / 분산 감소: (가능할 때) \(X\)를 \(\mathbb{E}[X \mid Y]\)로 대체하면 라오–블랙웰 정리를 통해 추정량의 분산이 줄어든다.
• 강화학습: 벨만 방정식 \(V(s) = \mathbb{E}[R + \gamma V(s') \mid s]\)이 반복된 조건부 기댓값이다.

\(\mathbb{E}[X]\)는 \(\mathbb{E}[X^+]\)와 \(\mathbb{E}[X^-]\)로 나누어 정의한다(\(X^\pm\)는 양·음 부분).

코시분포가 세 번째 경우다. 밀도가 \(f(x) = 1/(\pi(1+x^2))\)이라 양쪽 꼬리 모두 \(x f(x) \sim 1/(\pi x)\)로 발산한다.

import numpy as np

rng = np.random.default_rng(0)
print("코시 표본평균의 궤적 — 세 번 반복")
for trial in range(3):
    x = rng.standard_cauchy(1_000_000)
    run = np.cumsum(x) / np.arange(1, len(x) + 1)
    print(f"  시행 {trial+1}: n=10^3 {run[999]:>9.3f}   10^4 {run[9999]:>9.3f}"
          f"   10^5 {run[99999]:>9.3f}   10^6 {run[-1]:>9.3f}")

출력:

코시 표본평균의 궤적 — 세 번 반복
  시행 1: n=10^3    -0.934   10^4    -0.144   10^5     0.297   10^6     0.471
  시행 2: n=10^3     0.093   10^4    -1.466   10^5    -0.554   10^6    -0.884
  시행 3: n=10^3    -1.187   10^4     0.217   10^5    -0.966   10^6    -2.630

표본평균이 정착하지 않는다. \(n\)이 \(100\)만이어도 시행마다 \(0.471\), \(-0.884\), \(-2.630\)으로 흩어진다. 큰수의 법칙이 적용되지 않기 때문이다.

더 놀라운 사실. 코시분포의 표본평균 \(\bar{X}_n\)은 \(n\)과 무관하게 원래 코시분포와 정확히 같은 분포를 갖는다. 관측을 \(100\)만 개 모으는 것이 하나만 보는 것과 똑같다.

"무한"과 "정의되지 않음"의 실무적 차이.

• \(\mathbb{E}[X] = +\infty\)(예: 파레토 \(\alpha \le 1\), 상트페테르부르크 게임)이면 표본평균이 \(\infty\)로 발산한다. 방향이 정해져 있다.
• 정의되지 않으면(코시) 표본평균이 어느 방향으로도 가지 않고 헤맨다. 어떤 값에도 수렴하지 않는다.

어디서 코시를 만나는가. 두 독립 정규의 비 \(Z_1/Z_2\)가 코시이며, 그래서 비를 다룰 때 조심해야 한다. 회귀계수의 비, 두 추정값의 비, 각도 측정에서 코시가 자연스럽게 나타난다.

처방. 평균이 없으면 중앙값을 쓴다. 코시분포의 중앙값은 잘 정의되고, 표본중앙값이 그것으로 수렴한다(2장 평균·중앙값 문서 연습문제 8에서 본 대로 코시에서 중앙값이 압도적으로 낫다). \(\square\)

정리. 모든 (가측) 함수 \(g\)에 대해

이다. 즉 \(\mathbb{E}[Y\mid X]\)가 제곱오차를 최소화하는 \(X\)의 함수다.

증명. \(m(X) = \mathbb{E}[Y\mid X]\)라 두고 전개하면

이다. 가운데 항은 \(X\)로 조건을 걸어 계산하면 \(\mathbb{E}[Y-m \mid X] = 0\)이므로 사라진다. 남는 것이 \(\mathbb{E}[(m-g)^2] \ge 0\)이다. \(\square\)

import numpy as np

rng = np.random.default_rng(0)
n = 400_000
x = rng.normal(0, 1, n)
y = x ** 2 + rng.normal(0, 0.5, n)             # 참 관계는 이차

predictors = {
    "E[Y|X] = X^2": x ** 2,
    "최적 선형 예측": np.polyval(np.polyfit(x, y, 1), x),
    "상수 (전체 평균)": np.full(n, y.mean()),
}
for label, pred in predictors.items():
    print(f"  {label:>16}: MSE {np.mean((y - pred) ** 2):.5f}")

print(f"\n잡음의 분산 = {0.5 ** 2:.5f}  ← E[Y|X] 의 MSE 가 여기에 닿는다")
print(f"Cov(X, X^2) = {np.cov(x, x ** 2)[0,1]:+.5f}  ← 그래서 선형 예측이 상수와 같다")

출력:

      E[Y|X] = X^2: MSE 0.24974
          최적 선형 예측: MSE 2.25532
        상수 (전체 평균): MSE 2.25536

잡음의 분산 = 0.25000  ← E[Y|X] 의 MSE 가 여기에 닿는다
Cov(X, X^2) = -0.00444  ← 그래서 선형 예측이 상수와 같다

\(\mathbb{E}[Y\mid X]\)의 MSE가 \(0.2497\)로 잡음 분산 \(0.25\)에 닿는다. 더 줄일 수 없는 하한이다.

최적 선형 예측이 상수 예측과 거의 똑같다(\(2.2553\) 대 \(2.2554\)). \(\operatorname{Cov}(X, X^2) = 0\)이므로 선형 회귀의 기울기가 \(0\)이 되기 때문이다. 이 절의 독립성과 무상관성의 차이에서 다룬 "무상관이지만 종속"이 예측의 언어로 나타난 것이다.

기하학적 해석. \(L^2\) 공간에서 \(\mathbb{E}[Y\mid X]\)는 \(Y\)를 "\(X\)의 함수들이 이루는 부분공간" 위로 직교사영한 것이다. 잔차 \(Y - \mathbb{E}[Y\mid X]\)가 그 부분공간의 모든 원소와 직교하며, 그것이 위 증명의 가운데 항이 사라지는 이유다.

0장의 최소제곱과 같은 구조다. 거기서는 \(\mathbf{y}\)를 \(\mathbf{X}\)의 열공간에 사영했고, 여기서는 \(Y\)를 \(X\)의 함수공간에 사영한다. 선형회귀는 함수공간을 선형함수로 제한한 특수한 경우이며, 그래서 위 예에서 참 관계를 놓친다.

함의. 회귀분석의 목표가 \(\mathbb{E}[Y\mid X]\)를 추정하는 것이며, 모형의 유연성이 곧 함수공간을 얼마나 넓게 잡느냐다. 1장에서 본 편향–분산 절충이 이 선택의 다른 이름이다. \(\square\)

\(N\)이 확률변수이고 \(X_1, X_2, \ldots\)가 i.i.d.이며 \(N\)과 독립일 때

이다. 분산은 조금 더 복잡하다.

증명은 전기댓값의 법칙이다. \(N\)으로 조건을 걸면 \(\mathbb{E}[S \mid N=n] = n\mathbb{E}[X]\)이므로

이다. 분산도 총분산의 법칙으로 같은 방식으로 나온다. \(\square\)

import numpy as np

rng = np.random.default_rng(0)
lam, mu, sigma = 4.0, 3.0, 1.0

N = rng.poisson(lam, 50_000)
S = np.array([rng.normal(mu, sigma, k).sum() if k else 0.0 for k in N])

print(f"E[S]   모의 {S.mean():.4f}   이론 E[N]E[X] = {lam * mu:.4f}")
print(f"Var(S) 모의 {S.var():.4f}   이론 E[N]Var(X) + Var(N)E[X]^2 = "
      f"{lam * sigma ** 2 + lam * mu ** 2:.4f}")

출력:

E[S]   모의 12.0062   이론 E[N]E[X] = 12.0000
Var(S) 모의 40.1686   이론 E[N]Var(X) + Var(N)E[X]^2 = 40.0000

두 공식이 모두 맞는다.

분산 공식의 두 항이 각각 무엇인가.

• \(\mathbb{E}[N]\operatorname{Var}(X)\): 항의 개수가 고정이었다면 있었을 변동
• \(\operatorname{Var}(N)\mathbb{E}[X]^2\): 개수가 흔들려서 추가로 생기는 변동

포아송이면 \(\operatorname{Var}(N) = \mathbb{E}[N]\)이라 둘째 항이 \(\lambda\mu^2\)로 커진다. 위 예에서 \(4 + 36 = 40\) 중 \(36\)이 개수의 변동에서 온다. 개수의 불확실성이 지배적이다.

어디에 쓰이는가.

복합 포아송 모형이 이 구조의 표준 이름이며, 보험 수리의 기본 도구다.

주의: \(N\)과 \(X_i\)가 독립이어야 한다. 독립이 아니면 공식이 깨진다. 순차 검정처럼 \(N\)이 관측에 의존해 정해지는 경우에는 정지시각이라는 조건이 필요하며, 그것이 왈드 항등식의 원래 형태다. \(\square\)

쿠폰 수집가. \(n\)종의 쿠폰을 모두 모으려면 몇 번 뽑아야 하는가?

\(T_k\)를 "\(k-1\)종을 모은 뒤 새 종류가 나올 때까지의 뽑기 수"라 하면 \(T_k \sim \text{Geometric}(p_k)\)이고 \(p_k = (n-k+1)/n\)이다. 선형성으로

import numpy as np

rng = np.random.default_rng(0)
print("쿠폰 수집가")
print(f"{'n':>5}{'모의':>11}{'이론 n H_n':>13}{'근사 n(ln n + 0.577)':>22}")
for n in (5, 10, 50):
    H = sum(1 / i for i in range(1, n + 1))
    trials = []
    for _ in range(20_000):
        seen, count = set(), 0
        while len(seen) < n:
            seen.add(rng.integers(n)); count += 1
        trials.append(count)
    print(f"{n:>5}{np.mean(trials):>11.3f}{n * H:>13.3f}"
          f"{n * (np.log(n) + 0.5772):>22.3f}")

출력:

쿠폰 수집가
    n         모의     이론 n H_n    근사 n(ln n + 0.577)
    5     11.377       11.417                10.933
   10     29.244       29.290                28.798
   50    225.207      224.960               224.461

모의와 이론이 맞는다. 마지막 한 종류를 얻는 데만 평균 \(n\)번이 걸리므로, 전체의 상당 부분이 막바지에 소요된다.

매칭 문제. \(n\)명이 모자를 무작위로 다시 가져갈 때 자기 모자를 받는 사람의 기대 수는?

\(I_i\)를 "\(i\)번이 자기 모자를 받음"의 지시변수라 하면 \(\mathbb{E}[I_i] = 1/n\)이므로

\(n\)과 무관하게 항상 \(1\)이다.

import numpy as np

rng = np.random.default_rng(1)
print("\n매칭 문제 (자기 모자를 받는 사람 수)")
print(f"{'n':>5}{'평균':>10}{'분산':>10}{'0명일 확률':>13}{'1/e':>9}")
for n in (5, 20, 100):
    counts = np.array([np.sum(rng.permutation(n) == np.arange(n))
                       for _ in range(100_000)])
    print(f"{n:>5}{counts.mean():>10.4f}{counts.var():>10.4f}"
          f"{np.mean(counts == 0):>13.4f}{1 / np.e:>9.4f}")

출력:

매칭 문제 (자기 모자를 받는 사람 수)
    n        평균        분산       0명일 확률      1/e
    5    1.0005    1.0008       0.3667   0.3679
   20    1.0055    1.0046       0.3653   0.3679
  100    1.0012    1.0016       0.3679   0.3679

평균과 분산이 모두 \(1\)에 가깝고 \(n\)과 무관하다. 그리고 아무도 자기 모자를 받지 못할 확률이 \(1/e \approx 0.368\)로 수렴한다.

지시변수 기법의 힘은 독립을 요구하지 않는다는 것이다. 매칭 문제에서 \(I_i\)들은 명백히 종속이다(\(n-1\)명이 자기 모자를 받으면 나머지 한 명도 반드시 받는다). 그런데 기댓값의 선형성은 독립과 무관하게 성립하므로 합의 기댓값을 그대로 계산할 수 있다.

이것이 이 절의 핵심 도구다. 복잡한 확률변수를 지시변수의 합으로 쪼개면, 각 지시변수의 기댓값은 단순한 확률 하나이고 선형성이 나머지를 해 준다. 분산은 종속성 때문에 더 어렵지만(공분산 항이 필요하다), 기댓값만큼은 언제나 쉽다. \(\square\)