적률생성함수¶

이다. 두 번 미분해 보면

이고 \(\text{Var}(X) = (\sigma^2 + \mu^2) - \mu^2 = \sigma^2\)이다. 적분 한 번 하지 않고 평균과 분산을 얻었다.

(1) 해석적으로. \(M\) 을 \(0\) 근방에서 전개한다. 정리 1 의 증명이 보인 대로 \(M\) 은 \(0\) 근방에서 해석적이므로 전개가 가능하다.

1계. 빼면 짝수 차수 항이 모두 지워진다.

\(2h\) 로 나누면

2계. 더하면 홀수 차수 항이 지워진다.

\(h^2\) 으로 나누면

둘 다 절단오차가 \(h^2\) 에 비례한다. 그런데 계수에 \(M\) 의 고계도함수가 들어 있고, 정리 1 에 따르면 \(M^{(n)}(0) = E[X^n]\) 이므로 그것은 곧 적률이다. \(N(3,4)\) 에서는

이므로 절단오차의 예측값은

이다. 수치미분의 오차를 예측하는 데 그 자신이 꺼내려는 적률이 쓰인다는 점이 재미있다.

(2) 반올림이 반대 방향으로 민다. \(h\) 를 줄이면 절단오차는 \(h^2\) 으로 줄지만 분자가 거의 상쇄되어 유효숫자를 잃는다. 배정밀도에서 \(M(\pm h)\) 는 상대오차 \(\varepsilon \approx 2.2 \times 10^{-16}\) 을 달고 다니므로 분자에 생기는 절대오차는 \(h\) 와 무관하게 \(\varepsilon M(0)\) 규모다. 그것을 무엇으로 나누느냐가 두 차분을 가른다.

\(h = 10^{-6}\) 이면 앞의 것은 \(\varepsilon\) 을 \(10^{6}\) 배 키우지만 뒤의 것은 \(10^{12}\) 배 키운다. 같은 분자 오차를 \(h^2\) 으로 나누기 때문에 2계 차분이 훨씬 먼저 무너진다. 그래서 최적의 \(h\) 가 2계 쪽에서 더 크고, 코드가 1계에 \(10^{-4}\), 2계에 \(10^{-3}\) 을 쓴 이유가 이것이다.

(3) 수치적으로. 두 예측을 실제 오차와 나란히 잰다.

import numpy as np

def mgf_normal(t, mu, sigma2):
    """정규분포 N(mu, sigma2)의 적률생성함수.

        M(t) = exp(mu*t + sigma2*t^2/2)
    """
    return np.exp(mu * t + sigma2 * t**2 / 2)


def d1(f, x, h=1e-4):
    """중심차분으로 1계 도함수를 근사한다:  (f(x+h) - f(x-h)) / (2h)"""
    return (f(x + h) - f(x - h)) / (2 * h)


def d2(f, x, h=1e-3):
    """중심차분으로 2계 도함수를 근사한다:  (f(x+h) - 2f(x) + f(x-h)) / h^2

    h를 1계보다 크게 잡은 이유가 있다. 2계 차분은 h^2 으로 나누므로
    h가 너무 작으면 분자의 반올림 오차가 크게 증폭된다.
    """
    return (f(x + h) - 2 * f(x) + f(x - h)) / h**2


# "적률생성함수"라는 이름 그대로, t=0 에서 k번 미분하면 k번째 적률이 나온다.
#   M'(0)  = E[X]
#   M''(0) = E[X^2]
mu, sigma2 = 3.0, 4.0
M = lambda t: mgf_normal(t, mu, sigma2)

E_X = d1(M, 0.0)
E_X2 = d2(M, 0.0)
Var_X = E_X2 - E_X**2

print(f"E[X] = {E_X:.4f} (theoretical: {mu})")
print(f"E[X²] = {E_X2:.4f} (theoretical: {sigma2 + mu**2})")
print(f"Var(X) = {Var_X:.4f} (theoretical: {sigma2})")

# 유도한 절단오차 예측과 실제 오차를 h마다 나란히 잰다.
#   1계:  h^2 * M'''(0) / 6,    M'''(0)  = E[X^3]
#   2계:  h^2 * M''''(0) / 12,  M''''(0) = E[X^4]
m3 = mu**3 + 3 * mu * sigma2
m4 = mu**4 + 6 * mu**2 * sigma2 + 3 * sigma2**2
print(f"\nE[X^3] = {m3:.0f},  E[X^4] = {m4:.0f}")
print(f"{'h':>7}{'d1 err':>12}{'pred':>12}{'d2 err':>12}{'pred':>12}")
for k in range(1, 10):
    h = 10.0 ** (-k)
    e1 = abs(d1(M, 0.0, h) - mu)
    e2 = abs(d2(M, 0.0, h) - (sigma2 + mu**2))
    print(f"{h:>7.0e}{e1:>12.3e}{h**2 * m3 / 6:>12.3e}"
          f"{e2:>12.3e}{h**2 * m4 / 12:>12.3e}")

출력:

E[X] = 3.0000 (theoretical: 3.0)
E[X²] = 13.0000 (theoretical: 13.0)
Var(X) = 4.0000 (theoretical: 4.0)

E[X^3] = 63,  E[X^4] = 345
      h      d1 err        pred      d2 err        pred
  1e-01   1.067e-01   1.050e-01   2.912e-01   2.875e-01
  1e-02   1.050e-03   1.050e-03   2.875e-03   2.875e-03
  1e-03   1.050e-05   1.050e-05   2.875e-05   2.875e-05
  1e-04   1.050e-07   1.050e-07   2.874e-07   2.875e-07
  1e-05   1.047e-09   1.050e-09   3.460e-08   2.875e-09
  1e-06   8.627e-11   1.050e-11   4.548e-05   2.875e-11
  1e-07   8.627e-11   1.050e-13   7.116e-04   2.875e-13
  1e-08   3.972e-09   1.050e-15   3.227e-01   2.875e-15
  1e-09   2.618e-08   1.050e-17   1.300e+01   2.875e-17

