적분과 변수변환

다중적분

2차원 이상의 확률분포에서 확률·기댓값·주변분포를 구하려면 여러 변수에 대한 적분이 필요하다. 다중적분을 계산하는 실제 방법은 거의 항상 한 변수씩 차례로 적분하는 것이고, 이때 영역의 경계를 정확히 기술하는 것이 핵심이다.

정의. 영역 $\mathcal{R} \subset \mathbb{R}^2$ 위의 함수 $f$의 적분 $\iint_{\mathcal{R}} f\,dA$는 영역을 작은 칸으로 나누고 (칸 넓이) × (칸 안의 함수값)을 더한 합(Riemann 합)의 극한이다. 1차원 적분이 곡선 아래 넓이라면, 2차원 적분은 곡면 아래 부피다. $f$가 확률밀도이면 그 영역에 떨어질 확률이다.

반복적분과 Fubini 정리. $f$가 영역 위에서 연속이거나, 더 일반적으로 $\iint \lvert f \rvert\,dA < \infty$이면 다중적분은 반복적분(iterated integral)과 같고, 적분 순서를 바꿔도 값이 같다(Fubini 정리). 직사각형 영역에서는 한계가 상수이므로 순서 교환이 단순하다.

영역이 직사각형이 아닐 때. 안쪽 적분의 한계가 바깥 변수에 의존한다. 삼각형 $0 \le y \le x \le 1$에서 $x$를 바깥에 두면 $y$는 $0$부터 $x$까지, $y$를 바깥에 두면 $x$는 $y$부터 $1$까지 움직인다.

순서를 바꿀 때는 한계를 기계적으로 옮기지 말고, 영역을 그린 뒤 바깥 변수를 고정한 선분이 영역에 들어가고 나가는 위치를 읽는다. $f(x, y) = xy^2$이면 앞의 순서로 $\int_0^1 x \cdot x^3/3\,dx = 1/15$, 뒤의 순서로 $\int_0^1 y^2(1 - y^2)/2\,dy = 1/15$다.

import numpy as np
from scipy import integrate

f = lambda x, y: x * y**2
# dblquad(func(안쪽, 바깥), 바깥 하한, 바깥 상한, 안쪽 하한(바깥), 안쪽 상한(바깥))
I_dy_dx, _ = integrate.dblquad(lambda y, x: f(x, y), 0, 1, lambda x: 0, lambda x: x)
I_dx_dy, _ = integrate.dblquad(lambda x, y: f(x, y), 0, 1, lambda y: y, lambda y: 1)

n = 2000                                            # 격자 Riemann 합
c = (np.arange(n) + 0.5) / n                        # (n,) 칸의 중심
Xg, Yg = np.meshgrid(c, c, indexing="ij")           # (n, n)
riemann = np.sum(f(Xg, Yg) * (Yg <= Xg)) / n**2
print(f"dy 먼저 {I_dy_dx:.6f}   dx 먼저 {I_dx_dy:.6f}   Riemann {riemann:.6f}   1/15 = {1/15:.6f}")

dy 먼저 0.066667   dx 먼저 0.066667   Riemann 0.066729   1/15 = 0.066667

두 순서의 반복적분이 같고, 격자 Riemann 합도 가깝다. Riemann 합이 네 번째 자리에서 어긋나는 것은 대각선 경계에 걸친 칸을 통째로 넣거나 빼기 때문이다. 영역의 경계를 칸으로 근사하는 오차는 칸 크기에 비례해 줄어 내부 오차보다 느리게 줄어든다. 수치적분 방법은 9.5장에서 다룬다.

변수변환

적분 영역이 원이나 회전된 타원이면 직교좌표의 한계가 복잡해진다. 영역이나 피적분함수의 대칭에 맞는 좌표로 바꾸면 적분이 단순해지는데, 그 대가로 부피 배율을 곱해야 한다.

선형 변환의 경우. 선형사상 $\mathbf{x} \mapsto \mathbf{A}\mathbf{x}$는 단위 정사각형을 두 열벡터가 만드는 평행사변형으로 보내고, 그 넓이는 $\lvert\det\mathbf{A}\rvert$다. 모든 영역은 작은 정사각형으로 채울 수 있으므로, 선형사상은 모든 영역의 넓이를 똑같이 $\lvert\det\mathbf{A}\rvert$배 한다. 행렬식의 부호는 방향(시계·반시계)을 뒤집는지를 나타낼 뿐 넓이와는 관계없다.

