Redução de Variância por Variáveis Antitéticas

Existem várias técnicas para reduzir a variância de um estimador de Monte Carlo.

Uma técnica clássica de redução de variância em simulação de Monte Carlo é baseada em variáveis antitéticas (antithetic variates). Vamos visitar esta técnica através de um exemplo. A ideia será utilizar simulação de Monte Carlo para estimar esse valor e comparar o método convencional com o método de variáveis antitéticas, mostrando seu efeito.

Exemplo

Suponha que queremos aproximar a integral \[ I = \int_0^1 e^u\,du \] por meio do Método de Integração de Monte Carlo.

De fato, \[ I = E[e^U], \qquad \text{com} \qquad U \sim \text{U}(0,1) \] e sabemos que o valor verdadeiro da integral é \[ I=e-1\approx1,71828. \]

Monte Carlo tradicional

Considere

\[ U_1,U_2,\ldots,U_n \overset{iid}{\sim}\text{U}(0,1). \]

O estimador tradicional de Monte Carlo é

\[ \hat{I}_{MC} = \frac{1}{n} \sum_{i=1}^n e^{U_i}. \]

Pela Lei dos Grandes Números,

\[ \hat{I}_{MC} \longrightarrow E[e^U] = I \qquad \text{quando} \qquad n\rightarrow\infty. \] Além disso,

\[ E[\hat{I}_{MC}] = E\left[\frac{1}{n} \sum_{i=1}^n e^{U_i}\right] = \frac{1}{n}\sum_{i=1}^nE\left[ e^{U_i}\right] = \frac{1}{n}\sum_{i=1}^nE\left[ e^{U}\right] = \frac{1}{n}nE\left[ e^{U}\right] =E[I]=I, \]

portanto o estimador é não-viesado, e sua variância é

\[ \operatorname{Var}(\hat{I}_{MC}) = \operatorname{Var}\left(\frac{1}{n} \sum_{i=1}^n e^{U_i}\right) = \frac{1}{n^2}\left( \sum_{i=1}^n \operatorname{Var}(e^{U_i})\right) = \frac{1}{n^2}\left( n \operatorname{Var}(e^{U})\right) = \dfrac{\operatorname{Var}(e^U)}{n}. \]

Assim, quanto maior o número de simulações, menor será a variância do estimador.

Entretanto, o objetivo das técnicas de redução de variância é obter maior precisão sem simplesmente aumentar o número de simulações.

Variáveis antitéticas

A técnica de variáveis antitéticas explora a seguinte propriedade:

Se \(U\sim \text{U}(0,1)\), então \(1-U\sim \text{U}(0,1)\).

Portanto, em vez de gerar dois números aleatórios independentes, podemos gerar apenas um número \(U\) e utilizar o par

\[ U,\qquad 1-U. \]

Para cada \(U_i\), definimos \[ Y_i= \frac{e^{U_i}+e^{1-U_i}}{2}. \]

O estimador antitético será

\[ \hat{I}_A = \frac{1}{n} \sum_{i=1}^n \frac{e^{U_i}+e^{1-U_i}}{2}. \]

De fato, \[ E[\hat{I}_{A}] = E\left[ \frac{1}{n} \sum_{i=1}^n \frac{e^{U_i}+e^{1-U_i}}{2} \right ] = \frac{1}{n} \sum_{i=1}^n \frac{E[e^{U_i}]+E[e^{1-U_i}]}{2} = \frac{1}{n} \sum_{i=1}^n \frac{E[e^{U_i}]+E[e^{U_i}]}{2} = \frac{1}{n} \sum_{i=1}^n E[e^{U_i}] =I, \] portanto, o estimador \(\hat{I}_A\) também é não viesado.

Agora vamos determinar a variância de \(\hat{I}_{A}\).

Lembre-se que, para duas v.a.’s \(X\) e \(Y\), \[ \operatorname{Var} \left( \frac{X+Y}{2} \right) = \frac{1}{4} \left[ \operatorname{Var}(X) + \operatorname{Var}(Y) + 2\operatorname{Cov}(X,Y) \right]. \]

Logo, \[\begin{align*} \operatorname{Var}(\hat{I}_{A}) = \operatorname{Var}\left(\frac{1}{n} \sum_{i=1}^n \frac{e^{U_i}+e^{1-U_i}}{2}\right) &= \frac{1}{n^2}\sum_{i=1}^n\operatorname{Var}\left( \frac{e^{U_i}+e^{1-U_i}}{2}\right) \\ &= \frac{1}{n^2}\sum_{i=1}^n\left[\frac{1}{4}\left(\operatorname{Var}( e^{U_i})+\operatorname{Var}(e^{1-U_i})+2\operatorname{Cov}(e^{U_i},e^{1-U_i})\right)\right]. \end{align*}\]

