• Using Inverse & Rejection Methods

• Multivariate Normal Distribution

• Copulas (Joint Distributions)

• Gaussian Copulas

고차원 확률분포에서 샘플링을 하는 것은 쉬운 일이 아니며, 이는 대표적으로 차원의 저주라는 문제 때문에 하기가 어렵습니다.

확률분포 내의 변수들이 독립적이라면 사실 그저 모든 변수들을 독립적으로 뽑으면 되지만, 상관관계가 있다면 그것을 고려한채로 뽑는 것이 어렵습니다.

Using Inverse & Rejection Methods

Multivariate Inverse Transform Method

어떠한 다변량 분포에 대한 CDF 를 $F_{1, \cdots, d} (x_1, \cdots , x_d)$ 가 있다면 이것을 다음과 같이 쪼개는 것이 가능합니다.

그러면 저희는 다음과 같이 Inverse Transform Method 를 적용하는 것이 가능합니다.

물론 이렇게 진행하면 각각에 대한 식을 구해야 된다는 단점이 존재하며, $d$ 크기가 커지면 커질수록 이 식을 구하는 것은 Intractable 해집니다.

Example

다변량 정규분포 특성상 조건부 분포는 다음과 같이 구할 수가 있습니다.

이것은 Generalize 될 수가 있으며, 이것을 다음과 같은 예시 분포에 적용해보겠습니다.

1. $x_1$ 에 대해서 :
   • $X_1 \sim \mathrm{N}(0, 1.0)$
   • $F_1(x_1) = \Phi\left( \frac{x_1 - 0}{\sqrt{1.0}} \right)$
   • $F_1^{-1}(u) = 0 + \sqrt{1.0} \cdot \Phi^{-1}(u) = x_1$

2. $x_2$ 에 대해서 :
   • $X_2 \mid x_1 \sim \mathrm{N}\left( 0+\frac{0.4}{1.0}\left(X_1-0\right), 1.5-\frac{0.4^2}{1.0} \right) = \mathrm{N}\left( 0.4 \cdot X_1, 1.34 \right)$
   • $F_{2\mid 1}(x_2\mid x_1) = \Phi\left( \frac{x_2 - 0.4 \cdot x_1}{\sqrt{1.34}} \right)$
   • $F_{2\mid 1}^{-1}(u\mid x_1) = 0.4 \cdot x_1 + \sqrt{1.34} \cdot \Phi^{-1}(u)$

3. $x_3$ 에 대해서 :
   • $\Sigma_{11:22} = \left[\begin{matrix} 1.0 & 0.4 \\ 0.4 & 1.5 \end{matrix} \right], \quad \Sigma_{13:23} = \left[\begin{matrix} 0.2 \\ 0.3 \end{matrix}\right], \quad \Sigma_{31:32} = \left[\begin{matrix} 0.2 & 0.3 \end{matrix}\right], \quad \Sigma_{33} = 2.0$
   • $X_3 \mid x_1, x_2 \sim \mathrm{N}\left( 0 + \Sigma_{13:23} \Sigma^{-1}_{11:22}\left( \left[\begin{matrix} x_1 \\ x_2 \end{matrix} \right] \right) , \; \Sigma_{33} - \Sigma_{31:32} \Sigma_{11:22}^{-1} \Sigma_{13:23} \right)\\ =: \mathrm{N}\left( m_{3\mid 12}, v_{3 \mid 12}\right)$
   • $F_{3\mid 1, 2}(x_3\mid x_1, x_2) = \Phi\left( \frac{x_3 - m_{3 \mid 12}}{\sqrt{v_{3 \mid 12}}} \right)$
   • $F^{-1}_{3\mid 12} (u \mid x_1, x_2) = m_{3 \mid 12} + \sqrt{v_{3 \mid 12}} \cdot \Phi^{-1}\left( u\right)$

이것을 코드로 진행하게 된다면 다음과 같이 나오게 됩니다.

mu    <- c(0, 0, 0)
Sigma <- matrix(c(
  1.0, 0.4, 0.2,
  0.4, 1.5, 0.3,
  0.2, 0.3, 2.0
), nrow = 3, byrow = TRUE)

