Última atualização: 22 de maio, 2025 às 17:23.

Integração de Monte Carlo

De uma forma geral, a razão para simular pode ser formulada como a avaliação de 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]\) exato por por propriedade da distribuição gama a uma aproximação simulação.

5/2           ## valor teórico
## [1] 2.5
y <- rgamma(500, shape = 5, rate = 2)
mean(y)      ## aproximação por simulação
## [1] 2.466447
mean(rgamma(5000, shape = 5, rate = 2))  ## outra aproximação com N maior 
## [1] 2.513681
mean(rgamma(50000, shape = 5, rate = 2)) ## e mais outra
## [1] 2.502895

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). \]

Exemplo 2:

Seja \(Y \sim {\rm G}(5, 2)\) e desejamos obter \(P[Y > 4]\).
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.

pgamma(4, 5, 2, lower = FALSE)        ## valor teórico
## [1] 0.0996324
y <- rgamma(500, shape = 5, rate = 2) ## aproximação por siumlação
mean(y > 4)
## [1] 0.104

OBS: pode-se escrever a quantidade desejada da forma: \[P[Y > 4] = E[I_A(y)] = \int I_A(y) f(y) dy , \] que é portanto um valor esperado \(E[\phi(Y)] = \int \phi(y) f(y) dy\) em que \(\phi(y) = I_A(y) = \{y : y > 4\}\).

Em resumo, para obter (ou melhor, estimar) o valor esperado de \(\phi(y)\) por simulação seguimos os passos:

  1. obter uma amostra de \(y\),
  2. aplicar a função \(\phi(y)\) no valor obtido,
  3. repetir 1. e 2. muitas vezes,
  4. calcular a média dos valores obtidos no passo anterior.

Exemplo 3:

Seja \(X = (x_1, x_2) \sim N_2(\mu, \sigma, \rho)\)
com \(\mu = (\mu_1 = 0,5; \mu_2=0)^\prime\), \(\sigma = (\sigma^2_1 = 1^2; \sigma^2_1 = 1,2^2)^\prime\) e \(\rho = 0,5\).

Obter por simulação as quantidades a seguir e analiticamente, se possível.

  1. \(E[X_1+X_2]\)
  2. \(P[X_1 < 1 ; X_2 < 1]\)
  3. \(E[2X_1/(1-X_2)]\)
  4. \((1-P)^2\) em que \(P = P[X_1 < 1,5 ; X_2 > 1,2]\)

Iniciamos obtendo uma amostra da normal bivariada.
Usamos aqui a amostragem baseada na fatoração \[ f(x_1, x_2) = f(x_1) \cdot f(x_2|x_1).\] O algoritmo é então:

  1. simule \(x_1\) de \(f(x_1) \sim {\rm N}(\mu_1, \sigma^2_1)\),
  2. simule \(x_2\) de \(f(x_2|x_1) \sim {\rm N}(\mu_2 + \rho \frac{\sigma_2}{\sigma_1}(x_1 - \mu_1), (1-\rho^2) \sigma^2_2)\),
bvnsim <- function(n, mus = c(0,0), sigmas = c(1,1), rho = 0){
    x1 <- rnorm(n, mus[1], sigmas[1])
    x2 <- rnorm(n,
                mus[2] + (rho*sigmas[2]/sigmas[1]) * (x1-mus[1]),
                sqrt((1-rho^2) * sigmas[2]))
    cbind(x1, x2)
}

sam <- bvnsim(1000, mus=c(0.5, 0), sigmas = c(1, 1.2), rho=0.5)
dim(sam)

Uma discussão mais detalhada sobre amostragem da distribuição normal é disponível neste script.

Os valores simulados podem ser vistos no gráfico a seguir.

plot(sam, pch=19, cex=0.5)

O item a. pode ser verificado com o resultado analítico \(E[X_1 + X_2] = E[X_1] + E[X_2] = 0,5 + 0 = 0,5\). Por simulação teríamos.

mean(rowSums(sam))  
## [1] 0.5209159

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.

IND <- (sam[,1] < 1 & sam[,2] < 1)
mean(IND)  
## [1] 0.614
plot(sam, pch=19, cex=0.5, col=ifelse(IND, 2, 1))

No item c. \(E[\phi(X)] = E[2X_1/(1-X_2)] \approx N^{-1} \sum_i 2x_1i/(1-x_2i)\). Em R definimos uma função para calcular a quantidade desejada, aplicamos a cada par simulado e tiramos a média.

mean(apply(sam, 1, function(x) 2*x[1]/(1-x[2])))  
## [1] -1.571452

