• Inverse-Transform Method

• Acceptance-Rejection Method

• Transformation of Random Variables

1장에서 얘기한 바와 같이, 저희는 이미 어떻게하면 $U\stackrel{\mathrm{i.i.d}}{\sim} \mathrm{Unif}(0, 1)$ 에서 뽑을지 알고 있습니다. 하지만 이것을 활용해서 저희가 원하는 분포에서 뽑으려면 어떻게 해야 할까요? 그것에 대한 가장 간단한 1-차원 수준에서 뽑는 것이 2장에서 얘기하는 바입니다. 이것을 해결해주는 아주 간단한 방법이 3가지가 있으며, 그것들은 바로 Inverse-Transform Method, Acceptance-Rejection Method, 그리고 Transformation of Random Variables 가 있습니다.

Inverse-Transform Method

어떠한 Non-decreasing 함수 $h : \mathcal{X} \rightarrow \mathcal{Y}$ 가 있다고 하겠습니다. 이 함수의 Inverse 함수의 General Form 은 다음과 같이 나오게 됩니다.

이게 의미하는 바는 무엇이냐면, 어떠한 Non-decreasing 함수 $h(x)$ 가 있으며 이것은 때로는 Jump 가 생길 수도 있습니다. 그 점들을 고려하기 위해서 $\inf$ 을 넣어주게 됩니다. 아래 사진을 보겠습니다.

여기에서 $z_2$, $z_3$, $z_4$ 에 해당되는 $x$ 값들은 다 똑같습니다. 즉 $h^{-1}(z_2) = h^{-1}(z_3) = h^{-1}(z_4)$ 입니다. 해석을 하자면, $z_3$ 보다 큰 $h(x)$ 값은 $z_4$ 이상 밖에 없습니다. 그렇기에 그것의 Local 최솟값 인 $\inf$ (물론 $\inf$ 의 정의는 다르긴 하지만) 으로 잡는다면 $z_3$ 에 대한 Inverse 값을 $h^{-1}(z_4)$ 으로 둘 수 있습니다.

이제는 이 $h$ 함수를 cdf 함수 $F$ 으로 바라보겠습니다. $F$ 의 경우 Non-decreasing 하다는 특징이 있으며, 만약에 $F$ 에 대한 역함수가 존재한다면, 다음과 같이 세팅해줄 수 있습니다.

즉 다시 말하자면 $X$ 에 대해서 뽑기 위해서 저희가 쉽게 뽑을 수 있는 균일 확률변수 $U\sim \mathrm{Unif}(0, 1)$ 에서 뽑은 뒤 그 값을 cdf 의 Inverse 함수에 넣으면 $X$ 에서 뽑는 것과 똑같아 지는 효과를 줍니다. 즉 다시 말해,
$X = F^{-1}(U)$ 으로 세팅하게 된다면 다음과 같이 나오게 됩니다.

그러면 직접 $X\sim F$ 를 뽑기 위해서는 어떻게 해야 할까요? 바로 다음을 따라가기만 하면 됩니다.

1. $U \sim \mathrm{Unif}(0, 1)$ 을 뽑는다.

2. $X = F^{-1}(U)$ 으로 세팅한다.

예시 (정규분포) :

정규분포의 Inverse cdf 함수는 $\texttt{qnorm}$ 입니다. 이것을 활용시켜서 $X\sim\mathrm{N}\left(1, \sqrt{3} \right)$ 에서 샘플하기 위해서 다음과 같은 코드를 따라가면 됩니다.

U_samples <- runif(10000)
X_samples <- qnorm(U_samples, mean = 1, sd = 3)
hist(X_samples, breaks = 100, col = 'lightblue', freq = FALSE)
curve(dnorm(x, mean = 1, sd = 3), add = TRUE,
      col = 'pink', lwd = 4)