sample_mvnorm_inverse <- function(mu, Sigma) {
  # X1 ~ N(mu1, Sigma[1,1])
  u1 <- runif(1); z1 <- qnorm(u1)
  x1 <- mu[1] + sqrt(Sigma[1,1]) * z1
  
  # X2 | X1=x1 ~ N(m21, v21)
  m21 <- mu[2] + (Sigma[2,1]/Sigma[1,1]) * (x1 - mu[1])
  v21 <- Sigma[2,2] - Sigma[2,1]^2 / Sigma[1,1]
  u2  <- runif(1);   z2 <- qnorm(u2)
  x2  <- m21 + sqrt(v21) * z2
  
  # X3 | (X1,X2) ~ N(m312, v312)
  S11  <- Sigma[1:2, 1:2]
  S21  <- Sigma[3, 1:2]
  S12  <- Sigma[1:2, 3]
  invS <- solve(S11)
  delta <- c(x1, x2) - mu[1:2]
  m312  <- mu[3] + as.numeric(S21 %*% invS %*% delta)
  v312  <- Sigma[3,3] - as.numeric(S21 %*% invS %*% S12)
  u3    <- runif(1); z3 <- qnorm(u3)
  x3    <- m312 + sqrt(v312) * z3
  
  c(x1, x2, x3)
}

손으로 연산만 해도 굉장히 힘들며 오래 걸립니다. 식을 매번 구해야 된다는 단점이 있기 때문이죠. 추가적으로 Inverse CDF 도 매번 구해야 되는 단점도 있으며 이게 존재하지 않을 수도 있습니다.

반면에 Rejection Method 를 다변량에 적용시키면 어떨까요? 똑같은 방법으로 적용시키면 되지만, 좋은 Proposal Density $q$ 를 구하는 것이 어려우며, 대부분의 데이터가 기각당하는 경우도 빈번합니다. 그래서 이를 다르게 활용한 방법으로 MCMC 가 존재하며, 이는 나중에 다루도록 하겠습니다.

Multivariate Normal Distribution

다변량 정규 분포 $\mathbf{X} \sim \mathrm{N}_n\left(\boldsymbol{\mu}, \mathbf{\Sigma}\right)$ 는 다음과 같은 pdf 를 가집니다.

하지만 정규분포들끼리는 서로 선형적인 관계를 가진다는 것을 저희는 알고 있습니다. 다음과 같은 독립적인 분포의 다변량 표준정규분포는 쉽게 뽑을 수가 있습니다. 그저 표준정규분포에서 i.i.d 하게 뽑은 다음에 그것을 벡터 형태로 바꿔주면 되죠.

$\mathbf{X}$ 에 대해서는 그러면 선형적인 관계임을 이용해서 다음과 같이 구할 수가 있습니다.

이때의 가정은 $\mathbf{AA}^\top$ 가 Positive Definite 이라는 가정이 들어갑니다.

그러면 만약에 $\mathbf{X}$ 에서 뽑고 싶으면 그저 $\mathbf{\Sigma} = \mathbf{AA}^\top$ 인 Matrix $\mathbf{A}$ 를 찾은 다음에 $\mathbf{X} := \mathbf{AZ}+\boldsymbol \mu$ 으로 세팅해준다면 $\mathbf{X}$ 을 뽑을 수가 있게 됩니다.

이때 Cholesky Decomposition 을 활용해서 $\mathbf{A}$ 를 찾으면 더 효율적이게 됩니다.

Cholesky Decomposition 이란 $\mathbf{\Sigma}$ 와 같은 Symmetric + Positive Definite Matrix 를 Decompose 하기 위한 굉장히 효율적인 방법이며, 이렇게 하게 되면 $\mathbf{A}$ 가 Unique 한 하삼각행렬로 $\mathbf{\Sigma} = \mathbf{AA}^\top$ 으로 구성하는 것이 가능해집니다.

이것을 위의 예시 문제에 다시 한번 적용시켜보겠습니다.

n <- 10000
d <- 3

z_vals <- qnorm(runif(n * d))
z_mat <- matrix(data = z_vals, nrow = 3)

mu    <- c(0, 0, 0)
Sigma <- matrix(c(
  1.0, 0.4, 0.2,
  0.4, 1.5, 0.3,
  0.2, 0.3, 2.0
), nrow = 3, byrow = TRUE)

A <- t(chol(Sigma))

X <- A %*% z_mat + mu

이렇게 진행하게 된다면 위의 Inverse CDF 방법보다 더 빠르게 샘플링을 할 수 있게 됩니다.

Copulas (Joint Distributions)

Copula 란 그저 $[0,1]^d$ 위에 있는 Joint Distribution 을 의미합니다. 쉽게 말하자면, CDF 의 다변량 버전이며, correlation 까지 알려주는 cdf 함수라고 생각하면 됩니다.

예시로 2-차원일 때의 Copula 는 $C(x,y), \; 0 \le x \le 1, \; 0\le y \le 1$ 이며, 다음과 같은 특성들을 가집니다.

• $C(x,0) = 0, \quad C(0, y) = 0,$

• $C(x,1) = x, \quad C(1,y) = y$

