6장 Variance Reduction Techniques

  • Antithetic Variables
  • Control Variables
  • Conditioning
  • Stratified Sampling
  • Importance Sampling

여태까지는 SRS 에서 뽑는 것이었다. 즉 다시 말해 뽑는 과정이 그저 i.i.d 하게 뽑는 것만 고려를 했었다. 하지만 이제는 i.i.d 하게 뽑기보다는 조금 다르게 뽑아줘서 추정량에 대한 분산을 줄여줘보겠습니다.

정확하게 그래서 어떤 얘기를 하느냐, 우선은 예시로 Xˉ\bar X 를 떠올려보겠습니다.

θ=E[X]\theta = \mathbb{E}[X] 을 추정하기 위해서는 Xˉ\bar X 를 사용합니다. 이때 Xˉ\bar X 는 Unbiased + Consistent 합니다. 근데 MSE 를 가지고 Xˉ\bar X 가 얼마나 좋은지에 대한 지표도 구할 수가 있습니다.

MSE=E[(Xˉ−θ)2]→Var(Xˉ)=Var(X)n\text{MSE} = \mathbb{E}\left[ \left( \bar X - \theta \right)^2 \right] \rightarrow \mathrm{Var}(\bar X) = \frac{\mathrm{Var}(X)}{n}

그저 i.i.d 하게 뽑는 것만 고려한다면 샘플 개수인 nn 을 올려서 추정량에 대한 분산, 즉 표준오차를 줄일 수가 있습니다. 하지만 시뮬레이션에서는 이것을 줄일 수 있는 여러 방법들이 존재합니다.

Antithetic Variables

예전에 분산을 줄이기 위해서 사용되었던 방법이며, Negative 하게 뽑아서 Negative Correlation 인 점을 활용합니다.

X1,X2X_1, X_2 를 i.i.d 하게 뽑아준다면 다음과 같이 나옵니다.

Var(X1+X22)=14[Var(X1)+Var(X2)+2⋅Cov(X1,X2)]\mathrm{Var}\left(\frac{X_1 + X_2}{2} \right) = \frac 1 4 \left[ \mathrm{Var}(X_1) + \mathrm{Var}(X_2) + 2 \cdot \mathrm{Cov}(X_1, X_2)\right]

하지만 i.i.d 하기 때문에 당연하게도 공분산이 00 입니다. 그렇기에 분산은 사실 다음과 같습니다.

Var(X1+X22)=14[Var(X1)+Var(X2)]\mathrm{Var}\left(\frac{X_1 + X_2}{2} \right) = \frac 1 4 \left[ \mathrm{Var}(X_1) + \mathrm{Var}(X_2) \right]

여기에서 분산을 그러면 어떻게 줄이냐, 바로 Cov(X1,X2)\mathrm{Cov}(X_1, X_2) 를 음수로 바꿔주는 것입니다. 가장 간단하게는 X1=h(U1,⋯ ,Um)X_1 = h(U_1, \cdots, U_m) 이면 X2=h(1−U1,1−U2,⋯ ,1−Um)X_2 = h\left(1-U_1, 1-U_2, \cdots, 1-U_m \right) 으로 바꿔주는 것입니다.

이렇게 진행한다면 Cov(X1,X2)<0\mathrm{Cov}(X_1, X_2) < 0 이기에 다음이 성립하게 됩니다.

Var(X1+X22)=14[Var(X1)+Var(X2)+2⋅Cov(X1,X2)]<14[Var(X1)+Var(X2)]\begin{aligned} \mathrm{Var}\left(\frac{X_1 + X_2}{2} \right) &= \frac 1 4 \left[ \mathrm{Var}(X_1) + \mathrm{Var}(X_2) + 2 \cdot \mathrm{Cov}(X_1, X_2)\right] \\&< \frac 1 4 \left[ \mathrm{Var}(X_1) + \mathrm{Var}(X_2) \right]\end{aligned}

간단한 예시로 생각해보겠습니다.

Example 1

θ=E[eU]\theta = \mathbb{E}\left[e^{U} \right] 를 추정하고 싶다고 하겠습니다. 우선은 분산, 공분산에 대한 값들을 다음과 같이 구하겠습니다.