legend('topright', legend = c('Samples from Inverse cdf', 'pdf of N(1, 3)'),
       col = c('lightblue', 'pink'), lwd = 2)

이렇게 하게 된다면 결과는 다음과 같은 Plot 처럼 나오게 됩니다.

예시 (지수 분포) :

지수 분포를 따르는 확률변수 $X$ 의 cdf 는 다음과 같이 생겼습니다.

이 함수의 Inverse 는 다음과 같이 생겼습니다.

이것을 활용해서 $X \sim \mathrm{Exp}(3)$ 에서 샘플하도록 하겠습니다.

U_samples <- runif(10000)
X_samples <- qnorm(U_samples, mean = 1, sd = 3)
hist(X_samples, breaks = 100, col = 'lightblue', freq = FALSE)
curve(dnorm(x, mean = 1, sd = 3), add = TRUE,
      col = 'pink', lwd = 4)
legend('topright', legend = c('Samples from Inverse cdf', 'pdf of N(1, 3)'),
       col = c('lightblue', 'pink'), lwd = 2)

결과는 다음과 같이 나오게 됩니다.

예시 (이산 확률변수) :

가방에서 공들이 똑같은 개수로 썪여 있으며, 각각의 색깔을 뽑을 확률이 다음 테이블을 따른다고 하겠습니다.

그러면 이때 이 확률변수의 cdf 는 다음과 같이 나오게 됩니다.

여기에서 그러면 이제 Inverse CDF 를 활용하도록 하겠습니다. $U \sim \mathrm{Unif}(0, 1)$ 에서 이미 뽑았다고 하면, 이 값에 따라서 $X$ 를 다음과 같이 부여시키면 됩니다.

만약에 저희가 $U^\prime \sim \mathrm{Unif}(0, 4)$ 에서 뽑았다면, 이것을 다음과 같이 바꿔줄 수 있습니다.

이때 $U \sim \mathrm{Unif}(0, 1)$ 이라고 하면 $X$ 를 $X = \left\lfloor4U\right\rfloor + 1$ 으로 바꿔줄 수 있게 됩니다.

• 예시로 만약에 $u = 0.2$ 으로 나왔다면 $\left\lfloor4 \cdot 0.2\right\rfloor + 1 = \left\lfloor0.8\right\rfloor + 1 = 0 + 1 = 1$

• 만약에 $u = 0.4$ 으로 나왔다면 $\left\lfloor4 \cdot 0.4\right\rfloor + 1 = \left\lfloor1.6\right\rfloor + 1 = 1 + 1 = 2$

이것을 코드로 짜게 된다면 다음과 같습니다.

U_samples <- runif(10000)
X_samples <- floor(U_samples * 4) + 1
# 여기서부터는 시각화라 무시해도 됩니다.
prop <- prop.table(table(X_samples))
barplot(prop,
        names.arg = 1:4,
        ylim     = c(0, 0.35),
        ylab     = "Empirical P(X = k)",
        xlab     = "k",
        main     = "Inverse Transform Method",
        col = 'lightblue'
        )
abline(h = 0.25, col = "pink", lty = 2, lwd = 2)
legend("topright",
       legend = "P = 0.25",
       col    = "pink",
       lty    = 2,
       lwd    = 2,
       bty    = "n")

이에 해당되는 결과는 다음과 같습니다.

Random Permutation

가끔은 순서를 랜덤하게 바꾸어야 합니다. 예시로 저희는 Permutation $y_1, y_2, \cdots, y_n$ 이 있다고 하겠습니다. 저희의 목표는 이 순서를 완전히 랜덤하게 바꾸는 것이며, 이는 간단하게 이해시켜드리기 위해 $y_1, y_2, y_3, y_4$ 으로 설명해드리겠습니다.

그전에 설명을 하자면, $x$ 는 색깔이 $k$ 개인 위의 예시에서의 이산형 확률변수 (정확하게는 Categorical 변수) 라고 하겠습니다.

