Integração de Monte Carlo
Paulo Justiniano
Integração de Monte Carlo
De uma forma geral, a razão para simular pode ser formulada como uma integral da forma \[ E[Y] = \int y f(y) dy ,\] que por sua vez pode ser aproximada por uma amostra \(y_1, y_2, \ldots, y_n\) como \[ E[Y] = \int y f(y) dy \approx n^{-1} \sum_{i=1}^n y_i .\]
Exemplo 1:
Seja \(Y \sim {\rm G}(5, 2)\). Vamos obter \(E[Y]\) por simulação e pelo valor teórico, que neste caso é conhecido por propriedade da distribuição gama.
O raciocínio pode ser generalizado para uma função qualquer \(\phi(X)\). \[ E[\phi(Y)] = \int \phi(y) f(y) dy \approx n^{-1} \sum_{i=1}^n \phi(y_i). \] Ou seja, para obter (ou melhor, estimar) o valor esperado de \(\phi(y)\) por simulação seguimos os passos:
- obter uma amostra de \(y\),
- aplicar a função \(\phi(y)\) no valor obtido,
- tirar a média dos valores obtido no passo anterior.
Exemplo 2:
Seja \(Y \sim {\rm G}(5, 2)\) e desejamos obter \(P[Y > 5]\).
Este valor pode se obtido de forma exata (por exemplo usando pgamma() do R).
Aqui vamos estimar o valor desta probabilidade obtendo uma aproximação por simulação.
y <- rgamma(500, shape = 5, rate = 2)
mean(y > 5)
## [1] 0.032
pgamma(5, 5, 2, lower = FALSE)
## [1] 0.02925269OBS: pode-se escrever a quantidade esperada como: \[P[Y > 5] = P[I_A(y)] = \int I_A(y) f(y) dy, I_A = \{y : y > 5\}.\]
Exemplo 3:
Seja \((x, y) \sim N_2(\mu, \sigma, \rho)\)
com \(\mu = (\mu_x = 0,5; \mu_y=0)^\prime\), \(\sigma = (\sigma^2_x = 1^2; \sigma^2_x = 1,2^2)^\prime\) e \(\rho = 0,5\).
Obter por simulação as quantidades a seguir e analiticamente, se possível.
- \(E[X+Y]\)
- \(P[X < 1 ; Y < 1]\)
- \(E[2X/(1-Y)]\)
- \((1-P)^2\) em que \(P[X < 1,5 ; Y > 1,2]\)
Iniciamos obtendo uma amostra da normal bivariada.
Usamos aqui a amostragem baseada na fatoração \[ f(x, y) = f(x) \cdot f(y|x).\] O algorítmo é então:
- simule \(x\) de \(f(x) \sim {\rm N}(\mu_x, \sigma^2_x)\),
- simule \(y\) de \(f(y|x) \sim {\rm N}(\mu_y + \rho \frac{\sigma_y}{\sigma_x}(x - \mu_x), (1-\rho^2) * \sigma^2_y)\),
bvnsim <- function(n, mus, sigmas, rho){
x <- rnorm(n, mus[1], sigmas[1])
y <- rnorm(n,
mus[2] + (rho*sigmas[2]/sigmas[1]) * (x-mus[1]),
sqrt((1-rho^2) * sigmas[2]))
cbind(x, y)
}
sam <- bvnsim(1000, mus=c(0.5, 0), sigmas = c(1, 1.2), rho=0.5)
dim(sam)Uma discussão mais detalhada osbre amostragem da distribuição normal é disponível (neste script)[MVnormal.R]
As valores simulados podem ser vistos no gráfico a seguir.
O item a. pode ser verificado com o resultado analítico \(E[X + Y] = E[X] + E[Y] = 0,5 + 0 = 0,5\). Por simulação teríamos.
Gerando várias vezes a amostra pode-se ver que este resultado é muito variável e portanto 1000 amostras não estimam bem a probabilidade.
Para o item b. obtemos o resultado a seguir também ilustrado na próxima figura.
Algoritmo NORTA
No Exemplo 3 foi apresentado um algoritmo de geração de dados de uma distribuição normal multivariada.
A vantagem desta distribuição é que o algoritmo de geração e a especificação dos parâmetros (médias, variâncias e correlações) são simples.
Por outro lado, em muitas situações não estamos interessados em variáveis com distribuição normal. Podemos ter interesse em gerar dados de váriveis multivariadas discretas, positivas ou ainda categóricas. Se fosse utilizado o mesmo princípio do Exemplo 3, isso levaria a um algoritmo de simulação extremamente complicado e com alto custo computacional - e certamente haveria dificuldades na definição dos parâmetros.
Como alternativa, oode-se utilizar o algoritmo NORmal To Anything (**NORTA*) nessas situações. Nesse algoritmo, geramos dados da distribuição normal multivariada e convertemos as marginais para a(s) distribuição(ões) de interesse. A vantagem desse algoritmo é preservar as boas propriedades da distibuição normal multivariada.
Por exemplo, imagine que desejamos gerar dados da distribuição conjunta de uma variável aleatória Poisson(2) e uma variável aleatória Geométrica(1/2) de tal forma que ambas tenham em torno de 0.7 de correlação. Utilizando o algoritmo NORTA, bastaria gerar dados de uma distribuição normal multivariada fixando o parâmetro de correlação em 0.7 e converter para as distribuições desejadas.
Algoritmo NORTA:
Gere um vetor aleatório \(X \sim N(\boldsymbol{\mu}; \boldsymbol{\Sigma})\);
Obtenha \(U = F_X(X)\) para cada marginal;
Os valores da distribuição de interesse são obtidos via \(F^{-1}_Y(U)\).
Exemplo 4
Para mostrar o funcionamento do algoritmo NORTA vamos gerar dados de um vetor aleatório \(Y = (X_1, X_2)^T\), sendo \(X_1 \sim N(0, 1)\), \(X_2 \sim \text{Poisson}(2.5)\) e correlação \(\text{\rho}(X_1, X_2) \approx 0.8\).
## Primeiro passo: gerar um vetor aleatório de uma normal bivariada com
## correlação igual a 0.8 (e cada marginal sendo uma distribuição normal padrão).
n <- 1000 ## Tamanho da amostra.
mu <- c(0, 0) ## Médias marginais.
sigmas <- c(1, 1) ## Variâncias marginais.
rho <- 0.8 ## Correlação.
X <- bvnsim(n, mu, sigmas, rho) ## Gera o vetor aleatório.
cor(X)
## x y
## x 1.000000 0.810112
## y 0.810112 1.000000
## Segundo passo: obter U = F_X(X). Qual é a distribuição de U?
u <- pnorm(X)
## Terceiro passo: Obter a distribuição de interesse via F^{-1}_Y(U).
X1 <- X[, 1] ## X_1 já é normal padrão.
X2 <- qpois(u[, 2], 2.5) ## qpois = F^{-1} da Poisson.
Y <- cbind(X1, X2)
plot(Y)
Exercício
- Crie uma função chamada NORTA para gerar dados de um vetor de tamanho \(n\), em que seja possível o usuário especificar cada uma das \(n\) distribuições de interesse e a correlação entre cada componente do vetor.
Dica: utilize a funçãormvnorm()do pacotemvtnorm()para facilitar o primeiro passo do algoritmo.