Var(eU)=E[(eU)2]−(E[eU])2=∫01e2udu−([ex]x=0x=1)2=[12⋅e2x]x=0x=1−(e−1)2=e2−12−(e−1)2≈0.242 Cov(eU,e1−U)=E[eU⋅e1−U]−E[eU]⋅E[e1−U]=E[eU⋅e1⋅e−U]−(e−1)⋅∫01e1−udu=E[e]−(e−1)⋅e⋅[−e−x]x=0x=1=e−(e−1)⋅e⋅(1−e−1)=e−(e−1)⋅e⋅(e−1)e=e−(e−1)2≈−0.234\begin{aligned} \mathrm{Var}\left(e^U\right) &= \mathbb{E}\left[\left(e^{U}\right)^2 \right] - \left(\mathbb{E}\left[e^U \right]\right)^2 = \int_0^1 e^{2u} du - \left(\left[ e^x\right]_{x=0}^{x=1} \right)^2 \\&= \left[\frac 1 2 \cdot e^{2x} \right]_{x=0}^{x=1} - \left( e - 1 \right)^2 = \frac{e^2 - 1}{2} - \left( e - 1 \right)^2 \approx 0.242 \\~\\ \mathrm{Cov}\left(e^U, e^{1-U} \right) &= \mathbb{E}\left[e^U \cdot e^{1-U} \right] - \mathbb{E}\left[e^U \right]\cdot \mathbb{E}\left[e^{1-U}\right] \\&= \mathbb{E}\left[ e^U \cdot e^1 \cdot e^{-U} \right] - (e-1)\cdot \int_0^1 e^{1-u} du \\&= \mathbb{E}\left[e \right] - (e-1) \cdot e \cdot \left[-e^{-x} \right]_{x=0}^{x=1} \\&= e - (e-1)\cdot e \cdot (1-e^{-1}) \\&= e - (e-1) \cdot e \cdot \frac{(e-1)}{e} \\&= e - (e-1)^2 \approx -0.234 \end{aligned}

일반적인 Monte-Carlo Estimate 방법에서는 그저 U1,U2,⋯ ,UnU_1, U_2, \cdots, U_n 까지 뽑은 다음에 이것을 eUie^{U_i} 에 넣어준 다음에 1n∑i=1neUi\frac{1}{n}\sum_{i=1}^n e^{U_i} 으로 평균을 내서 추정량을 만들어주겠죠.

22 개의 샘플에 대해서는 다음과 같이 분산을 구해주겠습니다.

Var(eU1+eU22)=12⋅Var(eU)≈0.121\mathrm{Var}\left(\frac{e^{U_1} + e^{U_2}}{2} \right) = \frac{1}{2} \cdot \mathrm{Var}(e^{U}) \approx 0.121

이제는 Antithetic Variable 를 활용해보겠습니다.

똑같이 U1,(1−U1)U_1, (1-U_1) 을 사용해보겠습니다. 그러면 eU1+e1−U12\frac{e^{U_1} + e^{1-U_1}}{2} 가 저희의 추정량이며 이것의 분산은 다음과 같습니다.

Var(eU1+e1−U12)=14⋅[Var(eU1)+Var(e1−U1)+2⋅Cov(eU1,e1−U1)]=12⋅[Var(eU1)+Cov(eU1,e1−U1)]=12⋅[0.242−0.234]≈0.004\begin{aligned} \mathrm{Var}\left(\frac{e^{U_1} + e^{1 - U_1}}{2} \right) &= \frac{1}{4} \cdot \left[ \mathrm{Var}\left(e^{U_1} \right) + \mathrm{Var}\left(e^{1-U_1}\right) + 2 \cdot \mathrm{Cov}\left(e^{U_1}, e^{1-U_1} \right) \right] \\&= \frac{1}{2}\cdot \left[\mathrm{Var}\left(e^{U_1} \right) + \mathrm{Cov}\left(e^{U_1}, e^{1-U_1} \right) \right] \\&= \frac 1 2 \cdot \left[ 0.242 - 0.234 \right] \approx 0.004 \end{aligned}

두 분산을 비교하면, 확실하게 이득을 본 모습입니다.

사실은 Monte-Carlo Estimate 에서 22 번 뽑을 것을 저희는 11 번만 뽑아서 사용하고 있습니다. 한 개로 두 개 효과를 보는 것이죠. 그렇기에 위 예시에서는 반만큼의 샘플들만 필요하다는 점도 있습니다.

이것에 대한 실제 코드와 함께 확인해보겠습니다.