위를 Generalize 시키면, 다음과 같습니다.

1. Permutation $A = \left[y_1, y_2, \cdots, y_n\right]$ 으로 시작하며 $k = n$ 으로 세팅한다.

2. $U \sim \mathrm{Unif}(0, 1)$ 을 뽑고, $I = \left\lfloor kU \right\rfloor + 1$ 을 구합니다.

3. $A[I] \leftrightarrow A[k]$ 즉 둘의 순서를 바꿉니다.

4. $k = k-1$ 으로 업데이트 시키며 $2, 3$ 을 $k = 1$ 이 도달할 때까지 반복합니다.

예시로 위의 Permutation 을 $(1, 2, 3, 4, 5, 6, 7,8,9)$, 즉 $n = 9$ 에서 진행하겠습니다. 그러면 Random Permutation 이 성공했다면 $1$-번째 자리에서 나올 값이 List 의 평균인 $4.5$ 와 동일하겠죠. 이를 한번 코드로 실행해보겠습니다.

random_permutation <- function(x){
  n = length(x)
  for(k in n:2){
    u <- runif(1)
    I <- floor(k * u) + 1
    I_th <- x[I]
    x[I] <- x[k]
    x[k] <- I_th
    k = k - 1
  }
  return(x)
}


x <- c(1, 2, 3, 4, 5, 6, 7, 8, 9)
N <- 10000
x_1st_list <- rep(NA, N)
for(i in 1:N){
  x_1st_list[i] = random_permutation(x)[1]
}
cat("X 의 평균 : ", mean(x),
    "\nPermutation 의 평균 : ", mean(x_1st_list))

그러면 결과는 다음과 같이 나옵니다.

Acceptance-Rejection Method

우선은 Discrete 인 상황을 고려해보겠습니다. 예시로 아까 보았던 이산형 확률변수를 생각해보겠습니다.

이때 ${}^{\forall}y\in \{1, 2, 3, 4 \},\; q(Y=y) = 0.25$ 이라는 비교적 쉬운 분포 (위에 보았던 예시) 를 고려해보겠습니다. 그러면 $X$ 를 뽑기 위해서는 다음과 같은 Rejection Method Algorithm 으로 뽑아주면 됩니다.

1. $U_1 \sim \mathrm{Unif}(0, 1)$

2. $U_2 \sim \mathrm{Unif}(0, 1)$, $Y = \mathrm{Int}(4 U_2) + 1$

3. If $U_1 \le p(X = Y)$, $X = Y$ 으로 세팅해준 뒤에 멈춘다.

이것을 그림과 같이 해석을 해보겠습니다.

그림에 보이시는 바와 같이 그저 파랑색 쪽에 대해서 샘플링을 하는 것이 목표입니다. 이것을 이루기 위해서 보다 더 쉬운 분포 $Y\sim q_Y$ 에서 샘플링 해준 다음에 그 값이 $p_X$ 위에서는 파랑색 영역에 들어가는지 아닌지에 따라서 샘플링을 해줍니다.

위와 같은 알고리즘으로 진행하게 된다면 당연하게도 굉장히 비효율적입니다. for-loop 을 각 $Y$ 샘플에 대해서 사용해야하기 때문입니다. 심지어는 가장 높은 Threshold 까지 뽑아도 되지만 굳이 그 위까지 샘플링하고 있죠.

그렇기에 저희는 다음과 같은 상수를 규명하겠습니다.

사실 위의 알고리즘의 3 번은 $U_1 \le \frac{p(X=Y)}{c\cdot q(Y=Y)}$ 와 똑같습니다. $q_Y$ 가 그저 $\frac 1 4$ 의 균일분포이기에 모든 $Y$ 값에 대해서 $0.25$ 가 나오며, 이 값에 $c = 4$ 를 세팅해준다면 위의 알고리즘과 같은 조건이 되며, 이것을 $c = \max_x\frac{p(X=x)}{q(Y=x)}$ 으로 세팅한다면 $q_Y$ 분포를 최대한 파랑색에 당기게끔 만들어줍니다. 예시로 위의 똑같은 시뮬레이션을 $c = 1.8$ 일 때 다시 해보겠습니다.