Como \(e^{U_i}\) e \(e^{1-U_i}\) possuem a mesma distribuição, \[ \operatorname{Var}(\hat{I}_{A}) = \frac{1}{n^2}\sum_{i=1}^n\left[\frac{1}{4}\left(2\operatorname{Var}( e^{U_i})+2\operatorname{Cov}(e^{U_i},e^{1-U_i})\right)\right] = \frac{1}{n^2}\sum_{i=1}^n\left[\frac{1}{2}\left(\operatorname{Var}( e^{U_i})+\operatorname{Cov}(e^{U_i},e^{1-U_i})\right)\right]. \]

Cálculo da variância de \(e^U\):

\[ \operatorname{Var}(e^U) = E[e^{2U}] - (E[e^U])^2 = \left[\int_0^1 e^{2u}\,du\right] - (e-1)^2 = \left[\left.\frac{e^{2u}}{2}\right|_0^1\right] - (e-1)^2 = \frac{e^2-1}{2} - (e-1)^2. \] Numericamente, \[ \operatorname{Var}(e^U) \approx0,242036 . \]

Cálculo da covariância antitética:

\[ \operatorname{Cov}(e^U,e^{1-U}) = E[e^Ue^{1-U}] - E[e^U]E[e^{1-U}] = E[e]-E[e^U]E[e^U] = e-(e-1)^2. \] Numericamente,

\[ \operatorname{Cov}(e^U,e^{1-U}) \approx-0,234211 . \]

Comparação dos dois estimadores

Suponha que utilizemos \(U_1,U_2\) independentes. O estimador baseado em duas observações é

\[ \hat{\theta}_{I} = \frac{e^{U_1}+e^{U_2}}{2}. \]

Como as observações são independentes,

\[ \operatorname{Cov}(e^{U_1},e^{U_2})=0. \]

Assim,

\[ \operatorname{Var}(\hat{\theta}_I) = \frac{\operatorname{Var}(e^U)}{2}= \frac{0,242036}{2} \approx0,121018. \]

Por outro lado, considerando as variáveis dependentes \(U\) e \(1-U\), de forma que

\[ \hat{\theta}_A = \frac{e^U+e^{1-U}}{2}. \]

Temos

\[ \operatorname{Var}(\hat{\theta}_A) = \frac{ \operatorname{Var}(e^U) + \operatorname{Cov}(e^U,e^{1-U}) }{2} = \frac{ 0,242036-0,234211 }{2}\approx0,003912. \]

o que significa que a técnica reduziu significativamente a variância do estimador.

Implementação

Vamos implementar o Método de Monte Carlo para estimar a integral do exemplo utilizando

o método tradicional; e
variável de controle.

set.seed(123)

# resultado exato
I = exp(1) - 1
I

## [1] 1.718282

M = 10000       # número de experimentos independentes
n = 1000        # número de pares

estimativas_ind = numeric(M)
estimativas_ant = numeric(M)

for (j in 1:M) {

  # --------------------------------------------------------
  # Método tradicional
  # --------------------------------------------------------
  u_ind = runif(2*n)
  estimativas_ind[j] = mean(exp(u_ind))

  # --------------------------------------------------------
  # Método antitético
  # --------------------------------------------------------
  u = runif(n)
  y_ant = (exp(u) + exp(1 - u)) / 2
  estimativas_ant[j] = mean(y_ant)
}

# ----------------------------------------------------------
# Resultados empíricos
# ----------------------------------------------------------

resultados_empiricos = data.frame(
  Medida = c( 
     "Valor verdadeiro", 
     "Média independente", 
     "Variância independente", 
     "Média antitética",
     "Variância antitética"
     ),
  Valor = c( 
       I, 
       mean(estimativas_ind), 
       var(estimativas_ind),
       mean(estimativas_ant), 
       var(estimativas_ant)
     )
   ) 

resultados_empiricos

##                   Medida        Valor
## 1       Valor verdadeiro 1.718282e+00
## 2     Média independente 1.718279e+00
## 3 Variância independente 1.231625e-04
## 4       Média antitética 1.718272e+00
## 5   Variância antitética 3.885850e-06

Generalizando a ideia