rm(list=ls())
N = 10000
mc_estimator = exp(runif(N))
u_vals = runif(N / 2)
antithetic_estimator = c((exp(u_vals) + exp(1 - u_vals)) / 2)
cat("MC Estimator Mean : ", mean(mc_estimator), "Var : ", var(mc_estimator),
    "\nAntithetic Estimator : ", mean(antithetic_estimator), "Var : ", var(antithetic_estimator))

Output :

Example 2

원에서 π\pi 구하는 예시로 돌아가보겠습니다. 이것을 저희는 이전에는 그저 (U1,U2)(U_1, U_2) 들을 NN 개 만큼 뽑은 다음에 1(U12+U22≤1)\mathbb{1}\left(U_1^2 + U_2^2 \le 1 \right) 를 추정량으로 세웠습니다. 이게 Naive 한 MC Estimator 이지만, 이제는 다르게 (1−U1,1−U2)(1-U_1, 1-U_2) 까지 활용해보겠으며, 코드로는 다음과 같습니다.

rm(list=ls())
# MC-Estimate
N = 10000
u1_vals = runif(N) 
u2_vals = runif(N)
circle_indicator = ifelse(u1_vals ** 2 + u2_vals ** 2 <= 1, 4, 0)
mc_est = mean(circle_indicator)
mc_var = var(circle_indicator)
mc_SE = sqrt(mc_var / N)

# Antithetic Estimator
M = N / 2
u1_vals = runif(M)
u2_vals = runif(M)
u1_u2_est = ifelse(u1_vals ** 2 + u2_vals ** 2 <= 1, 4, 0)
anti_u1_u2_est = ifelse((1 - u1_vals) ** 2 + (1 - u2_vals) ** 2 <= 1, 4, 0)

anti_circle_est = (u1_u2_est + anti_u1_u2_est) / 2

anti_est = mean(anti_circle_est)
anti_var = var(anti_circle_est)
anti_SE = sqrt(anti_var / M)

cat("MC Estimate :\n", "Mean : ", mc_est, "\tVar : ", mc_var, "\tSE : ", mc_SE,
    "\nAntithetic Estimate :\n", "Mean : ", anti_est, "\tVar : ", anti_var, "\tSE : ", anti_SE)

드라마틱한 감소는 없지만, 표준오차에 대한 감소가 있기는 합니다.

Control Variables

똑같이 θ=E[X]\theta = \mathbb{E}\left[X\right] 을 추정하는게 목표지만, 이번에는 다음과 같은 YY 라는 새로운 변수를 활용하겠습니다.

X+c(Y−μY)X + c\left(Y - \mu_Y\right)
  • E[X+c(Y−μY)]=E[X]=θ\mathbb{E}\left[X + c\left(Y - \mu_Y\right) \right] = \mathbb{E}[X] = \theta, Unbiased 하다.
  • Var(X+c(Y−μY))=Var(X+cY)=Var(X)+c2⋅Var(Y)+2c⋅Cov(X,Y)\mathrm{Var}\left(X + c\left(Y - \mu_Y\right) \right) = \mathrm{Var}\left(X + cY \right) = \mathrm{Var}(X) + c^2 \cdot \mathrm{Var}(Y) + 2 c \cdot \mathrm{Cov}(X, Y)

분산 부분을 살피면, 저희는 cc 를 어떻게 고르냐에 따라서 추정량의 분산을 줄일지 올릴지 정해줄 수가 있습니다.

c∗:=arg min⁡c[Var(X)+c2⋅Var(Y)+2c⋅Cov(X,Y)]c^* := \argmin_c\left[ \mathrm{Var}(X) + c^2 \cdot \mathrm{Var}(Y) + 2 c \cdot \mathrm{Cov}(X, Y) \right]

이것을 구하기 위해서는 cc 에 대해서 미분을 시켜주면 됩니다.

2c⋅Var(Y)+2⋅Cov(X,Y)=set0c∗=−Cov(X,Y)Var(Y)\begin{aligned} 2c \cdot \mathrm{Var}(Y) &+ 2 \cdot \mathrm{Cov}(X, Y) \stackrel{\text{set}}{=} 0 \\ c^* & = -\frac{\mathrm{Cov}(X, Y)}{\mathrm{Var}(Y)} \end{aligned}

cc 를 위와 같이 대입시킨다면, 다음과 같이 분산 Term 이 바뀝니다.