비선형 변환의 경우. 8.1장에서 본 것처럼 미분 가능한 사상 $g$는 한 점 근처에서 아핀 사상 $g(\mathbf{u} + \mathbf{h}) \approx g(\mathbf{u}) + \mathbf{J}_g(\mathbf{u})\mathbf{h}$처럼 행동한다. 그래서 점 $\mathbf{u}$ 근처의 작은 칸은 넓이가 $\lvert\det\mathbf{J}_g(\mathbf{u})\rvert$배가 된다. 배율이 점마다 다를 뿐 원리는 같다.

변수변환 공식. $g: \mathcal{U} \to \mathcal{X}$가 가역이고 연속적으로 미분 가능하며 역사상도 그렇다면

이다. 왼쪽은 $\mathbf{x}$ 공간의 칸으로 더한 합, 오른쪽은 $\mathbf{u}$ 공간의 칸으로 더한 합이다. $\mathbf{u}$ 공간의 칸 $d\mathbf{u}$는 $\mathbf{x}$ 공간에서 넓이 $\lvert\det\mathbf{J}_g\rvert d\mathbf{u}$인 칸에 대응하므로 그 배율을 곱한다. 1차원에서는 익숙한 치환적분 $dx = \lvert g'(u) \rvert du$다.

극좌표. $(r, \theta) \mapsto (r\cos\theta, r\sin\theta)$의 Jacobian 행렬식은 8.1장에서 계산했듯 $r$이다. 그래서

반지름 $r$에서 각도 $d\theta$가 만드는 호의 길이가 $r\,d\theta$라는 기하와 맞는다. 원점에서 멀수록 같은 $(dr, d\theta)$ 칸이 넓은 영역에 대응한다.

두 사실을 계산으로 확인한다. 첫째, 선형사상 $\mathbf{A}$가 단위 정사각형을 보낸 평행사변형의 넓이를 Monte Carlo로 재서 $\lvert\det\mathbf{A}\rvert$와 비교한다. 둘째, 극좌표 격자로 단위원의 넓이를 더할 때 $r$을 곱한 경우와 빠뜨린 경우를 비교한다.

rng = np.random.default_rng(0)
A = np.array([[2.0, 1.0], [0.5, 1.5]])
corners = np.array([[0, 0], [1, 0], [0, 1], [1, 1]]) @ A.T   # (4, 2) 정사각형 꼭짓점의 상
lo, hi = corners.min(0), corners.max(0)
box = rng.uniform(lo, hi, size=(200_000, 2))        # (N, 2) 경계 상자 위 균등 표본
pre = np.linalg.solve(A, box.T).T                    # (N, 2) 역사상으로 되돌린 점
frac = np.mean(np.all((pre >= 0) & (pre <= 1), axis=1))   # 평행사변형 안의 비율
print(f"평행사변형 넓이 (MC) {frac * np.prod(hi - lo):.3f}   |det A| = {abs(np.linalg.det(A)):.3f}")

m = 400                                              # 극좌표 격자로 단위원의 넓이
r = (np.arange(m) + 0.5) / m                         # (m,)
dr, dth = 1 / m, 2 * np.pi / m
print(f"원의 넓이: r 포함 {np.sum(r) * m * dr * dth:.4f}   r 빠뜨림 {m * m * dr * dth:.4f}   pi = {np.pi:.4f}")

평행사변형 넓이 (MC) 2.500   |det A| = 2.500
원의 넓이: r 포함 3.1416   r 빠뜨림 6.2832   pi = 3.1416

$r$을 빠뜨리면 $(r, \theta)$ 직사각형 $[0, 1] \times [0, 2\pi]$의 넓이 $2\pi$를 계산하게 된다. 좌표를 바꿨다는 사실만으로 넓이가 두 배로 틀린다.

가우시안 적분

정규분포의 밀도 앞에 붙는 $1/\sqrt{2\pi\sigma^2}$은 어디서 오는가. 정규화 상수를 직접 유도해 보면 변수변환이 왜 강력한지 보이고, 다변량 가우시안의 $\det(\boldsymbol{\Sigma})^{-1/2}$이 무엇을 뜻하는지도 이해된다.