Suponha que estamos interessados em utilizar simulação de Monte Carlo para estimar a integral \(I = E[X]\).

Simulando duas v.a. \(X_1\) e \(X_2\) identicamente distribuídas, \[ \operatorname{Var} \left( \frac{X_1+X_2}{2} \right) = \frac{1}{4} \left[ \operatorname{Var}(X_1) + \operatorname{Var}(X_2) + 2\operatorname{Cov}(X_1,X_2) \right] = \frac{1}{2} \left[ \operatorname{Var}(X_1) + \operatorname{Cov}(X_1,X_2) \right]. \] Neste caso, se \(X_1\)e \(X_2\) forem negativamente correlacionadas, há a possibilidade de reduzir a variância de estimadores baseados em \(X_1\) e \(X_2\).

De fato, tomando \(U_1,...,U_m\) v.a. i.i.d. \(\text{U}(0,1)\), \[X_1 = h(U_1,...,U_m) \qquad \text{e} \qquad X_2 = h(1-U_1,...,1-U_m)\] serão v.a. identicamente distribuídas, uma vez que \(1-U_1,...,1-U_m\) também são i.i.d. \(\text{U}(0,1)\). Adicionalmente, se \(h\) é uma função monótona de cada um de seus argumentos, \(X_1\) e \(X_2\) também serão v.a. negativamente correlacionadas\(^{(\star)}\).

Utilizando este argumento, o método antitético propõe que, após a geração de \(X_1\) baseada na geração de \(U_1,...,U_m\), \(X_2\) seja baseado em \(1-U_1,...,1-U_m\) ao invés de \(m\) novas variáveis aletórias independentes. Esta estratégia poupa esforço computacional (pois evita a geração de um segundo conjunto de valores), além de reduzir a variância do estimador de Monte Carlo (ao menos quando \(h\) é uma função monótona).

\(^{(\star)}\)A demonstração deste resultado pode ser vista no Apêndice do Capítulo 9 de Ross (2022).

Redução de Variância por Variáveis de Controle

Suponha novamente que queremos utilizar simulação de Monte Carlo para estimar a integral \(I = E[X]\).

Para qualquer variável \(Y\) com média \(\mu_Y\) conhecida, e qualquer constante \(c\), temos que o estimador \[ \hat{I}_{C} = X + c(Y-\mu_Y) \] é não viesado e \[ \operatorname{Var}(\hat{I}_{C}) = \operatorname{Var}(\; X + c\,(Y-\mu_Y) \;) = \operatorname{Var}(\; X + c\,Y \;) = \operatorname{Var}(X) + c^2\,\operatorname{Var}(Y) + 2\,c\,\operatorname{Cov}(X,Y). \] É fácil provar que o valor de \(c\) que minimiza \(\operatorname{Var}(\hat{I})\) é \[ c^* = -\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(Y)}; \] e, neste caso, \[\begin{align*} \operatorname{Var}(\hat{I}_{C}) &= \operatorname{Var}(X) + \left(c^*\right)^2\,\operatorname{Var}(Y) + 2\,c^*\,\operatorname{Cov}(X,Y) \\ &= \operatorname{Var}(X) + \left(-\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(Y)}\right)^2\,\operatorname{Var}(Y) - 2\,\frac{\operatorname{Cov}(X,Y)}{\operatorname{Var}(Y)}\,\operatorname{Cov}(X,Y) \\ &= \operatorname{Var}(X) - \frac{\left(\operatorname{Cov}(X,Y)\right)^2}{\operatorname{Var}(Y)}. \end{align*}\]

De fato, \[ \frac{\operatorname{Var}(\hat{I}_{C})}{\operatorname{Var}(X)} = \frac{ \operatorname{Var}(\; X + c\,(Y-\mu_Y) \;)}{\operatorname{Var}(X)} = 1 - \left(\;Corr(X,Y)\;\right)^2; \] o que indica que o uso da variável de controle \(Y\) leva a uma redução de 100\(Corr(X,Y)\)% na variância do estimador (quando comparado ao estimador de Monte Carlo tradicional), desde que \(X\) e \(Y\) sejam correlacionadas.

Observação:

Na prática, se \(\operatorname{Cov}(X,Y)\) e \(\operatorname{Var}(Y)\) são desconhecidos, estima-se \(c^*\) utilizando os próprios valores simulados:

\[ \widehat{c^*} = -\frac{\widehat{Cov}(X,Y)}{\widehat{\operatorname{Var}(Y)}} = -\frac{\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})}{\sum_{i=1}^n(Y_i-\bar{Y})^2}. \]