Var(X−Cov(X,Y)Var(Y)⋅(Y−μY))=Var(X)+[Cov(X,Y)Var(Y)]2⋅Var(Y)−2⋅[Cov(X,Y)Var(Y)]⋅Cov(X,Y)=Var(X)−[Cov(X,Y)]2Var(Y)=Var(X)−ρ2⋅Var(X)=Var(X)⋅(1−ρ2)≤Var(X)\begin{aligned} \mathrm{Var}\left(X - \frac{\mathrm{Cov}(X, Y)}{\mathrm{Var}(Y)} \cdot \left(Y - \mu_Y\right) \right) &=\mathrm{Var}(X) + \left[\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(Y)}\right]^2 \cdot \mathrm{Var}(Y) \\& - 2 \cdot \left[\frac{\mathrm{Cov}(X,Y)}{\mathrm{Var}(Y)} \right] \cdot \mathrm{Cov}(X, Y) \\ & = \mathrm{Var}(X) - \frac{\left[\mathrm{Cov}(X,Y) \right]^2}{\mathrm{Var}(Y)} \\&= \mathrm{Var}(X) - \rho^2 \cdot \mathrm{Var}(X) \\&= \mathrm{Var}(X) \cdot \left( 1 - \rho^2\right) \le \mathrm{Var}(X) \end{aligned}

Example 3

이전과 똑같은 예시를 가지겠습니다. θ=E[eU]\theta = \mathbb{E}\left[e^{U} \right] 를 Y=UY = U 으로 세팅하겠습니다.

rm(list=ls())
N = 10000
u_vals = runif(N)
exp_u_vals = exp(u_vals)
c = - cov(u_vals, exp_u_vals) / var(u_vals)
mu_y = 0.5
control_estimators = exp_u_vals + c * (u_vals - mu_y)
mc_mean = mean(exp_u_vals)
mc_var = var(exp_u_vals)
mc_SE = sqrt(mc_var / N)
control_mean = mean(control_estimators)
control_var = var(control_estimators)
control_SE = sqrt(control_var / N)
cat("MC Estimators :\nMean :", mc_mean, "\tVar :", mc_var, "\tSE :", mc_SE,
    "\nControl Estimators : \nMean :", control_mean, "\tVar :", control_var, "\tSE :", control_SE)

표준오차가 줄어든 모습입니다.

Conditioning

EVVE 를 활용하면 다음과 같이 분산을 적어줄 수가 있습니다.

Var(X)=E[Var(X∣Y)]+Var(E[X∣Y])≥Var(E[X∣Y])\begin{aligned} \mathrm{Var}(X) & =\mathbb E[\mathrm{Var}(X \mid Y)]+\mathrm{Var}(\mathbb E[X \mid Y]) \\ & \geq \mathrm{Var}(\mathbb E[X \mid Y]) \end{aligned}

이 점을 활용해서 XX 대신에 E[X∣Y]\mathbb{E}\left[X \mid Y\right] 으로 추정을 해주면 됩니다.

  • E[E[X∣Y]]=E[X]=θ\mathbb{E}\left[\mathbb{E}\left[X \mid Y \right] \right] = \mathbb{E}\left[ X \right] = \theta, Unbiased
  • Var(E[X∣Y])≤Var(X)\mathrm{Var}\left(\mathbb{E}[X \mid Y] \right) \le \mathrm{Var}(X)

그렇기에 저희는 사실은 추정량을 E[X∣Y]=g(Y)\mathbb{E}[X \mid Y] = g(Y) 으로 사용하고 있습니다.

Example 4

(U1,U2)∼Unif⁡([0,1]2)\left(U_1, U_2\right) \sim \operatorname{Unif}\left([0,1]^2\right) 이며 X=4⋅1{U12+U22≤1}X=4 \cdot \mathbf{1}\left\{U_1^2+U_2^2 \leq 1\right\} 일 때, E[X]=π\mathbb{E}[X] = \pi 입니다.

기본적인 MC-Estimate 을 사용한다면 분산은 다음과 같이 나옵니다.

Var⁡(X)=16Var⁡(I)=16p(1−p)=16(π4)(1−π4)=4π−π2\operatorname{Var}(X)=16 \operatorname{Var}(I)=16 p(1-p)=16\left(\frac{\pi}{4}\right)\left(1-\frac{\pi}{4}\right)=4 \pi-\pi^2

대신에 다음과 같은 Condition 을 활용하겠습니다.

E[X∣U1=u1]=4⋅P(U2≤1−u12)=4⋅1−u12\mathbb{E}\left[X \mid U_1=u_1\right]=4\cdot \mathrm{P}\left(U_2 \leq \sqrt{1-u_1^2}\right)=4 \cdot \sqrt{1-u_1^2}