O item d. fica como exercício para o leitor.

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, pode-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 distribuição de interesse. A vantagem desse algoritmo é preservar as boas propriedades da distribuição normal multivariada.
Por exemplo, imagine que desejamos gerar dados da distribuição conjunta de uma variável aleatória com distribuição Poisson(3) e uma outra variável aleatória com distribuição Geométrica(1/5), de tal forma que ambas tenham em torno de 0.7 de correlação. Utilizando o algoritmo NORTA, geramos dados de uma distribuição normal multivariada fixando o parâmetro de correlação em 0.7 e convertemos os valores para as distribuições desejadas.

Algoritmo NORTA:

  1. Gere um vetor aleatório \(X \sim N(\boldsymbol{\mu}; \boldsymbol{\Sigma})\);
  2. Obtenha \(U = F_X(X)\) para cada marginal;
  3. 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 = (Y_1, Y_2)^T\), sendo \(Y_1 \sim N(0, 1)\), \(Y_2 \sim \text{Poisson}(2.5)\) e correlação \(\rho(Y_1, Y_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.
rho <- 0.8                      ## Correlação na normal bivariada.
X <- bvnsim(N, rho = rho) ## Gera matrix de pares (x1, x2)
dim(X)
## [1] 1000    2
head(X)
##              x1         x2
## [1,] -1.1955706 -0.2078853
## [2,] -0.3396138  0.4678280
## [3,]  1.4419015  0.8965611
## [4,] -0.3155906 -0.6801533
## [5,] -1.5737056 -0.8254713
## [6,] -0.2897338 -0.5566116
cor(X)
##           x1        x2
## x1 1.0000000 0.8068934
## x2 0.8068934 1.0000000

## Segundo passo: obter U = F_X(X). Qual é a distribuição de U? U ~ U[0,1]
u <- pnorm(X)

## Terceiro passo: Obter a distribuição de interesse das marginais via F^{-1}_(U).
Y1 <- X[, 1]                    ## Y_1 já é normal padrão.
Y2 <- qpois(u[, 2], 2.5)        ## qpois é a F^{-1} para Poisson.
Y <- cbind(Y1, Y2)
cor(Y)
##           Y1        Y2
## Y1 1.0000000 0.7837694
## Y2 0.7837694 1.0000000

Os gráficos das distribuições conjunta da das marginais é estão ne figura a seguir.

par(mfrow=c(1,3), mar=c(3,3,1,1), mgp = c(1.8,0.8,0))
plot(Y, pch = 19, cex = 0.3)
hist(Y1, freq = FALSE, nclass = 25, ylab = "f(y1)", main = "", col = "transparent"); lines(density(Y1))
curve(dnorm(x, 0, 1), add = TRUE, col = "blue", lwd = 2)
M2 <- max(Y2)
plot(0:M2, prop.table(table(factor(Y2, levels = 0:M2, labels = 0:M2))), type = "h", xlab = "Y2", ylab = "f(y2)")
lines((0:M2)+0.1, dpois(0:M2, lambda = 2.5), type = "h", col = "blue")

Exercícios

  1. Seja \(Y \sim {\rm Exp}(\lambda = 1/5)\). Obtenha a aproximação por simulação das quantidades solicitadas, notando que as quantidades podem ser escritas como valores esperados. Compare com valores teóricos.
  1. O valor esperado de \(Y\).
  2. \(E[Y^2]\)
  3. A variância de \(Y\)
  4. \(P[Y < 3]\)
  5. \(E[Y_1 + Y_2 + \ldots + Y_8]\) em que os \(Y_i\)’s são independentes e cada um tem distribuição igual a de \(Y\).
  6. \(Var[Y_1 + Y_2 + \ldots + Y_8]\) em que os \(Y_i\)’s são independentes e cada um tem distribuição igual a de \(Y\).
  7. \(P[Y_1 + Y_2 > 6]\) em que os \(Y_i\)’s são independentes e cada um tem distribuição igual a de \(Y\).

  1. Deseja-se gerar dados da distribuição conjunta de uma variável aleatória com distribuição Poisson(3) e uma outra variável aleatória com distribuição Geométrica(1/5), de tal forma que ambas tenham em torno de 0.7 de correlação. Obtenha as simulações utilizando o algoritmo NORTA. Faça gráficos das simulações da conjunta e marginais e compare com os gráficos da distribuições desejadas.

  2. Crie uma função chamada NORTA para gerar dados de um vetor de tamanho \(p\), em que seja possível o usuário especificar cada uma das \(p\) distribuições de interesse e a correlação entre cada componente do vetor.
    Dica: utilize a função rmvnorm() do pacote mvtnorm() para facilitar o primeiro passo do algoritmo.

Licença Creative Commons 4.0

Este conteúdo está disponível por meio da Licença Creative Commons 4.0