절단오차 예측이 그대로 맞는다. \(h = 10^{-2}\) 부터 \(10^{-4}\) 까지 "실제 오차" 열과 "예측" 열이 유효숫자 세 자리까지 같고, \(h\) 를 열 배 줄일 때마다 오차가 백 분의 일로 줄어 \(h^2\) 비례도 확인된다. \(h = 10^{-1}\) 에서만 \(1.067\) 대 \(1.050\) 으로 조금 어긋나는데, 그때는 버린 \(O(h^4)\) 항이 아직 보일 만큼 크기 때문이다.

무너지는 자리가 다르다. 1계 차분은 \(h = 10^{-6}\) 에서 오차 \(8.6\times 10^{-11}\) 로 바닥을 치고 그 아래로 내려가도 \(10^{-8}\) 수준에 머문다. 2계 차분은 \(h = 10^{-5}\) 에서 \(3.5\times 10^{-8}\) 로 바닥을 친 뒤 급격히 올라가 \(h = 10^{-9}\) 에서는 오차가 \(13.0\), 곧 답 자체보다 크다. 그 자리에서 \(D_2\) 가 돌려준 값은 \(0\) 근처의 잡음일 뿐이다.

코드가 고른 \(h\) 는 둘 다 절단오차가 지배하는 안전한 쪽에 있다. \(h = 10^{-4}\) 의 1계 오차가 \(1.05\times 10^{-7}\) 이고 \(h = 10^{-3}\) 의 2계 오차가 \(2.875\times 10^{-5}\) 이므로, 소수 넷째 자리까지 찍는 출력에서는 둘 다 보이지 않고 \(3.0000\), \(13.0000\), \(4.0000\) 이 깨끗하게 나왔다. 2계 쪽을 \(h = 10^{-6}\) 으로 "더 정밀하게" 바꾸면 오차가 \(4.5\times 10^{-5}\) 로 오히려 커진다. 수치미분에서 \(h\) 는 작을수록 좋은 것이 아니다.

이다. 이것은 \(N(\mu_1+\mu_2,\ \sigma_1^2+\sigma_2^2)\)의 적률생성함수이므로, 유일성(정리 2)에 의해

이다.

(1) 둘 다 아무 분포에나 성립한다. 먼저 \(t = 0\) 을 넣으면 \(e^{0 \cdot X} = 1\) 이므로

이다. 분포가 무엇이든 그렇다. 볼록성도 마찬가지다. 기댓값 안으로 미분을 넣으면

인데 \(X^2 \ge 0\) 이고 \(e^{tX} > 0\) 이므로 피적분함수가 음이 될 수 없다. 적률생성함수는 유한한 곳에서 언제나 아래로 볼록하다. \(X\) 가 상수가 아니면 부등호가 엄격해 진짜로 휜다.