그러면 이를 활용해서 g(U1)=4⋅1−U12g(U_1) = 4\cdot \sqrt{1-U_1^2} 를 사용하면 됩니다. 이때 U2U_2 를 한번 더 샘플링 할 필요도 없으며, 분산은 다음과 같은 이유로 작아지게 됩니다.

E[g(U)]=4∫011−u2du=4⋅π4=π\mathbb{E}[g(U)]=4 \int_0^1 \sqrt{1-u^2} d u=4 \cdot \frac{\pi}{4}=\pi
E[g(U)2]=16∫01(1−u2)du=16[u−u33]01=16(1−13)=323\mathbb{E}\left[g(U)^2\right]=16 \int_0^1\left(1-u^2\right) d u=16\left[u-\frac{u^3}{3}\right]_0^1=16\left(1-\frac{1}{3}\right)=\frac{32}{3}
Var(E[X∣U1])=E[g(U)2]−(E[g(U)])2=323−π2\mathrm{Var}\left(\mathbb{E}\left[X \mid U_1 \right]\right) = \mathbb{E}\left[g(U)^2 \right] - \left(\mathbb{E}\left[g(U) \right]\right)^2 = \frac{32}{3} - \pi^2

위의 것과 비교하면 4π>32/34\pi > 32/3 이기에 분산이 줄여진 모습을 볼 수가 있습니다.

Stratified Sampling

Conditioning 개념을 이제는 Group 들로 나누어서 가는 개념입니다. 즉 만약에 YY 가 Discrete 한 확률변수이면, pi=Pr{Y=yi},i=1,⋯ ,kp_i = \mathrm{Pr}\{Y=y_i\}, i = 1, \cdots, k 를 저희는 알고 있다고 하겠습니다. 그러면 다음과 같이 바꿔줄 수가 있습니다.

E[X]=E[E[X∣Y]]=∑i=1kE[X∣Y=yi]⋅pi\mathbb{E}[X] = \mathbb{E}[\mathbb{E}[X\mid Y]] = \sum_{i=1}^k \mathbb{E}\left[X \mid Y = y_i \right] \cdot p_i

저희의 최종 추정량은 그럼 θ^=∑i=1kpiXˉi\hat \theta = \sum_{i=1}^k p_i \bar X_i 이며, 분산은 다음과 같습니다.

Var⁡(θ^)=∑i=1kpi2niVar⁡(X∣Y=yi)=1n∑i=1kpiVar⁡(X∣Y=yi)( if we let ni=npi)=1nE[Var⁡(X∣Y)]≤1nVar⁡(X)\begin{aligned}\operatorname{Var}(\hat{\theta}) & =\sum_{i=1}^k \frac{p_i^2}{n_i} \operatorname{Var}\left(X \mid Y=y_i\right) \\& =\frac{1}{n} \sum_{i=1}^k p_i \operatorname{Var}\left(X \mid Y=y_i\right) \quad\left(\text { if we let } n_i=n p_i\right) \\& =\frac{1}{n} E[\operatorname{Var}(X \mid Y)] \leq \frac{1}{n} \operatorname{Var}(X)\end{aligned}

Importance Sampling

샘플링하기 어려운 것들을 더 쉽게 샘플링하는 방법입니다. 예시로 θ=Ef[h(X)]=∫h(x)f(x)dx\theta = \mathbb{E}_f[h(X)] = \int h(x) f(x) dx 에서 샘플하고 싶은데, ff 에 대한 분포를 샘플하는게 어렵거나 h(x)h(x) 에 대한 분산이 매우 높다고 하겠습니다. 그러면 다른 쉬운 분포 gg 를 사용해서 뽑아줍니다.

θ=Ef[h(X)]=∫h(x)f(x)dx=∫h(x)g(x)g(x)f(x)dx=∫g(x)h(x)g(x)f(x)dx=Eg[h(x)g(x)⋅f(x)]≈1n∑h(xi)f(xi)g(xi)\begin{aligned} \theta &= \mathbb{E}_f[h(X)]= \int h(x) f(x) dx = \int h(x) \frac{g(x)}{g(x)} f(x) dx \\&= \int g(x) \frac{h(x)}{g(x)} f(x) dx = \mathbb{E}_g\left[\frac{h(x)}{g(x)} \cdot f(x)\right] \approx \frac{1}{n}\sum \frac{h(x_i)f(x_i)}{g(x_i)} \end{aligned}