이렇게 하게 된다면 버리는 샘플들을 최소화 시킬 수가 있게 됩니다. 이때 $c$ 의 값이 $1$ 에 가까워질 수록 Rejection Method 가 더 빨라지며, $c$ 만큼 뽑으면 $1$번 $p_X$ 에서 뽑는 효과를 가집니다.

$q_Y$ 를 일종의 $p_X$ 를 둘러싸는 쉬운 분포라고 생각해주시면 됩니다. 연속적일 때도 비슷하며, 이는 다음과 같은 알고리즘을 띄우게 됩니다.

목표 : $X \sim F$, Proposal $Y \sim G$

1. 먼저 $c = \sup_x \frac{f(x)}{g(x)}$ 세팅부터 한다.

2. $Y\sim G$, $U\sim \mathrm{Unif}(0, 1)$ 을 샘플한다.

3. 만약에 $U \le \frac{f(Y)}{c\cdot g(Y)}$ 이면, $X=Y$ 으로 세팅하면서 뽑으며, 그게 아니라면 2번으로 다시 돌아갑니다.

예시로 Mixture of 정규분포에서 뽑아보겠습니다.

이것의 Density Plot 은 다음과 같이 나옵니다.

이때 이 분포 내에서 한번 샘플을 뽑아 달라고 했을 때, Rejection Method 를 통해서 쉽게 수행하는 것이 가능합니다. Proposal Density 의 경우 $g(x) = \texttt{dnorm}(x; 0, 4)$ 으로 설정해두겠습니다. 그러면 다음과 같이 Rejection Method 를 거리에서 하면 다음과 같은 결과가 나옵니다.

1. $c = \max_x \frac{f(x)}{g(x)} = 4.418967$

2. $Y \sim G, \; U_2 \sim \mathrm{Unif}(0, 1), \text{ where } Y, U_2 \in \mathbb{R}^{N}$

3. $\texttt{ifelse}(U_2 \le \frac{f(Y)}{c \cdot g(Y)}, Y, \text{Skip})$

이것을 시각화 시키게 된다면 다음과 같이 나옵니다.

위에 보이시는 바와 같이 $\frac{1}{c}$ 정도의 샘플들을 버리긴 하지만 결국에는 파랑색 분포 (Target) 에서 샘플링하는 효과를 해주게 됩니다.

Transformation of Random Variables

마지막으로는 Transformation of R.V’s 입니다. 직관적인 예시를 위해서 $\mathcal{X}^2(3)$ 에서 뽑고 싶다고 하겠습니다. 저희는 $Z \sim N(0, 1)$ 에서 어떻게 뽑는지 Inverse Transform Method 를 통해서 이미 알고 있습니다. 그렇기에 표준정규분포와 카이제곱 분포가 가지는 관계를 이용해서 $\mathcal{X}^2 (3)$ 에서 뽑아보겠습니다.

• $Z^2 \sim \mathcal{X}(1)$

• $Z_1, Z_2, Z_3 \stackrel{\mathrm{i.i.d}}{\sim} N(0,1)$

• $V = \sum_{i=1}^3 Z_i^2 \sim \mathcal{X}(3)$

$V$ 를 시각화 시켜서 원래 $\mathcal{X}^2(3)$ 의 pdf 와 비교해서 그린다면 다음과 같이 나오게 됩니다.

이처럼 표준정규분포에서 쉽게 뽑을 수 있기에 함수들을 활용시켜서 다른 분포들에게서도 뽑을 수 있게 됩니다. 이것도 일종의 Inverse Transform Method 이라고 부르기도 합니다.