기울기는 정리 1 이 곧바로 준다. \(M_X'(0) = E[X]\) 이므로

이다. 네 곡선이 한 점에서 만난 뒤 갈라지는 각도가 곧 평균이다. 특히 \(N(0,1)\) 만 평균이 \(0\) 이라 \(t = 0\) 에서 접선이 수평이고, 볼록하므로 그 점이 곡선의 최솟값이다. 그림에서 파란 곡선만 \(U\) 자로 보이는 이유가 이것이다.

(2) 왼쪽 극한은 원점의 확률질량이다. \(X \ge 0\) 이면 \(t \to -\infty\) 일 때 \(e^{tX}\) 는 \(X = 0\) 에서 \(1\), \(X > 0\) 에서 \(0\) 으로 간다. \(0 \le e^{tX} \le 1\) 이므로 유계수렴정리로 극한을 기댓값 안에 넣을 수 있고

이다. 셋을 대입하면 \(\text{Bernoulli}(0.5) \to 1 - p = 0.5\), \(\text{Poisson}(3) \to e^{-3} = 0.0498\), \(\text{Exp}(1) \to 0\) 이다. 마지막 것이 \(0\) 인 이유는 지수분포가 연속이라 한 점에 질량이 없기 때문이다. \(N(0,1)\) 은 다르다. 음의 값을 가질 수 있으므로 \(t\) 가 음으로 크면 \(e^{tX}\) 가 왼쪽 꼬리에서 터져 \(M(t) = e^{t^2/2} \to \infty\) 가 된다.

\(\text{Exp}(1)\) 이 끊기는 까닭. 정의대로 적분하면

이고 이것은 \(1 - t > 0\), 곧 \(t < 1\) 일 때만 수렴해 \(1/(1-t)\) 를 준다. \(t \ge 1\) 에서는 \(e^{tx}\) 가 꼬리 \(e^{-x}\) 보다 빨리 커져 적분이 발산한다. 꼬리가 지수적으로 줄어드는 속도가 정의역의 오른쪽 끝을 정한다.

그래도 정리 2 는 그대로 쓸 수 있다. 정의역 \((-\infty, 1)\) 이 \(0\) 을 품은 열린구간이기 때문이다. 단서가 걸러 내려는 것은 이런 분포가 아니라 \(0\) 의 어느 근방에서도 적률생성함수가 무한한 분포다. 연습문제 5 의 코시분포는 \(t \ne 0\) 인 모든 \(t\) 에서, 연습문제 6 의 로그정규분포는 모든 \(t > 0\) 에서 적분이 발산하며, 그래서 로그정규는 적률이 모두 같은 다른 분포를 곁에 둘 수 있다.

(3) 수치적으로. 네 곡선을 그리고 그림에서 읽을 값들을 따로 찍어 둔다.

import numpy as np
import matplotlib.pyplot as plt

def plot_mgf_comparison():
    """분포 넷의 적률생성함수를 한 그림에 겹쳐 그린다.

    네 곡선이 모두 t=0 에서 1을 지난다는 점에 주목하라.
    M(0) = E[e^0] = E[1] = 1 이므로 어떤 분포든 반드시 그렇다.
    그리고 t=0 에서의 기울기가 곧 그 분포의 평균이다.
    """
    t = np.linspace(-1.5, 1.5, 300)

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

    # 정규분포 N(0,1):  M(t) = exp(t^2/2). 모든 t에서 유한하다.
    ax.plot(t, np.exp(t**2 / 2), label='N(0, 1)', lw=2)

    # 지수분포 Exp(1):  M(t) = 1/(1-t).  t < 1 에서만 존재한다.
    # 꼬리가 지수적으로 감소하는 속도보다 e^{tx} 가 빨리 커지면 적분이 발산하기 때문이다.
    # 그래서 t >= 1 인 부분을 아예 잘라 내고 그린다.
    t_exp = t[t < 1]
    ax.plot(t_exp, 1 / (1 - t_exp), label='Exp(1)', lw=2)

    # 포아송(3):  M(t) = exp(lambda*(e^t - 1))
    lam = 3
    ax.plot(t, np.exp(lam * (np.exp(t) - 1)), label='Poisson(3)', lw=2)

    # 베르누이(0.5):  M(t) = 1 - p + p*e^t.  값이 둘뿐이라 가장 단순한 형태다.
    p = 0.5
    ax.plot(t, 1 - p + p * np.exp(t), label='Bernoulli(0.5)', lw=2)

    ax.set_xlabel('t')
    ax.set_ylabel('M_X(t)')
    ax.set_title('Moment Generating Functions')
    ax.set_ylim(0, 15)
    ax.legend()
    ax.spines[['top', 'right']].set_visible(False)
    plt.tight_layout()
    plt.show()

plot_mgf_comparison()

# 그림에서 읽을 수치를 따로 찍어 둔다.
mgfs = {
    "N(0,1)":        lambda t: np.exp(t ** 2 / 2),
    "Exp(1)":        lambda t: 1 / (1 - t),
    "Poisson(3)":    lambda t: np.exp(3 * (np.exp(t) - 1)),
    "Bernoulli(.5)": lambda t: 0.5 + 0.5 * np.exp(t),
}
h = 1e-5
print(f"{'dist':>14}{'M(0)':>9}{'slope(0)':>10}{'M(-1.5)':>10}{'M(0.5)':>10}")
for name, M in mgfs.items():
    slope = (M(h) - M(-h)) / (2 * h)
    print(f"{name:>14}{M(0.0):>9.4f}{slope:>10.4f}{M(-1.5):>10.4f}{M(0.5):>10.4f}")

# t -> -무한 에서의 극한. X >= 0 이면 e^{tX} 가 1{X=0} 으로 가므로
# 극한이 원점의 확률질량 P(X=0) 이다. N(0,1) 은 음의 값을 가지므로 발산한다.
print("\nt -> -inf 에서의 극한  (X >= 0 이면 P(X=0))")
atom = {"Exp(1)": 0.0, "Poisson(3)": np.exp(-3), "Bernoulli(.5)": 0.5}
for name, a in atom.items():
    M = mgfs[name]
    print(f"{name:>14}  M(-5)={M(-5.0):.6f}  M(-50)={M(-50.0):.6f}  P(X=0)={a:.6f}")
print(f"{'N(0,1)':>14}  M(-5)={mgfs['N(0,1)'](-5.0):.1f}  -> 발산 (X 가 음수를 가질 수 있다)")

# 각 곡선이 M = 15 를 넘어 화면 밖으로 나가는 자리
print(f"\nM=15 를 넘는 t:  Poisson(3) {np.log(1 + np.log(15) / 3):.4f}"
      f"   Exp(1) {1 - 1 / 15:.4f}"
      f"   N(0,1) {np.sqrt(2 * np.log(15)):.4f}   Bernoulli(.5) {np.log(29):.4f}")

출력:

          dist     M(0)  slope(0)   M(-1.5)    M(0.5)
        N(0,1)   1.0000    0.0000    3.0802    1.1331
        Exp(1)   1.0000    1.0000    0.4000    2.0000
    Poisson(3)   1.0000    3.0000    0.0972    7.0018
 Bernoulli(.5)   1.0000    0.5000    0.6116    1.3244

t -> -inf 에서의 극한  (X >= 0 이면 P(X=0))
        Exp(1)  M(-5)=0.166667  M(-50)=0.019608  P(X=0)=0.000000
    Poisson(3)  M(-5)=0.050804  M(-50)=0.049787  P(X=0)=0.049787
 Bernoulli(.5)  M(-5)=0.503369  M(-50)=0.500000  P(X=0)=0.500000
        N(0,1)  M(-5)=268337.3  -> 발산 (X 가 음수를 가질 수 있다)

M=15 를 넘는 t:  Poisson(3) 0.6433   Exp(1) 0.9333   N(0,1) 2.3273   Bernoulli(.5) 3.3673

네 값이 모두 \(M(0) = 1\) 로 나오고 기울기가 \(0\), \(1\), \(3\), \(0.5\) 로 (1) 의 평균과 맞는다. 왼쪽 극한도 맞는다. \(t = -50\) 에서 \(\text{Poisson}(3)\) 은 \(0.049787 = e^{-3}\) 에 이미 소수 여섯째 자리까지 닿았고 \(\text{Bernoulli}(0.5)\) 는 정확히 \(0.5\) 다. \(\text{Exp}(1)\) 만 \(0.019608 = 1/51\) 로 아직 \(0\) 에서 멀다. 극한이 \(0\) 인 것은 맞지만 \(1/(1-t)\) 는 \(1/\lvert t \rvert\) 꼴로 천천히 내려가기 때문이며, 지수적으로 내려가는 앞의 둘과 속도가 다르다.

그림에서 읽을 것은 왼쪽과 오른쪽에서 순서가 뒤집힌다는 점이다. 오른쪽 \(t = 0.5\) 에서는 \(\text{Poisson}(3)\) 이 \(7.00\) 으로 가장 높고 \(N(0,1)\) 이 \(1.13\) 으로 가장 낮다. 왼쪽 \(t = -1.5\) 에서는 정반대로 \(N(0,1)\) 이 \(3.08\) 로 가장 높고 \(\text{Poisson}(3)\) 이 \(0.097\) 로 가장 낮다. 평균이 클수록 오른쪽에서 빨리 솟지만, 왼쪽에서는 원점에 질량이 얼마나 있는가가 높이를 정하기 때문이다.

화면 밖으로 나가는 자리도 계산과 맞는다. \(\text{Poisson}(3)\) 은 \(t = 0.6433\) 에서, \(\text{Exp}(1)\) 은 \(t = 0.9333\) 에서 \(M = 15\) 를 넘어 사라진다. \(N(0,1)\) 과 \(\text{Bernoulli}(0.5)\) 는 각각 \(t = 2.33\), \(t = 3.37\) 에서야 넘으므로 그림 안에 끝까지 남는다. \(\text{Exp}(1)\) 이 \(t = 1\) 직전에서 수직으로 서는 것과 \(\text{Poisson}(3)\) 이 그보다 먼저 사라지는 것은 성격이 다르다. 앞의 것은 정의역의 끝이고, 뒤의 것은 \(y\) 축 범위를 \(15\) 로 잘랐기 때문일 뿐이다. 포아송의 적률생성함수는 모든 \(t\) 에서 유한하다.

(1) 해석적으로. 표에서 정규분포의 적률생성함수를 가져와 정리 3 대로 곱한다.

지수의 어깨가 그냥 더해진다. 마지막 식이 \(N(7, 25)\) 의 적률생성함수이므로 유일성(정리 2)에 의해 \(S \sim N(7, 25)\) 이고 \(\operatorname{sd}(S) = 5\) 다. 보기 3 의 일반식에 \(\mu_1 = 2\), \(\mu_2 = 5\), \(\sigma_1 = 3\), \(\sigma_2 = 4\) 를 넣은 것과 같다.

더해지는 것은 분산이지 표준편차가 아니다. \(t^2\) 의 계수에 \(\sigma^2\) 이 들어 있기 때문이다. 그래서

이고, 이 관계는 \(3\)–\(4\)–\(5\) 직각삼각형과 같은 모양이다. 표준편차는 피타고라스로 합쳐진다.

\(7\) 이 되려면 독립이 깨져야 한다. 일반적으로

이므로 \(\rho = 1\) 일 때 \(9 + 16 + 24 = 49\), 곧 \(\operatorname{sd} = 7\) 이다. \(3 + 4 = 7\) 은 두 변수가 완전히 같이 움직일 때만 맞는 셈이며, 독립일 때 그보다 작은 것이 분산이 주는 이득이다.

(2) 이론값. 표본 \(n = 10^5\) 개로 재면 두 추정값 모두 흔들린다. 평균 쪽은

이고, 정규표본의 표본분산은 \(\operatorname{Var}(\hat\sigma^2) \approx 2\sigma^4/n\) 이므로

이다. 평균은 \(7 \pm 0.03\), 분산은 \(25 \pm 0.22\) 안쪽이면(두 표준오차) 아무 일도 아니다.

(3) 수치적으로.

import numpy as np

def verify_sum_of_normals(n_simulations=100_000):
    """독립인 정규분포의 합이 다시 정규분포임을 모의실험으로 확인한다.

    적률생성함수로 증명한 결과를 눈으로 확인하는 것이다.
    MGF가 곱해지면 지수의 어깨가 더해지므로
    평균은 mu1+mu2, 분산은 sigma1^2+sigma2^2 이 된다.
    """
    np.random.seed(42)
    mu1, sigma1 = 2, 3
    mu2, sigma2 = 5, 4

    X1 = np.random.normal(mu1, sigma1, n_simulations)
    X2 = np.random.normal(mu2, sigma2, n_simulations)
    S = X1 + X2

    # 분산은 더해지지만 **표준편차는 더해지지 않는다**는 점에 주의하라.
    # sqrt(3^2 + 4^2) = 5 이지 3 + 4 = 7 이 아니다.

    print(f"E[X1+X2] = {S.mean():.4f} (theoretical: {mu1 + mu2})")
    print(f"Var(X1+X2) = {S.var():.4f} (theoretical: {sigma1**2 + sigma2**2})")

    # 몬테카를로 표준오차. 이만큼은 벗어나는 것이 정상이다.
    #   평균:  sd(S)/sqrt(n)
    #   분산:  sd(S)^2 * sqrt(2/n)   (정규표본의 표본분산이 갖는 퍼짐)
    var_true = sigma1**2 + sigma2**2
    sd_true = np.sqrt(var_true)
    se_mean = sd_true / np.sqrt(n_simulations)
    se_var = var_true * np.sqrt(2 / n_simulations)
    print(f"  평균의 SE = {se_mean:.6f},  (관측 - 이론)/SE = {(S.mean() - 7) / se_mean:+.3f}")
    print(f"  분산의 SE = {se_var:.6f},  (관측 - 이론)/SE = {(S.var() - var_true) / se_var:+.3f}")

    # 표준편차는 더해지지 않는다. 더해지려면 상관계수가 1 이어야 한다.
    print(f"sd(S) = {S.std():.4f} (theoretical: {sd_true})"
          f"   sigma1 + sigma2 = {sigma1 + sigma2}")
    print(f"  rho=1 이면 Var = {sigma1**2 + sigma2**2 + 2 * sigma1 * sigma2},"
          f"  sd = {sigma1 + sigma2}")

verify_sum_of_normals()

출력:

E[X1+X2] = 7.0068 (theoretical: 7)
Var(X1+X2) = 25.1339 (theoretical: 25)
  평균의 SE = 0.015811,  (관측 - 이론)/SE = +0.432
  분산의 SE = 0.111803,  (관측 - 이론)/SE = +1.198
sd(S) = 5.0134 (theoretical: 5.0)   sigma1 + sigma2 = 7
  rho=1 이면 Var = 49,  sd = 7

모의값 \(7.0068\) 은 이론값 \(7\) 에서 \(0.43\) 표준오차, \(25.1339\) 는 \(25\) 에서 \(1.20\) 표준오차 떨어져 있다. 둘 다 몬테카를로 오차 안이며 어긋남이 아니다. 표본표준편차 \(5.0134\) 도 이론값 \(5\) 와 맞고, \(\sigma_1 + \sigma_2 = 7\) 과는 확연히 다르다.

되풀이를 \(10\) 만 번으로 늘려도 분산 쪽이 셋째 자리에서 흔들리는 것은 어쩔 수 없다. 분산의 표준오차가 \(n\) 의 제곱근에 반비례하므로, 소수 셋째 자리까지 맞추려면 되풀이를 \(10^4\) 배로 늘려야 한다. 모의실험은 유도한 식을 확인해 주지만 그 식을 대신하지는 못한다.

(a) \(t < \lambda\)에 대해 \(M_X(t) = \int_0^\infty e^{tx} \lambda e^{-\lambda x} dx = \lambda \int_0^\infty e^{-(\lambda - t) x} dx = \lambda/(\lambda - t)\)이다.

(b) \(M_X'(t) = \lambda/(\lambda - t)^2\)이므로 \(M_X'(0) = 1/\lambda = \mathbb{E}[X]\)이다. \(M_X''(t) = 2\lambda/(\lambda - t)^3\)이므로 \(M_X''(0) = 2/\lambda^2 = \mathbb{E}[X^2]\)이다.

(c) \(\mathrm{Var}(X) = 2/\lambda^2 - (1/\lambda)^2 = 1/\lambda^2\)이다. 평균과 표준편차가 모두 \(1/\lambda\)로 같다.

선형변환:

독립인 변수의 합: \(X\)와 \(Y\)가 독립이면 \(e^{tX}\)와 \(e^{tY}\)도 독립이므로(독립인 변수의 함수이므로)

이다.

두 번째 공식은 일반화된다. 상호독립인 \(S_n = X_1 + \cdots + X_n\)에 대해 \(M_{S_n}(t) = \prod_i M_{X_i}(t)\)이다.

이 두 성질 덕분에 적률생성함수가 확률변수의 변환과 합을 분석하는 데 그토록 유용하다.

포아송의 적률생성함수는 \(M_X(t) = e^{\lambda (e^t - 1)}\)이다.

합: \(M_{\sum X_i}(t) = \prod_i M_{X_i}(t) = \prod_i e^{\lambda_i (e^t - 1)} = e^{(\sum \lambda_i)(e^t - 1)}\).

이것이 \(\mathrm{Poisson}(\sum \lambda_i)\)의 적률생성함수다. 유일성에 의해 \(\sum X_i \sim \mathrm{Poisson}(\sum \lambda_i)\)이다.

함의: 포아송은 덧셈에 대해 "닫혀" 있다. 독립인 포아송들의 합은 비율을 합한 포아송이다. 이것이 계수 자료 모형의 토대다. 독립인 여러 원천에서 나온 총 사건 수가 결합된 비율의 포아송을 따르므로 모듈식 분석이 가능하다.

누율생성함수는 \(K_X(t) = \ln M_X(t)\)이며 그 테일러 전개는

이고 \(\kappa_n\)이 \(n\)차 누율이다. 처음 몇 개는 다음과 같다.

• \(\kappa_1 = \mathbb{E}[X] = \mu\) (평균)
• \(\kappa_2 = \mathrm{Var}(X) = \sigma^2\)
• \(\kappa_3 = \mathbb{E}[(X - \mu)^3]\) (3차 중심적률)
• \(\kappa_4 = \mathbb{E}[(X - \mu)^4] - 3\sigma^4\) (4차 중심적률에서 \(3\sigma^4\)을 뺀 값)

그러면 왜도 \(= \kappa_3/\kappa_2^{3/2}\)이고 초과첨도 \(= \kappa_4/\kappa_2^2\)이다. 둘 다 누율에서 나오며, 누율은 \(M_X(t)\)가 아니라 \(K_X(t)\)의 도함수다.

누율이 적률보다 나은 이유: 독립인 합의 누율은 더해진다. 독립인 \(X, Y\)에 대해 \(\kappa_n(X + Y) = \kappa_n(X) + \kappa_n(Y)\)이지만 적률은 그렇지 않다. 예를 들어 \(\mathrm{Var}(X + Y) = \mathrm{Var}(X) + \mathrm{Var}(Y)\)가 \(\kappa_2\)의 덧셈이고, 더 높은 누율에서도 마찬가지다. 그래서 누율이 중심극한정리와 에지워스 전개의 자연스러운 언어가 된다.

코시 밀도는 \(f(x) = 1/(\pi(1 + x^2))\)이다.

\(M(t) = \int_{-\infty}^\infty e^{tx} / (\pi (1 + x^2)) dx\)인데, \(t \ne 0\)이면 \(e^{tx}\)가 한쪽 방향으로 지수적으로 커져 피적분함수가 적분 가능하지 않다. 적률생성함수는 \(t = 0\)(자명하게 1)을 제외하면 정의되지 않는다.

특성함수(CF): \(\phi_X(t) = \mathbb{E}[e^{itX}]\). \(|e^{itX}| = 1\)이므로 \(\mathbb{E}[|e^{itX}|] = 1 < \infty\)가 되어 언제나 존재한다.

코시분포의 경우 적률생성함수가 없는데도 \(\phi(t) = e^{-|t|}\)라는 깔끔한 닫힌 형태를 갖는다.

실무적 함의: 꼬리가 두꺼운 분포를 다룰 때는 적률생성함수에서 특성함수로 갈아탄다. 대부분의 이론적 결과(레비 연속성 정리, 역변환 공식)는 더 넓은 부류를 포괄하는 특성함수로 진술된다. 적률생성함수는 모든 적률이 존재하는 분포, 즉 응용통계의 "얌전한" 다수에는 여전히 유용하다.

\(x > 0\)에서 확률밀도함수가 \(f_1(x) = \frac{1}{x\sqrt{2\pi}} e^{-(\ln x)^2/2}\)인 로그정규분포와 그것을 변형한

을 생각하자. 둘 다 타당한 확률밀도함수다(변형의 진폭을 음이 아니도록 잡았다).

적률: \(\mathbb{E}_2[X^k] = \int_0^\infty x^k f_1(x)(1 + \sin(2\pi \ln x)) dx = \mathbb{E}_1[X^k] + \int x^k f_1(x) \sin(2\pi \ln x) dx\).

\(u = \ln x\)로 치환하면 두 번째 적분이 \(\int e^{(k+1)u - u^2/2}\sin(2\pi u)/\sqrt{2\pi}\, du\)가 된다. 이 적분은 모든 정수 \(k\)에 대해 0이다. 피적분함수가 가우시안 인자와 적분값이 0이 되는 사인 변조를 결합하기 때문이다. 따라서 모든 적률이 일치하지만 \(f_1 \ne f_2\)이다.

왜 이렇게 되는가: 로그정규분포의 적률생성함수는 0의 어떤 근방에서도 수렴하지 않는다(\(t > 0\)이면 적분 \(\int e^{tx} f_1(x) dx\)가 발산한다). 그래서 적률생성함수가 존재하지 않고 적률의 열이 분포를 결정하지 못한다.

교훈: 적률생성함수가 열린구간에서 존재하면 적률이 분포를 유일하게 결정한다. 적률생성함수가 존재하지 않으면(두꺼운 꼬리, 로그정규) 적률만으로는 부족할 수 있으며 여러 분포가 모든 적률을 공유할 수 있다. 이것이 스틸체스 적률 문제다. 카를레만형 조건이 성립할 때에 한해 적률이 분포를 유일하게 결정한다.

임의의 \(t>0\)에 대해 \(\{X \ge a\} = \{e^{tX} \ge e^{ta}\}\)이므로 마르코프 부등식을 \(e^{tX}\)에 적용하면

이다. 이것이 모든 \(t>0\)에서 성립하므로 가장 좋은 \(t\)를 골라 최적화한다.

지수 안의 것이 \(\ln M_X\)의 르장드르 변환이며, 대편차이론에서 속도함수라 부른다.

import numpy as np
from scipy import optimize, stats

n, p = 100, 0.5
mu, var = n * p, n * p * (1 - p)
M = lambda t: (1 - p + p * np.exp(t)) ** n          # 이항분포의 적률생성함수

print(f"{'a':>5}{'정확한 값':>13}{'마르코프':>12}{'체비쇼프':>12}{'체르노프':>12}")
for a in (60, 70, 80):
    exact = stats.binom.sf(a - 1, n, p)
    markov = min(mu / a, 1.0)
    cheb = min(var / (a - mu) ** 2, 1.0)
    r = optimize.minimize_scalar(lambda t: -t * a + np.log(M(t)),
                                 bounds=(1e-9, 10), method="bounded")
    print(f"{a:>5}{exact:>13.3e}{markov:>12.3e}{cheb:>12.3e}{np.exp(r.fun):>12.3e}")

출력:

    a        정확한 값        마르코프        체비쇼프        체르노프
   60    2.844e-02   8.333e-01   2.500e-01   1.335e-01
   70    3.925e-05   7.143e-01   6.250e-02   2.670e-04
   80    5.580e-10   6.250e-01   2.778e-02   4.258e-09

격차가 극적이다. \(a=80\)에서 정확한 값이 \(5.6\times10^{-10}\)인데

마르코프와 체비쇼프는 \(a\)가 커져도 거의 줄지 않는다(각각 \(1/a\), \(1/a^2\)). 체르노프만이 지수적으로 줄어든다.

왜 이렇게 강한가. 마르코프는 \(X\) 자체를, 체비쇼프는 \(X^2\)을 쓴다. 체르노프는 \(e^{tX}\)를 쓰므로 모든 차수의 적률을 한꺼번에 동원한다. 그래서 적률생성함수가 존재하는 분포에서만 쓸 수 있고, 코시나 로그정규에서는 쓸 수 없다(연습문제 \(5\), \(6\)).

닫힌 형태. 이항분포에서는 최적화를 손으로 풀 수 있다. \(q=a/n\)일 때 최적 \(t^\ast = \ln\frac{q(1-p)}{p(1-q)}\)이고 한계가

가 된다. \(a=70\)이면 \(\mathrm{KL}(0.7\|0.5)=0.08228\)이라 한계가 \(e^{-8.228}=2.670\times10^{-4}\)로, 위 수치최적화 결과와 정확히 일치한다.

지수의 계수가 곧 정보량이다. 꼬리확률의 감소 속도가 쿨백–라이블러 발산으로 표현된다는 것은 우연이 아니며, 나중에 가설검정의 검정력과 대편차이론에서 다시 만난다. \(\square\)

음이 아닌 정수값 \(X\)에 대해 \(G_X(s)=\sum_{k\ge0}P(X=k)s^k\)이다. 적률생성함수와는 \(G_X(s)=M_X(\ln s)\)로 연결된다.

마지막 합성 성질이 확률생성함수만의 강점이다. 기댓값과 선형성 연습문제 \(9\)의 왈드 항등식이 여기서는 분포 전체로 강화된다.

분지과정. 한 개체가 자손을 \(X\)명 낳고(\(\mathbb{E}[X]=\lambda\)), 각 자손이 독립적으로 같은 일을 반복한다. \(n\)세대 개체 수의 확률생성함수는 \(G\)의 \(n\)중 합성 \(G^{\circ n}\)이다. 소멸확률 \(q=\lim_n P(Z_n=0)\)은

의 \([0,1]\)에서의 가장 작은 근이다. 그리고 \(\lambda \le 1\)이면 \(q=1\), \(\lambda>1\)이면 \(q<1\)이다.

import numpy as np
from scipy import optimize, stats

rng = np.random.default_rng(0)
print(f"{'lam':>6}{'소멸확률 (s=G(s))':>20}{'모의 (30세대)':>18}")
for lam in (0.8, 1.0, 1.5, 2.0):
    G = lambda s: np.exp(lam * (s - 1))           # 자손수 ~ Poisson(lam)
    q = optimize.brentq(lambda s: G(s) - s, 0, 1 - 1e-12) if lam > 1 else 1.0
    extinct = 0
    # 개체 수가 200 을 넘으면 더 따라가지 않고 "살아남았다" 로 친다.
    # 개체 m 개에서 출발한 가지가 모두 끊길 확률은 q^m 이므로, lam = 1.5 면
    # q^200 = 1e-76 이고 lam = 2 면 1e-139 다. 끊는 데서 생기는 치우침은
    # 20,000 회의 몬테카를로 오차 0.0035 보다 70 자릿수 작다.
    # 상한을 100,000 으로 두면 같은 답을 얻는 데 140 배 오래 걸린다.
    for _ in range(20_000):
        pop = 1
        for _ in range(30):
            if pop == 0 or pop > 200:
                break
            pop = rng.poisson(lam, pop).sum()
        extinct += (pop == 0)
    print(f"{lam:>6.1f}{q:>20.4f}{extinct / 20_000:>18.4f}")

print("\n합성 성질 확인: N ~ Poisson(3), X ~ Bernoulli(0.4) 이면 S ~ Poisson(1.2)")
N = rng.poisson(3.0, 200_000)
S = rng.binomial(N, 0.4)
for k in range(4):
    print(f"  P(S={k}) 모의 {np.mean(S == k):.4f}   Poisson(1.2) "
          f"{stats.poisson.pmf(k, 1.2):.4f}")

출력:

   lam       소멸확률 (s=G(s))         모의 (30세대)
   0.8              1.0000            0.9994
   1.0              1.0000            0.9377
   1.5              0.4172            0.4173
   2.0              0.2032            0.2033

합성 성질 확인: N ~ Poisson(3), X ~ Bernoulli(0.4) 이면 S ~ Poisson(1.2)
  P(S=0) 모의 0.3004   Poisson(1.2) 0.3012
  P(S=1) 모의 0.3619   Poisson(1.2) 0.3614
  P(S=2) 모의 0.2169   Poisson(1.2) 0.2169
  P(S=3) 모의 0.0866   Poisson(1.2) 0.0867

임계점 \(\lambda=1\)이 특이하다. 이론적으로 소멸확률이 \(1\)인데 \(30\)세대 모의에서는 \(0.9377\)밖에 되지 않는다. 임계 분지과정은 확률 \(1\)로 소멸하지만 매우 느리게 소멸한다 — \(P(Z_n>0)\sim 2/(n\sigma^2)\)로 \(1/n\) 속도라, 모든 유한 세대에서는 상당수가 살아남아 있다.

합성 성질도 확인된다. \(G_S(s)=G_N(G_X(s))=e^{3(0.4s+0.6-1)}=e^{1.2(s-1)}\)이므로 \(S\sim\text{Poisson}(1.2)\)이고, 모의 분포가 정확히 그와 일치한다. "포아송 개수의 베르누이 성공"이 다시 포아송이라는 이 사실이 포아송 씨닝(thinning)이라 불리며 대기행렬과 점과정에서 기본 도구다.

\(\lambda=1\)이 임계점인 이유. \(G\)는 볼록증가함수이고 \(G(1)=1\)이다. \(G'(1)=\lambda\)가 \(1\)보다 크면 \(s=1\) 직전에서 \(G\)가 대각선 아래로 내려가므로 \([0,1)\) 안에 교점이 하나 더 생긴다. 평균 자손 수가 \(1\)을 넘느냐가 전부다 — 감염병의 재생산지수 \(R_0\)가 정확히 이 \(\lambda\)다. \(\square\)

\(K_X(t)=\ln M_X(t)\)가 누율생성함수이고 \(\kappa_n = K^{(n)}(0)\)이 \(n\)번째 누율이다. 독립이면 \(M_{X+Y}=M_XM_Y\)이므로 로그를 취해

적률은 이런 성질이 없다(\(\mathbb{E}[(X+Y)^2]\ne\mathbb{E}[X^2]+\mathbb{E}[Y^2]\)). 누율이 적률보다 자연스러운 좌표계인 이유다.

import numpy as np

def cumulants(x):                       # 처음 네 개의 표본누율
    c = x - x.mean()
    m2, m3, m4 = (c ** 2).mean(), (c ** 3).mean(), (c ** 4).mean()
    return x.mean(), m2, m3, m4 - 3 * m2 ** 2

rng = np.random.default_rng(0)
n = 4_000_000
X = rng.poisson(3.0, n)                 # 누율이 모두 3
Y = rng.gamma(2.0, 1.0, n)              # 누율이 (n-1)! * 2

for name, v, theory in [("X ~ Poisson(3)", X, [3, 3, 3, 3]),
                        ("Y ~ Gamma(2,1)", Y, [2, 2, 4, 12]),
                        ("X + Y", X + Y, [5, 5, 7, 15])]:
    k = cumulants(v)
    print(f"{name:>14}: k1 {k[0]:>7.3f}  k2 {k[1]:>7.3f}  k3 {k[2]:>7.3f}"
          f"  k4 {k[3]:>7.3f}   이론 {theory}")

출력:

X ~ Poisson(3): k1   2.999  k2   3.001  k3   3.009  k4   2.984   이론 [3, 3, 3, 3]
Y ~ Gamma(2,1): k1   2.000  k2   1.998  k3   3.995  k4  11.995   이론 [2, 2, 4, 12]
         X + Y: k1   4.999  k2   4.999  k3   7.009  k4  15.020   이론 [5, 5, 7, 15]

네 개 모두 자리별로 더해진다. \(3+2=5\), \(3+2=5\), \(3+4=7\), \(3+12=15\).

정규분포의 특징화. \(N(\mu,\sigma^2)\)의 누율생성함수는

로 \(t\)의 이차식이다. 따라서 \(\kappa_1=\mu\), \(\kappa_2=\sigma^2\), 그리고 \(\kappa_n=0\) for all \(n\ge3\)이다. 왜도와 첨도가 0인 것(연습문제 \(4\))이 이 사실의 처음 두 항일 뿐이다.

이것이 정규분포의 지위를 설명한다.

• 누율이 더해지므로, 독립합을 표준화하면 \(\kappa_n\)이 \(n\ge3\)에서 \(m^{1-n/2}\) 비율로 사라진다. 남는 것은 \(\kappa_1,\kappa_2\)뿐이고 그것이 곧 정규분포다 — 중심극한정리의 누율 버전이다.
• 고차 누율을 조금씩 되살려 보정한 것이 에지워스 전개이며, 3.5절 베리–에센 한계와 함께 "얼마나 정규에 가까운가"를 정량화한다.

실무적 함의. 자료가 비정규임을 진단할 때 왜도(\(\kappa_3/\sigma^3\))와 초과첨도(\(\kappa_4/\sigma^4\))를 보는 것은 정규분포에서 \(0\)이어야 하는 양을 재는 것이다. 그리고 표본 크기를 \(m\)배 늘려 평균을 내면 이 두 값이 각각 \(\sqrt{m}\), \(m\)배로 줄어든다. \(\square\)

import numpy as np

rng = np.random.default_rng(1)
print("X ~ N(0,1),  n = 10000,  참값 M(t) = exp(t^2/2)")
print(f"{'t':>5}{'참값':>10}{'추정 (같은 조건 5회 반복)':>46}")
for t in (0.5, 1.0, 2.0, 3.0):
    est = [np.mean(np.exp(t * rng.normal(0, 1, 10_000))) for _ in range(5)]
    print(f"{t:>5.1f}{np.exp(t ** 2 / 2):>10.3f}   "
          + " ".join(f"{e:>8.3f}" for e in est))

print("\nX ~ Exp(1),  M(t) = 1/(1-t) 는 t < 1 에서만 유한하다")
for t in (0.5, 0.9, 1.5):
    x = rng.exponential(1.0, 100_000)
    true = 1 / (1 - t) if t < 1 else float("inf")
    print(f"  t={t}:  추정 {np.mean(np.exp(t * x)):>14.3f}   참값 {true:>10.3f}")

출력:

X ~ N(0,1),  n = 10000,  참값 M(t) = exp(t^2/2)
    t        참값                              추정 (같은 조건 5회 반복)
  0.5     1.133      1.127    1.124    1.128    1.128    1.131
  1.0     1.649      1.620    1.661    1.669    1.631    1.642
  2.0     7.389      8.411    6.678    7.169    7.353    7.376
  3.0    90.017     59.410   86.072  106.658   67.707   85.592

X ~ Exp(1),  M(t) = 1/(1-t) 는 t < 1 에서만 유한하다
  t=0.5:  추정          1.981   참값      2.000
  t=0.9:  추정          6.746   참값     10.000
  t=1.5:  추정   14997922.650   참값        inf

\(t\)가 커질수록 추정이 무너진다. \(t=0.5\)에서는 \(1.124\)–\(1.131\)로 안정적인데, \(t=3\)에서는 참값 \(90.0\)에 대해 \(59.4\)부터 \(106.7\)까지 흩어진다.

원인은 추정량 자체의 분산이다. \(\operatorname{Var}(\hat M(t)) = \frac{1}{n}\left(M(2t)-M(t)^2\right)\)인데, 정규분포에서 이 값이

로 \(t^2\)의 지수로 폭발한다. \(t=3\)이면 \(e^{18}\approx 6.6\times10^7\)이라 \(n=10^4\)로는 어림도 없다. 그리고 \(\hat M(t)\)는 실질적으로 표본 중 가장 큰 몇 개가 결정한다 — 지수가 나머지를 전부 눌러버리기 때문이다.

더 나쁜 경우: 존재하지 않는 값을 추정한다. \(\mathrm{Exp}(1)\)에서 \(M(t)\)는 \(t\ge1\)이면 무한대다. 그런데 \(t=1.5\)를 넣어도 프로그램은 조용히 \(1.5\times10^{7}\)이라는 유한한 수를 내놓는다. 표본에는 유한 개의 값뿐이므로 \(\hat M(t)\)는 언제나 유한하다. 발산을 자료만 보고 알아낼 방법이 없다. \(t=0.9\)에서 이미 추정값 \(6.75\)가 참값 \(10\)에 크게 못 미치는 것이 그 전조다.

그래서 실무에서는 무엇을 쓰는가.

특성함수는 사정이 다르다. \(|e^{itX}|=1\)이라 항상 유계이고, 경험적 특성함수 \(\frac1n\sum e^{itX_j}\)는 모든 \(t\)에서 분산이 \(1/n\) 이하로 잘 행동한다. 실제로 특성함수 기반의 적합도 검정이 존재하는 이유이며, 적률생성함수 기반의 것은 없는 이유다. \(\square\)