A simulação de Monte Carlo é utilizada para estimar esperanças matemáticas 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.
Queremos estimar uma integral \(\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} \]
Se conseguirmos construir um conjunto de amostras correlacionadas negativamente com a amostra de \(X\)’s, podemos reduzir essa variância.
A ideia principal das variáveis antitéticas é gerar pares de amostras \(X_i\) e \(Y_i\) que sejam correlacionadas negativamente. Se tivermos uma segunda amostra \(Y_1, ..., Y_n\), não necessariamente independente de \(X\), com \[ E[Y] = E[X] = \theta \] e \[ V[Y] = V[X] = \sigma^2, \] então podemos definir um novo estimador: \[ \hat{\theta}_A = \frac{1}{2n} \sum_{i=1}^{n} (X_i + Y_i) \] cuja esperança é \[ \text{E}(\hat{\theta}_A)=\theta \] e 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) \] onde \(\rho = \text{Corr}(X, Y)\). Se \(\rho < 0\), então \(\text{Var}(\hat{\theta}_A) < \text{Var}(\hat{\theta})\), tornando o estimador mais eficiente.
Uma maneira 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.
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 <- (theta.hatx + theta.haty) / 2
# 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 = (theta_hatx + theta_haty) / 2
# 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.686052168748484)]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.1439204126724436e-05), np.float64(1.1832216036912732e-05), np.float64(1.978147697118848e-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 uma variável \(Y\) para a qual sabemos simular \(\mathbb{E}[X \mid Y]\). Podemos usar a seguinte propriedade da expectativa condicional:
\[ \mathbb{E}[\mathbb{E}[X \mid Y]] = \mathbb{E}[X] = \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{E}[X \mid Y]\) pode ser usado como estimador para \(\theta\), e ele tem uma variância menor ou igual à de \(X\).
Na prática, o algoritmo funciona da seguinte forma:
Voltando à estimativa de \(\pi\), considere o seguinte pseudo-algoritmo:
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. \]
Sabemos que \(\mathbb{E}[Z_i] = \pi / 4\). Vamos agora usar \(\mathbb{E}[Z_i \mid X]\) para reduzir a variância.
A probabilidade condicional é dada por:
\[ \mathbb{E}[Z_i \mid X = x] = \mathbb{P}(X^2 + Y^2 \leq 1 \mid X = x). \]
Como \(Y^2 \leq 1 - x^2\), obtemos:
\[ \mathbb{E}[Z_i \mid X = x] = \int_{-\sqrt{1-x^2}}^{\sqrt{1-x^2}} \frac{1}{2} \, dy = \sqrt{1 - x^2}. \]
Logo:
\[ \mathbb{E}[\mathbb{E}[Z_i \mid X]] = \mathbb{E}[\sqrt{1 - X^2}] = \pi / 4. \]
set.seed(42) # Para reprodutibilidade
# Função para estimar pi 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)
return(4 * mean(Z))
}
# Função para estimar pi com condicionamento
estimate_pi_conditioning <- function(n) {
X <- runif(n, -1, 1) # Gera X uniformemente em [-1,1]
g_Y <- sqrt(1 - X^2) # Esperança condicional E[Z | X]
return(4 * mean(g_Y))
}
# Número de simulações
n <- 100000
# Estimativas
pi_simple <- estimate_pi_simple(n)
pi_conditioning <- estimate_pi_conditioning(n)
# Resultados
cat("Estimativa de pi sem condicionamento:", pi_simple, "\n")Estimativa de pi sem condicionamento: 3.1378 Estimativa de pi com condicionamento: 3.142653 import numpy as np
np.random.seed(42) # Para reprodutibilidade
# Função para estimar pi 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)
return 4 * np.mean(Z)
# Função para estimar pi com condicionamento
def estimate_pi_conditioning(n):
X = np.random.uniform(-1, 1, n) # Gera X uniformemente em [-1,1]
g_Y = np.sqrt(1 - X**2) # Esperança condicional E[Z | X]
return 4 * np.mean(g_Y)
# Número de simulações
n = 100000
# Estimativas
pi_simple = estimate_pi_simple(n)
pi_conditioning = estimate_pi_conditioning(n)
# Resultados
print("Estimativa de pi sem condicionamento:", pi_simple)Estimativa de pi sem condicionamento: 3.14412Estimativa de pi com condicionamento: 3.1401670733896316Seja \(Y \sim \text{Exp}(1)\) e suponha que \(X \mid Y = y \sim \mathcal{N}(y, 4)\). Desejamos estimar \(\mathbb{P}(X > 1)\).
Definimos a variável:
\[ Z_i = \begin{cases} 1, & \text{se } X_i > 1, \\ 0, & \text{caso contrário}. \end{cases} \]
Assim, temos \(\mathbb{E}[\bar{Z}] = \mathbb{P}(X > 1)\) e \(\mathbb{E}[\mathbb{E}[Z_i \mid Y]] = \mathbb{P}(X > 1)\).
A probabilidade condicional é:
\[ \mathbb{E}[Z_i \mid Y = y] = \mathbb{P}(X_i > 1 \mid Y = y). \]
Sabendo que \(X \mid Y = y \sim \mathcal{N}(y, 4)\), temos:
\[ \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). \]
set.seed(42) # Para reprodutibilidade
# Função para estimar P(X > 1) sem condicionamento
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)
return(mean(Z))
}
# Função para estimar P(X > 1) com condicionamento
estimate_prob_conditioning <- function(n) {
Y <- rexp(n, rate = 1) # Y ~ Exp(1)
g_Y <- 1 - pnorm((1 - Y) / 2) # E[Z | Y] usando a CDF da normal
return(mean(g_Y))
}
# Número de simulações
n <- 100000
# Estimativas
prob_simple <- estimate_prob_simple(n)
prob_conditioning <- estimate_prob_conditioning(n)
# Resultados
cat("Estimativa de P(X > 1) sem condicionamento:", prob_simple, "\n")Estimativa de P(X > 1) sem condicionamento: 0.4905 Estimativa de P(X > 1) com condicionamento: 0.4903664 import numpy as np
from scipy.stats import norm
np.random.seed(42) # Para reprodutibilidade
# Função para estimar P(X > 1) sem condicionamento
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)
return np.mean(Z)
# Função para estimar P(X > 1) com condicionamento
def estimate_prob_conditioning(n):
Y = np.random.exponential(1, n) # Y ~ Exp(1)
g_Y = 1 - norm.cdf((1 - Y) / 2) # E[Z | Y] usando a CDF da normal
return np.mean(g_Y)
# Número de simulações
n = 100000
# Estimativas
prob_simple = estimate_prob_simple(n)
prob_conditioning = estimate_prob_conditioning(n)
# Resultados
print("Estimativa de P(X > 1) sem condicionamento:", prob_simple)Estimativa de P(X > 1) sem condicionamento: 0.48977Estimativa de P(X > 1) com condicionamento: 0.4900730050240222Exercício 1. Implemente a técnica das variáveis antitéticas para estimar \(E[\sin(X)]\) 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)\).