A simulação de Monte Carlo é utilizada para estimar esperanças e calcular integrais numéricas. No entanto, a precisão do método pode ser limitada pela variância das amostras geradas. As técnicas de redução de variância visam minimizar essa variabilidade sem aumentar significativamente o número de amostras.
Assuma que queremos estimar uma integral que pode ser expressa como a esperança de uma variável aleatória, \(\theta = E[X]\). O estimador de Monte Carlo usual é: \[ \hat{\theta} = \frac{1}{n} \sum_{i=1}^{n} X_i \] onde \(X_1, ..., X_n\) são variáveis aleatórias iid com média \(\theta\). A variância desse estimador é dada por \[ \text{Var}(\hat{\theta}) = \frac{\text{Var}(X)}{n} \]
Variáveis antitéticas visam reduzir essa variância. Para tanto, a ideia principal é gerar pares de amostras \(X_i\) e \(Y_i\) que sejam correlacionadas negativamente. Concretamente, se tivermos uma segunda amostra \(Y_1, ..., Y_n\) i.i.d entre si, com \(E[Y] = E[X] = \theta\) e \(V[Y] = V[X] = \sigma^2\), com \(Corr(X_i,Y_j)=\rho I(i=j)\) (portanto não necessariamente independente da sequência \(X_i\)), então podemos definir um novo estimador: \[ \hat{\theta}_A = \frac{1}{2n} \sum_{i=1}^{n} (X_i + Y_i) \] cuja esperança é justamente \[ \text{E}(\hat{\theta}_A)=\theta.\] Além disso, sua variância é \[ \text{Var}(\hat{\theta}_A) = \frac{1}{4n^2} \sum_{i=1}^{n} \left( \text{Var}(X) + \text{Var}(Y) + 2\text{Cov}(X,Y) \right) \] ou, simplificando, \[ \text{Var}(\hat{\theta}_A) = \frac{\text{Var}(X)}{n} \left( \frac{1 + \rho}{2} \right). \] Assim, se \(\rho < 0\), então \(\text{Var}(\hat{\theta}_A) < \text{Var}(\hat{\theta})\), tornando este estimador mais eficiente que o estimador de Monte Carlo usual com uma amostra de tamanho \(2n\).
Uma maneira (mas não a única!) de gerar variáveis negativamente correlacionadas é a partir de distribuições uniformes. Para isso, vamos assumir que \(X \sim h(U)\), para algum \(h\) conhecido, em que \(U \sim U(0,1)\). Sabemos, pelo método da transformação inversa, que sempre é possível encontrar \(h\) com essa propriedade (embora nem sempre seja simples!).
A ideia chave é então usar que \(1-U\) também tem distribuição \(U(0,1)\). Assim \(h(1-U)\) tem a mesma distribuição de \(X\) e, em particular, a mesma média \(\theta\). Além disso, pode-se mostrar que se \(h\) é monótono, \(h(U)\) tem correlação negativa com \(h(1-U)\).
O método das variáveis antitéticas pode ser implementado da seguinte forma:
Se estamos estimando a média de uma função de \(X\), \(g(X)\), o passo 4. é trocado para
Queremos estimar a seguinte integral: \[ I = \int_0^\infty \log(1 + x^2) e^{-x}\, dx \] Sabemos que podemos reescrever essa integral como uma esperança matemática: \[ I = E[\log(1 + X^2)], \quad X \sim \textrm{Exp}(1) \]
O algoritmo é:
\[ Z_i = \frac{1}{2} \left( \log(1 + X_i^2) + \log(1 + Y_i^2) \right) \]
set.seed(42)
n <- 50000
# Geração das amostras
u <- runif(n, 0, 1)
x <- -log(u)
y <- -log(1 - u)
# Cálculo das estimativas
theta.hatx <- log(1 + x^2)
theta.haty <- log(1 + y^2)
theta.hatz <- c(theta.hatx, theta.haty)
# Média das estimativas
mean_values <- c(mean(theta.hatx), mean(theta.haty), mean(theta.hatz))
cat("Médias das estimativas:\n")Médias das estimativas:[1] 0.6906937 0.6919855 0.6913396import numpy as np
np.random.seed(42)
n = 50000
# Geração das amostras
u = np.random.uniform(0, 1, n)
x = -np.log(u)
y = -np.log(1 - u)
# Cálculo das estimativas
theta_hatx = np.log(1 + x**2)
theta_haty = np.log(1 + y**2)
theta_hatz = np.concatenate((theta_hatx, theta_haty))
# Média das estimativas
mean_values = [np.mean(theta_hatx), np.mean(theta_haty), np.mean(theta_hatz)]
print("Médias das estimativas:")Médias das estimativas:[np.float64(0.6885117805461132), np.float64(0.683592556950855), np.float64(0.6860521687484841)]A seguir, repetiremos esse processo várias vezes para avaliar as variâncias dos diferentes estimadores.
set.seed(42)
n <- 50000
B <- 1000 # Número de repetições
theta.hatx <- numeric(B)
theta.haty <- numeric(B)
theta.hatz <- numeric(B)
for (b in 1:B) {
u <- runif(n, 0, 1)
x <- -log(u)
y <- -log(1 - u)
theta.hatx[b] <- mean(log(1 + x^2))
theta.haty[b] <- mean(log(1 + y^2))
theta.hatz[b] <- mean((log(1 + x^2) + log(1 + y^2)) / 2)
}
# Cálculo das variâncias das estimativas
var_values <- c(var(theta.hatx), var(theta.haty), var(theta.hatz))
cat("Variâncias das estimativas:\n")Variâncias das estimativas:[1] 1.227314e-05 1.160393e-05 2.119955e-06import numpy as np
np.random.seed(42)
n = 50000
B = 1000 # Número de repetições
theta_hatx = np.zeros(B)
theta_haty = np.zeros(B)
theta_hatz = np.zeros(B)
for b in range(B):
u = np.random.uniform(0, 1, n)
x = -np.log(u)
y = -np.log(1 - u)
theta_hatx[b] = np.mean(np.log(1 + x**2))
theta_haty[b] = np.mean(np.log(1 + y**2))
theta_hatz[b] = np.mean((np.log(1 + x**2) + np.log(1 + y**2)) / 2)
# Cálculo das variâncias das estimativas
var_values = [np.var(theta_hatx), np.var(theta_haty), np.var(theta_hatz)]
print("Variâncias das estimativas:")Variâncias das estimativas:[np.float64(1.1439204126724435e-05), np.float64(1.1832216036912739e-05), np.float64(1.9781476971188576e-06)]Ao estimar o valor esperado de uma variável aleatória \(X\), muitas vezes queremos melhorar a precisão da estimativa reduzindo a variância associada. Para isso, podemos introduzir uma variável auxiliar \(Y\), conhecida como variável de controle, que possui um valor esperado conhecido, \(\mathbb{E}[Y] = \mu\), e uma relação com \(X\) que nos ajuda a diminuir a variância.
Dado que queremos estimar \(\theta = \mathbb{E}[X]\), definimos uma nova estimativa ajustada:
\[ Z = X + c(Y - \mu), \]
onde \(c\) é uma constante a ser escolhida para minimizar a variância de \(Z\) e \(\mathbb{E}[Y]=\mu\). Ao calcular a expectativa de \(Z\), temos:
\[ \mathbb{E}[Z] = \mathbb{E}[X + c(Y - \mu)] = \mathbb{E}[X] + c \cdot (\mathbb{E}[Y] - \mu) = \mathbb{E}[X] = \theta. \]
Assim, a nova estimativa é não viesada. Agora, calculamos a variância:
\[ \text{Var}[Z] = \text{Var}[X + c(Y - \mu)] = \text{Var}[X] + c^2 \text{Var}[Y] + 2c \cdot \text{Cov}(X, Y). \]
Para minimizar essa variância, derivamos em relação a \(c\) e encontramos o valor ótimo:
\[ c^* = -\frac{\text{Cov}(X, Y)}{\text{Var}(Y)}. \]
Substituindo \(c^*\) na expressão da variância, obtemos:
\[ \text{Var}[Z] = \text{Var}[X] - \frac{[\text{Cov}(X, Y)]^2}{\text{Var}(Y)}. \]
Essa fórmula mostra que o uso da variável de controle reduz a variância de forma significativa.
Usando a correlação \(\text{Corr}(X, Y) = \text{Cov}(X, Y)/\sqrt{\text{Var}(X)\text{Var}(Y)}\), podemos reescrever a variância como:
\[ \text{Var}[Z] = \text{Var}[X] \cdot \left(1 - [\text{Corr}(X, Y)]^2\right). \]
A redução percentual da variância é dada por:
\[ \frac{\text{Var}[Z]}{\text{Var}[X]} = 1 - [\text{Corr}(X, Y)]^2. \]
Quanto maior a correlação entre \(X\) e \(Y\), maior será a redução da variância.
Na prática, as quantidades \(\text{Cov}(X, Y)\) e \(\text{Var}(Y)\) podem ser estimadas a partir de uma amostra gerada via Monte Carlo. Suponha que temos \(n\) amostras de \(X\) e \(Y\): \(X_1, \dots, X_n\) e \(Y_1, \dots, Y_n\). Então:
\[ \widehat{\text{Cov}}(X, Y) = \frac{1}{n-1} \sum_{i=1}^n (X_i - \bar{X})(Y_i - \bar{Y}), \]
\[ \widehat{\text{Var}}(Y) = \frac{1}{n-1} \sum_{i=1}^n (Y_i - \bar{Y})^2. \]
Com esses valores, aproximamos \(c^*\) por:
\[ \widehat{c}^* = -\frac{\widehat{\text{Cov}}(X, Y)}{\widehat{\text{Var}}(Y)}. \]
O estimador final é dado por:
\[ Z = X + \widehat{c}^* (Y - \bar{Y}). \]
Considere \(X = e^U\), onde \(U \sim \text{Unif}(0, 1)\). Vamos usar \(Y = U\) como variável de controle. Sabemos que \(\mathbb{E}[U] = 0.5\) e \(\textrm{Var}(U)=1/12\). Calculamos:
\[ \text{Cov}(e^U, U) = \mathbb{E}[e^U U] - \mathbb{E}[U]\mathbb{E}[e^U] = \int_0^1 xe^xdx - \frac{e-1}{2}.. \]
Fazendo substituição na primeira integral com \(u=x, v=e^x,\), encontramos:
\[ \text{Cov}(e^U, U) = 0.14086. \]
A variância de \(e^U + c^*(U - 0.5)\) é:
\[ \text{Var}[e^U + c^*(U - 0.5)] = 0.2420 - 12(0.14086^2) = 0.0039. \]
A redução de variância é:
\[ 1 - \frac{0.0039}{0.2420} = 0.9838 \quad \text{(98.38%)}. \]
Isso demonstra como o uso de variáveis de controle pode reduzir drasticamente a variância na estimativa.
Queremos estimar \(\theta = \mathbb{E}[X]\). Considere que conhecemos uma variável \(Y\) da qual sabemos simular, e que sabemos calcular \(\mathbb{E}[X \mid Y]\).
Dada uma amostra \(Y_1,\ldots,Y_n\), a técnica do condicionamento estima \(\theta\) via \[\frac{1}{n} \sum_{i=1}^n g(Y_i),\] em que \(g(Y_i)=E[X_i|Y_i]\).
Este estimador funciona pois, da Teoria das Probabilidades,
\[
\mathbb{E}[\mathbb{E}[X \mid Y]] = \mathbb{E}[X] = \theta.
\] Assim, \[\mathbb{E}\left[\frac{1}{n} \sum_{i=1}^n g(Y_i)\right] =\mathbb{E}\left[g(Y_1)\right]=\mathbb{E}[\mathbb{E}[X \mid Y_1]] = \theta.\]
Além disso, usando a fórmula da variância condicional, temos: \[ \text{Var}(X) = \mathbb{E}[\text{Var}(X \mid Y)] + \text{Var}(\mathbb{E}[X \mid Y]). \] Como consequência: \[ \text{Var}(X) \geq \text{Var}(\mathbb{E}[X \mid Y]). \]
Isso implica que \[\mathbb{V}\left[\frac{1}{n} \sum_{i=1}^n g(Y_i)\right]\leq \mathbb{V}\left[\frac{1}{n} \sum_{i=1}^n X_i \right],\] ou sejam este estimador tem uma variância menor ou igual à do método de Monte Carlo sem condicionamento.
Considere o seguinte pseudo-algoritmo para estimar \(\pi\), visto no capítulo de Monte Carlo:
No passo \(i\), para \(1 \leq i \leq n\):
Calcule: \[ \hat{\theta}_n = \frac{1}{n} \sum_{i=1}^n Z_i. \]
A estimativa final para \(\pi\) será: \[ \hat{\pi}_n = 4\hat{\theta}_n. \]
Vamos agora usar \(\mathbb{E}[Z_i \mid X]\) para reduzir a variância. Esta esperança condicional é dada por \[ \mathbb{E}[Z_i \mid X = x] = \mathbb{P}(X^2 + Y^2 \leq 1 \mid X = x)= \int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \frac{1}{2} \, dy = \sqrt{1 - x^2}. \] Assim, o estimador de Monte Carlo com condicionamento em \(x\) é dado por \[ 4 \cdot\frac{1}{n} \sum_{i=1}^n \sqrt{1 - X_i^2}. \]
set.seed(42) # Para reprodutibilidade
# Estima pi e o erro padrão - sem condicionamento
estimate_pi_simple <- function(n) {
U1 <- runif(n, -1, 1)
U2 <- runif(n, -1, 1)
Z <- ifelse(U1^2 + U2^2 <= 1, 1, 0)
p_hat <- mean(Z)
pi_hat <- 4 * p_hat
se <- 4 * sqrt(p_hat * (1 - p_hat) / n) # EP por binomial
list(estimate = pi_hat, se = se)
}
# Estima pi e o erro padrão - com condicionamento em X
estimate_pi_conditioning <- function(n) {
X <- runif(n, -1, 1)
g <- sqrt(1 - X^2) # E[Z | X]
pi_hat <- 4 * mean(g)
se <- 4 * sd(g) / sqrt(n) # EP via amostra de g
list(estimate = pi_hat, se = se)
}
# Número de simulações
n <- 10000
# Estimativas
res_simple <- estimate_pi_simple(n)
res_cond <- estimate_pi_conditioning(n)
# Resultados
cat("Sem condicionamento: est =", res_simple$estimate, " | EP =", res_simple$se, "\n")Sem condicionamento: est = 3.1272 | EP = 0.01652096 Com condicionamento: est = 3.130729 | EP = 0.009053532 import numpy as np
np.random.seed(42) # Para reprodutibilidade
# Estima pi e o erro padrão - sem condicionamento
def estimate_pi_simple(n):
U1 = np.random.uniform(-1, 1, n)
U2 = np.random.uniform(-1, 1, n)
Z = (U1**2 + U2**2 <= 1).astype(int)
p_hat = Z.mean()
pi_hat = 4 * p_hat
se = 4 * np.sqrt(p_hat * (1 - p_hat) / n) # EP por binomial
return pi_hat, se
# Estima pi e o erro padrão - com condicionamento em X
def estimate_pi_conditioning(n):
X = np.random.uniform(-1, 1, n)
g = np.sqrt(1 - X**2) # E[Z | X]
pi_hat = 4 * g.mean()
se = 4 * g.std(ddof=1) / np.sqrt(n) # EP via amostra de g
return pi_hat, se
# Número de simulações
n = 10000
# Estimativas
pi_simple, se_simple = estimate_pi_simple(n)
pi_conditioning, se_conditioning = estimate_pi_conditioning(n)
# Resultados
print(f"Sem condicionamento: est = {pi_simple} | EP = {se_simple}")Sem condicionamento: est = 3.1348 | EP = 0.016468846225525333Com condicionamento: est = 3.1540094743355818 | EP = 0.008885501290884038Seja \(Y \sim \text{Exp}(1)\) e suponha que \(X \mid Y = y \sim \mathcal{N}(y, 4)\). Desejamos estimar \(\mathbb{P}(X > 1)\).
Para isso, definimos a variável: \[ Z_i = \begin{cases} 1, & \text{se } X_i > 1, \\ 0, & \text{caso contrário}. \end{cases} \]
Por construção, \(\mathbb{E}[\mathbb{E}[Z_i \mid Y]] = \mathbb{P}(X > 1)\). Assim, é natural fazer um estimador Monte Carlo com condicionamento em \(X\). Para isso, vamos calcular a probabilidade condicional \[ \mathbb{E}[Z_i \mid Y = y] = \mathbb{P}(X_i > 1 \mid Y = y). \] Como \(X \mid Y = y \sim \mathcal{N}(y, 4)\), \[ \mathbb{P}(X_i > 1 \mid Y = y) = \mathbb{P}\left(\frac{X_i - y}{2} > \frac{1 - y}{2} \mid Y = y\right). \] Usando a função de distribuição normal padrão \(\Phi\), obtemos \[ \mathbb{E}[Z_i \mid Y = y] = 1 - \Phi\left(\frac{1 - y}{2}\right). \] Assim, o estimador é dado por \[\frac{1}{n}1 - \Phi\left(\frac{1 - Y_i}{2}\right),\] em que \(Y_i \sim Exp(1).\)
set.seed(42) # Para reprodutibilidade
# Estima P(X > 1) e o EP - sem condicionamento (Bernoulli)
estimate_prob_simple <- function(n) {
Y <- rexp(n, rate = 1) # Y ~ Exp(1)
X <- rnorm(n, mean = Y, sd = 2) # X | Y ~ N(Y, 4)
Z <- ifelse(X > 1, 1, 0) # Indicador
p_hat <- mean(Z)
se <- sqrt(p_hat * (1 - p_hat) / n) # EP do mean(Z)
list(estimate = p_hat, se = se)
}
# Estima P(X > 1) e o EP - com condicionamento (g_i = E[Z|Y_i])
estimate_prob_conditioning <- function(n) {
Y <- rexp(n, rate = 1)
g <- 1 - pnorm((1 - Y) / 2) # E[Z | Y]
est <- mean(g)
se <- sd(g) / sqrt(n) # EP via amostra de g
list(estimate = est, se = se)
}
# Número de simulações
n <- 100000
# Estimativas
res_simple <- estimate_prob_simple(n)
res_cond <- estimate_prob_conditioning(n)
# Resultados
cat("Sem condicionamento: est =", res_simple$estimate, " | EP =", res_simple$se, "\n")Sem condicionamento: est = 0.4905 | EP = 0.001580853 Com condicionamento: est = 0.4903664 | EP = 0.0005183444 import numpy as np
from scipy.stats import norm
np.random.seed(42) # Para reprodutibilidade
# Estima P(X > 1) e o EP - sem condicionamento (Bernoulli)
def estimate_prob_simple(n):
Y = np.random.exponential(1, n) # Y ~ Exp(1)
X = np.random.normal(Y, 2, n) # X | Y ~ N(Y, 4)
Z = (X > 1).astype(int)
p_hat = Z.mean()
se = np.sqrt(p_hat * (1 - p_hat) / n) # EP do mean(Z)
return p_hat, se
# Estima P(X > 1) e o EP - com condicionamento (g_i = E[Z|Y_i])
def estimate_prob_conditioning(n):
Y = np.random.exponential(1, n)
g = 1 - norm.cdf((1 - Y) / 2) # E[Z | Y]
est = g.mean()
se = g.std(ddof=1) / np.sqrt(n) # EP via amostra de g
return est, se
# Número de simulações
n = 100000
# Estimativas
prob_simple, se_simple = estimate_prob_simple(n)
prob_conditioning, se_conditioning = estimate_prob_conditioning(n)
# Resultados
print(f"Sem condicionamento: est = {prob_simple} | EP = {se_simple}")Sem condicionamento: est = 0.48977 | EP = 0.0015808078539152062Com condicionamento: est = 0.4900730050240221 | EP = 0.0005167915603236426Exercício 1. Implemente a técnica das variáveis antitéticas para estimar \(E[X^2]\) onde \(X \sim U(0, \pi)\). Compare a variância do estimador com o estimador usual.
Exercício 2. Considere a integral:
\[ \theta = \int_{0}^{\pi/4} \int_{0}^{\pi/4} x^2 y^2 \sin(x + y) \log(x + y) \, dx \, dy. \]
Exercício 3. Seja \(Y \sim \text{Exp}(1)\) e \(X \mid Y = y \sim \mathcal{N}(y, 4)\). Desejamos estimar \(\mathbb{P}(X > 1)\).
Exercício 4. Seja \(X\mid B\sim N(\mu_B,1)\) com \(B\sim\text{Bern}(p)\) e \(\mu_1=2\), \(\mu_0=-1\). Estime \(\theta=\mathbb P(X>1)\).
Exercício 5 Seja \(U\sim \mathrm{Unif}(0,1)\) e \(\theta=\mathbb E[U^p]\) com \(p\in\{1,2,3,4\}\).
Exercício 6 Seja \(U\sim\mathrm{Unif}(0,1)\).
Exercício 7 Seja \(X=-\log U\) e \(Y=-\log(1-U)\) com \(U\sim\mathrm{Unif}(0,1)\).
Exercício 8 Seja \(X=Y+\varepsilon\) com \(Y\sim\mathrm{Unif}(0,1)\) e \(\varepsilon\sim N(0,\sigma^2)\) independentes.