이 예시를 이어가서 만약에 Joint Distribution 이 $(X,Y) \sim H$ 이라고 하며 Marginal Distribution 들이 $X\sim F$, $Y\sim G$ 이라고 한다면, $H$ 의 Copula 인 $C_H$ 는 다음과 같이 표현이 가능합니다.

이를 Sklar’s 정리에 의해서 Copula 도 알고 Marginal 도 안다면 다음과 같이 Joint Distribution 인 $H(x,y)$ 도 다음과 같이 구할 수가 있습니다.

이게 무슨 말이냐면, $(U_1, U_2)$ 를 $C_H$ 에서 뽑을 수가 있다면, 그러면 $U_1$ 과 $U_2$ 는 상관관계가 있는 확률 값들이며, 이들을 Inverse Transform Method 으로 $X = F^{-1}(U_1), \; Y = G^{-1}(U_2)$ 으로 세팅하게 된다면 Joint Distribution 에서 뽑을 수 있게 됩니다.

반대로 얘기하자면, $(X,Y) \sim H$ 에서 뽑을 수 있게 되면, $U_1 = F(X), \; U_2 = G(Y)$ 으로 세팅하게 된다면 $(U_1, U_2) \sim C_H$ 에서 뽑게 되는 셈이랑 똑같게 된다는 것입니다.

예시로 Gaussian Copula 를 보겠습니다.

Gaussian Copula

Copula 는 그저 $[0,1]^d$ 위에 있는 Joint Distribution 를 의미한다고 했습니다. 그러면 Gaussian Copula 란, 다음과 같은 $(X,Y)$ 가 있다고 하겠습니다.

그러면 $\Phi(X)$ 와 $\Phi(Y)$ 의 Joint Distribution 을 Gaussian Copula 이라고 부릅니다. 예시로 $\rho = 0.5$ 일 때의 Gaussian Copula 를 뽑겠습니다.

이런 방식으로 Copula 에서 뽑는 것이 가능합니다. 위의 코드에서는 사실 다음과 같이 흘러갑니다.

• 먼저 $X$ 를 뽑아준다.

• 그 다음에 $Y$ 를 뽑는데, 이때 $Y = \rho \sigma_2 Z_1 + \sigma_2 \sqrt{1-\rho^2}Z_2 + \mu_2$ 인 점을 활용합니다.

• $U_1 := F(X)$ 으로 세팅합니다.

• $U_2 := G(Y)$ 으로 세팅합니다.

위의 방법이 사실상 Copula 부분에서 설명한 것이랑 정확히 똑같습니다.

이 방법을 활용시켜서 다음과 같은 예시에 적용해보겠습니다.

예시 문제

$X \sim \mathrm{Exp}(3), \; Y \sim \Gamma(2, 2)$ 이며 $\mathrm{Corr}(X, Y) = \rho = 0.5$ 일 때 $(X,Y)$ 를 뽑아서 Plot 을 그려주세요.

1. 우선은 $\rho = 0.5$ 인 Gaussian Copula 에서 뽑아줍니다. 그러면 $(U_1, U_2)$ 들이 있을 것이며 이들의 상관관계는 $0.5$ 일 것입니다.

2. $X:= F^{-1}(U_1), \; Y:= G^{-1}(U_2)$ 으로 세팅해주세요.

3. Plot 시켜줍니다.

Marginal 들을 그려주면 다음과 같습니다.

$X$ 와 $Y$ 에 대한 Plot 을 그려주면 다음과 같습니다.

N = 10000
rho = -0.5

x = qnorm(runif(N))
z_2 = qnorm(runif(N))
y = x * rho + sqrt(1 - rho^2) * z_2

plot(x,y, cex = 0.2, pch = 20) # Joint Distribution of Gaussian

u_1 <- pnorm(x)
u_2 <- pnorm(y)

plot(u_1, u_2, cex = 0.2, pch = 20) # Gaussian Copula

new_x <- qexp(u_1, 3)
new_y <- qgamma(u_2, 2, scale = 2)

plot(new_x, new_y, cex = 0.2, pch = 20) # Joint Distribution

hist(new_x, breaks = 100)
hist(new_y, breaks = 100)

이것이 가능한 이유는 사실은 다음과 같이 Copula 에 Monte-Carlo 를 적용시키는 것이 가능하기 때문입니다.

만약에 Copula $C$ 가 Gaussian Copula 와 같이 뽑기가 쉽다면 저희가 위에서 한 것과 같이 진행하면 됩니다.

• 번외

  Marshall-Olkin Copula 란 그저 Copula 를 Multivariate Exponential 에 적용시킨 Copula 입니다.