1차원. $I = \int_{-\infty}^{\infty} e^{-x^2/2}\,dx$는 초등함수로 된 부정적분이 없다. 대신 $I^2$을 2차원 적분으로 쓰고 극좌표로 바꾼다.

따라서 $I = \sqrt{2\pi}$다. 여기서 Jacobian $r$이 결정적이다. $r e^{-r^2/2}$는 $-e^{-r^2/2}$의 도함수라서 적분이 닫힌 형태로 풀린다. 이어서 $x = \mu + \sigma u$로 치환하면 $dx = \sigma\,du$이므로 $\int e^{-(x - \mu)^2/(2\sigma^2)}dx = \sigma\sqrt{2\pi}$이고, 정규분포의 밀도는

가 된다.

$d$차원. 공분산 $\boldsymbol{\Sigma}$(양의 정부호)를 Cholesky 분해 $\boldsymbol{\Sigma} = \mathbf{L}\mathbf{L}^\top$하고 $\mathbf{x} = \mathbf{L}\mathbf{u}$로 바꾸면 $\mathbf{x}^\top\boldsymbol{\Sigma}^{-1}\mathbf{x} = \lVert\mathbf{u}\rVert^2$이고 $d\mathbf{x} = \lvert\det\mathbf{L}\rvert\,d\mathbf{u} = \det(\boldsymbol{\Sigma})^{1/2}d\mathbf{u}$다. 남는 적분은 1차원 가우시안 적분 $d$개의 곱이다.

그래서 다변량 정규분포의 밀도 앞에는 $(2\pi)^{-d/2}\det(\boldsymbol{\Sigma})^{-1/2}$이 붙는다. $\det\boldsymbol{\Sigma}$는 분포가 차지하는 부피의 제곱이고, 넓게 퍼진 분포일수록 밀도의 높이가 낮아져야 전체 확률이 1로 유지된다. 다변량 가우시안의 성질은 11.4장에서 자세히 다룬다.

I1, _ = integrate.quad(lambda x: np.exp(-x**2 / 2), -np.inf, np.inf)
print(f"1차원: {I1:.6f}   sqrt(2 pi) = {np.sqrt(2 * np.pi):.6f}")

Sigma = np.array([[2.0, 0.8], [0.8, 1.0]])           # (2, 2) 공분산
P = np.linalg.inv(Sigma)
g = np.linspace(-10, 10, 801)
dx = g[1] - g[0]
G = np.stack(np.meshgrid(g, g, indexing="ij"), axis=-1)   # (801, 801, 2)
q = np.einsum("...i,ij,...j->...", G, P, G)               # (801, 801) 이차형식
I2 = np.exp(-0.5 * q).sum() * dx**2
print(f"2차원: {I2:.6f}   2 pi sqrt(det) = {2 * np.pi * np.sqrt(np.linalg.det(Sigma)):.6f}")

1차원: 2.506628   sqrt(2 pi) = 2.506628
2차원: 7.327390   2 pi sqrt(det) = 7.327390

정규화 상수는 학습에서 사라지지 않는다. 분산을 학습하는 모델(12.3장의 이분산 회귀, 13장의 베이지안 모델)에서 음의 로그우도에 $\log\sigma$ 또는 $\tfrac{1}{2}\log\det\boldsymbol{\Sigma}$ 항이 남는 것이 이 정규화 상수 때문이다. 이 항이 없으면 모델은 분산을 무한히 키워 오차 항을 0으로 만들 수 있다.

확률밀도의 변환

측정값의 단위를 바꾸거나, 로그를 취하거나, 신경망으로 잠재변수를 데이터로 보낼 때마다 확률변수가 변환된다. 확률밀도는 단위 부피당 확률이므로, 부피가 늘어나는 곳에서는 밀도가 그만큼 낮아져야 한다.

1차원 유도. $X$의 밀도가 $p_X$이고 $Y = g(X)$, $g$는 증가하는 가역 함수라 하자. 누적분포함수로 쓰면

이다. 감소 함수이면 부호가 바뀌므로 일반적으로는 절댓값을 쓴다. 역함수의 도함수는 $1/g'(x)$이므로 $p_Y(y) = p_X(x)/\lvert g'(x)\rvert$($x = g^{-1}(y)$)로도 쓴다.