Exemplo

Vamos voltar ao exemplo em que deseja-se aproximar a integral \[ I = \int_0^1 e^u\,du. \]

Assumindo a variável de controle \(Y=U \sim \text{U}(0,1)\), temos

\[\begin{align*} \operatorname{Var}(e^U) &= \frac{e^2-1}{2} - (e-1)^2 \approx 0,242036\\ \operatorname{Cov}(e^U,U) &= E[Ue^U] - E[e^U]E[U] = \int_0^1 ue^udu - \frac{1}{2}(e-1) = 1-\frac{e-1}{2} \approx 0,140859\\ \operatorname{Var}\left(\,e^U + c^*\,\left(u-\frac{1}{2}\right)\,\right) &= \operatorname{Var}(e^U) - \frac{\left(\, \operatorname{Cov}(e^U,U) \,\right)^2}{Var(U)} \approx 0,242036 - 12\,(0,140859)^2 = 0,003940; \end{align*}\]

portanto, \[ \frac{\operatorname{Var}(e^U)}{\operatorname{Var}\left(\,e^U + c^*\,\left(u-\frac{1}{2}\right)\,\right)} = \frac{0,003940}{0,242036} = 0,016279; \] ou seja, há uma redução de 98,37% na variância utilizando o estimador baseado na variável de controle \(U\).

Implementação

Vamos implementar o Método de Monte Carlo para estimar a integral do exemplo utilizando

o método tradicional;
variáveis antitéticas; e
variável de controle.

set.seed(123)

# resultado exato
I = exp(1) - 1
I

## [1] 1.718282

M = 10000       # número de experimentos independentes
n = 1000        # número de pontos

estimativas_ind = numeric(M)
estimativas_ant = numeric(M)
estimativas_con = numeric(M)

for (j in 1:M) {

  # --------------------------------------------------------
  # Método tradicional
  # --------------------------------------------------------
  u_ind = runif(n)
  estimativas_ind[j] = mean(exp(u_ind))

  # --------------------------------------------------------
  # Método antitético
  # --------------------------------------------------------
  u = runif(n)
  y_ant = (exp(u) + exp(1 - u)) / 2
  estimativas_ant[j] = mean(y_ant)
  
  # --------------------------------------------------------
  # Variável de controle
  # --------------------------------------------------------
  u = runif(n)
  x = exp(u)
  y = u
  
  x_con = x - (cov(x,y)/var(y))*(y-0.5)
  estimativas_con[j] = mean(x_con)  
}

# ----------------------------------------------------------
# Resultados empíricos
# ----------------------------------------------------------

resultados_empiricos = data.frame(
  Medida = c( 
     "Valor verdadeiro", 
     "Média independente", 
     "Variância independente", 
     "Média antitética",
     "Variância antitética",
     "Média controle",
     "Variância controle"     
     ),
  Valor = c( 
       I, 
       mean(estimativas_ind), 
       var(estimativas_ind),
       mean(estimativas_ant), 
       var(estimativas_ant),
       mean(estimativas_con), 
       var(estimativas_con)
     )
   ) 

resultados_empiricos

##                   Medida        Valor
## 1       Valor verdadeiro 1.718282e+00
## 2     Média independente 1.718366e+00
## 3 Variância independente 2.432425e-04
## 4       Média antitética 1.718302e+00
## 5   Variância antitética 3.872065e-06
## 6         Média controle 1.718218e+00
## 7     Variância controle 3.929707e-06

Referências

ROSS, Sheldon M. Simulation. 6. ed. London: Academic Press, 2022.

Integração via Método de Monte Carlo
Redução de Variância

Disciplina: GET00204 - Estatística Computacional 2

Redução de Variância por Variáveis Antitéticas

Exemplo

Monte Carlo tradicional

Variáveis antitéticas

Comparação dos dois estimadores

Implementação

Generalizando a ideia

Redução de Variância por Variáveis de Controle

Observação:

Exemplo

Implementação

Referências

Integração via Método de Monte CarloRedução de Variância

Disciplina: GET00204 - Estatística Computacional 2

Redução de Variância por Variáveis Antitéticas

Exemplo

Monte Carlo tradicional

Variáveis antitéticas

Comparação dos dois estimadores

Implementação

Generalizando a ideia

Redução de Variância por Variáveis de Controle

Observação:

Exemplo

Implementação

Referências

Integração via Método de Monte Carlo
Redução de Variância