다변량. $\mathbf{Y} = g(\mathbf{X})$이고 $g$가 가역이며 미분 가능하면, 변수변환 공식에서 바로

를 얻는다. 임의의 영역 $\mathcal{B}$에 대해 $P(\mathbf{Y} \in \mathcal{B}) = P(\mathbf{X} \in g^{-1}(\mathcal{B}))$가 성립해야 하고, 오른쪽 적분을 $\mathbf{y}$로 변수변환하면 이 식이 나온다. 로그를 취하면

예: 로그정규분포. $X \sim \mathcal{N}(0, 1)$, $Y = e^X$이면 $g^{-1}(y) = \log y$, $\lvert d g^{-1}/dy\rvert = 1/y$이므로 $p_Y(y) = p_X(\log y)/y$다. 표본의 구간별 비율과 공식의 적분을 비교하고, Jacobian 항 $1/y$를 빠뜨린 "밀도"의 전체 적분도 계산한다.

x = rng.normal(size=1_000_000)                       # X ~ N(0, 1)
y = np.exp(x)                                        # Y = exp(X), 로그정규분포
p_x = lambda t: np.exp(-t**2 / 2) / np.sqrt(2 * np.pi)
p_y = lambda s: p_x(np.log(s)) / s                   # |d g^{-1}/dy| = 1/y
p_wrong = lambda s: p_x(np.log(s))                   # Jacobian 항을 빠뜨림

edges = np.array([0.25, 0.5, 1.0, 2.0, 4.0])
counts, _ = np.histogram(y, bins=edges)
for a, b, cnt in zip(edges[:-1], edges[1:], counts):
    exact, _ = integrate.quad(p_y, a, b)
    print(f"P({a:4.2f} < Y < {b:4.2f})  표본 {cnt / y.size:.4f}  공식 {exact:.4f}")
print(f"전체 적분: 올바른 밀도 {integrate.quad(p_y, 0, np.inf)[0]:.4f}"
      f"   Jacobian 누락 {integrate.quad(p_wrong, 0, np.inf)[0]:.4f}")

P(0.25 < Y < 0.50)  표본 0.1616  공식 0.1613
P(0.50 < Y < 1.00)  표본 0.2556  공식 0.2559
P(1.00 < Y < 2.00)  표본 0.2565  공식 0.2559
P(2.00 < Y < 4.00)  표본 0.1613  공식 0.1613
전체 적분: 올바른 밀도 1.0000   Jacobian 누락 1.6487

Jacobian을 빠뜨린 함수의 전체 적분은 $\int p_X(\log y)\,dy = \int p_X(x)\,e^x\,dx = e^{1/2} \approx 1.6487$이다. 확률밀도가 아니다. 이 함수로 최대우도 학습을 하면 손실이 "확률"이 아닌 양을 키우는 방향으로 움직인다.

normalizing flow. 이 공식은 생성모델의 한 갈래인 normalizing flow의 출발점이다. 다루기 쉬운 분포(보통 표준정규) $\mathbf{z} \sim p_Z$에서 시작해 학습 가능한 가역 변환 $\mathbf{x} = f_\theta(\mathbf{z})$로 데이터 분포를 표현한다. 데이터 $\mathbf{x}$의 로그우도는

이고, 이것을 최대화하는 것이 학습이다. $f_\theta = f_K \circ \cdots \circ f_1$처럼 여러 층을 합성하면 chain rule에 의해 Jacobian은 층별 Jacobian의 곱이고, 곱의 행렬식은 행렬식의 곱이므로 로그 행렬식은 층별로 더해진다.

설계의 핵심 제약은 계산량이다. $d \times d$ 행렬의 행렬식은 일반적으로 $O(d^3)$ 연산이 들어 고차원 데이터에서 매 step 계산할 수 없다. 그래서 flow의 층은 Jacobian이 삼각행렬이 되도록 설계하는 경우가 많다. 삼각행렬의 행렬식은 대각 성분의 곱이므로 $O(d)$로 계산된다. 실습에서 쓰는 변환이 그 가장 작은 예이고, 구체적인 구조는 19.4장에서 다룬다.