Inferência Bayesiana - GET/IME/UFF

Métodos de Monte Carlo

1 Método de Monte Carlo Simples

Seja \(h(\theta)\) uma função qualquer. Suponha que quero estimar \[I = \int_{a}^{b}{h(\theta)} d\theta.\] Podemos reescrever essa integral da seguinte forma \[I = \int_{a}^{b}{\frac{h(\theta)}{p(\theta)}p(\theta) d\theta} = E\left[g(\theta)=\frac{h(\theta)}{p(\theta)}\right],\] sendo \(p(\theta)\) a função de densidade de probabilidade de uma variável aleatória \(\theta\).

O método de Monte Carlo simples estima isto da seguinte forma

• Gera-se uma amostra aleatória simples \((\theta_1, \ldots, \theta_n)\) da distribuição \(p(\theta)\);
• Calcula-se a função g para cada valor gerado, ou seja, \(g(\theta_1), \ldots, g(\theta_n)\);
• Calcula-se um valor aproximado para a integral desejada da seguinte forma: \[\hat{I} = \bar{g} = \frac{1}{n}\sum_{i=1}^n{ g(\theta_i) }\]

1.1 Convergência

O método de Monte Carlo simples é baseado na geração de amostras aleatórias de uma distribuição e no cálculo da média das avaliações da função nessas amostras. A precisão da estimativa aumenta à medida que o tamanho da amostra aumenta, de acordo com a Lei dos Grandes Números. Quanto maior o número de amostras, menor será a variação e maior será a precisão da estimativa.

No entanto, é importante ter em mente que o método de Monte Carlo simples não garante a convergência rápida ou a precisão para todas as funções ou tipos de integrais. Em alguns casos, pode ser necessário um número muito grande de amostras para obter uma estimativa precisa, especialmente quando a função tem um comportamento complexo ou está definida em um espaço de alta dimensão.

Além disso, o método de Monte Carlo simples pode ter uma alta variabilidade quando a função tiver picos ou regiões de baixa probabilidade. Isso ocorre porque as amostras são geradas aleatoriamente e podem não capturar adequadamente essas regiões. Nesses casos, podem ser necessárias técnicas avançadas, como amostragem estratificada, amostragem por rejeição ou métodos de redução de variância, para melhorar a precisão da estimativa.

Em resumo, o método de Monte Carlo simples pode ser uma técnica poderosa para estimar valores esperados de funções complexas, mas sua precisão pode variar dependendo de vários fatores. A escolha adequada do tamanho da amostra, técnicas de amostragem apropriadas e métodos de redução de variância podem ajudar a melhorar a precisão das estimativas.

• Como \(\theta_1, \theta_2, \ldots, \theta_n\) é uma amostra aleatória simples da distribuição \(p(\theta)\), i.e., \(\theta_i \stackrel{iid}{\sim} p(\theta)\), temos que \(g(\theta_i)\) é independente de \(g(\theta_j)\) quando \(i \neq j\) e que são v.a.s identicamente distribuídas.

• O estimador \(\hat{I} = \frac{1}{n} \sum_{i=1}^n{ g(\theta_i)}\) é um estimador não viesado para \(I = E[g(\theta)]\) pois \[\begin{eqnarray*} E[\hat{I}] &=& E\left[ \frac{1}{n} \sum_{i=1}^n{ g(\theta_i)}\right] = \frac{1}{n}\sum_{i=1}^n{ E[g(\theta_i)]} \stackrel{iid}{=} \frac{1}{n} n E[g(\theta)] = E[g(\theta)] = I. \end{eqnarray*}\]

• Pela Lei Forte dos Grandes Números, temos que \[\hat{I} \stackrel{q.c.}{\rightarrow} I\] quando \(n\) for suficientemente grande.

• A variância do estimador \(\hat{I}\) é \[\tau = Var[\hat{I}] \stackrel{indep}{=} \frac{1}{n^2} \sum_{i=1}^n{Var[g(\theta_i)]} \stackrel{ident. distr.}{=} \frac{1}{n} Var[g(\theta)].\] Note que um bom estimador para \(Var[g(\theta)]\) é a variância amostral \(\frac{\sum_{i=1}^n{(g(\theta_i) - \bar{g})^2}}{n-1}\) sendo \(\bar{g} = \hat{I}\). Logo, podemos estimar a variância \(\tau\) usando o seguinte estimador \[\hat{\tau} = \frac{1}{n(n-1)} \sum_{i=1}^n{(g(\theta_i) - \hat{I})^2}.\]

• Pelo Teorema Central do Limite, temos que, quando o tamanho da amostra \(n\) é suficientemente grande, \[\hat{I} \rightarrow W\] sendo \(W\) uma variável aleatória com distribuição \(N\left( I, \hat{\tau}\right)\).

• Podemos usar o resultado acima para testar a convergência: se tivermos uma amostra de \(\hat{I}\), a distribuição desta amostra tem que convergir para a distribuição normal quando o tamanho da amostra for suficientemente grande. Caso não convirja, pode ser que o tamanho da amostra não seja suficientemente grande. Também podemos estimar o erro deste estimador a partir da amostra: \(z_{\alpha/2} \sqrt{\hat{\tau}}\) sendo \(z_{\alpha/2}\) o quantil (\(1-\alpha/2\)) da distribuição normal padrão.

1.2 Exemplo 1

Suponha que quero estimar \(I = \int_{1}^3{\exp(-\theta)} d\theta\).

Posso reescrever isto da seguinte forma \[\begin{eqnarray*} I &=& \int_{1}^3{\exp(-\theta) (3-1) \frac{1}{3-1}} d\theta = E[2\exp(-\theta)] \end{eqnarray*}\] sendo \(\theta\) uma v.a. com distribuição Unif(1,3). Logo, pelo método de Monte Carlo simples, estimo \(I\) da seguinte forma:

• obtenho uma amostra de tamanho \(n\) da distribuição Unif(1,3).
• Para cada valor amostrado, calculo \(g(\theta_i) = 2\exp(-\theta_i), i=1, \ldots, n\).
• Depois calculo \(\hat{I} = \bar{g} = \frac{1}{n}\sum_{i=1}^n{ g(\theta_i)}\).

A integral acima é conhecida e assume o seguinte valor: \[I = -\exp(-3) + \exp(-1) = 0,3181.\] Logo, para \(n\) suficientemente grande, teremos \(\hat{I} \rightarrow 0,3181\).

n       = 1000 #tamanho da amostra
a       = 1 #definindo os parâmetros da distribuição da v.a.
b       = 3 #definindo os parâmetros da distribuição da v.a.
theta   = runif(n, a, b) #Gerando uma amostra da distribuição uniforme
gtheta  = 2*exp(-theta)
Ichapeu  = mean(gtheta)
#comparando o valor estimado com o verdadeiro
c(Ichapeu, -exp(-b)+exp(-a))

## [1] 0.3113276 0.3180924

Avaliando a convergência desta aproximação:

#definindo as variáveis das repetições do algoritmo
nite      = 1000
Ichapeu    = NULL

#executando o algoritmo
for(ite in 1:nite){
  theta   = runif(n, a, b) #Gerando uma amostra da distribuição uniforme
  gtheta  = 2*exp(-theta)
  Ichapeu[ite]  = mean(gtheta)
}

#calculando a média do estimador e comparando com o valor verdadeiro
c(mean(Ichapeu), -exp(-b)+exp(-a))

## [1] 0.3180816 0.3180924

#comparando a variancia do estimador com a estimativa tauchapeu
c(var(Ichapeu), sum((gtheta-mean(gtheta))^2)/(n*(n-1)) )

## [1] 3.150084e-05 3.268799e-05

#estimando o erro do estimador
alfa = 0.05
qnorm(alfa)*sqrt(var(Ichapeu))

## [1] -0.009231842

#Avaliando a convergência desta aproximação através de um histograma
hist(Ichapeu, freq=F, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     main="Histograma", ylab="densidade", xlab=expression(hat(I)))
curve(dnorm(x, mean(Ichapeu), sd(Ichapeu)), add=T, col="red", lwd=2)

#Testando normalidade
nivelsig = 0.05
T1 = shapiro.test(Ichapeu)
T2 = ks.test(Ichapeu, "pnorm", mean(Ichapeu), sd(Ichapeu))

c(T1$p.value, T2$p.value)

## [1] 0.01810871 0.40596555

if(T1$p.value<=nivelsig){print("Rejeito normalidade pelo teste de SW")}else{print("Não rejeito normalidade pelo teste de SW")}

## [1] "Rejeito normalidade pelo teste de SW"

if(T2$p.value<=nivelsig){print("Rejeito normalidade pelo teste de KS")}else{print("Não rejeito normalidade pelo teste de KS")}

## [1] "Não rejeito normalidade pelo teste de KS"

1.3 Exemplo 2

Suponha que quero estimar \(I = \int_{0}^1{\exp(-\theta)} d\theta\).

Posso reescrever isto da seguinte forma \[\begin{eqnarray*} I &=& \int_{0}^1{\frac{\exp(-\theta)}{p(\theta)} p(\theta)} d\theta = E\left[ \frac{\exp(-\theta)}{p(\theta)} \right] \end{eqnarray*}\] sendo \(p(\theta)\) a f.d.p. da distribuição Beta ou Uniforme, por exemplo.

Logo, pelo método de Monte Carlo simples, estimo \(I\) da seguinte forma:

• obtenho uma amostra de tamanho \(n\) da distribuição \(p(\theta)\).
• Para cada valor amostrado, calculo \(g(\theta_i) = \frac{\exp(-\theta_i)}{p(\theta_i)}\).
• Depois calculo \(\hat{I} = \bar{g} = \frac{1}{n}\sum_{i=1}^n{ g(\theta_i)}\).

Note que \[I = -\exp(-1) + \exp(0) = 0,6321.\] Logo, para \(n\) suficientemente grande, teremos \(\hat{I} \rightarrow 0,6321\).

1.3.1 Usando a distribuição uniforme

#Gerando uma amostra da distribuição Beta(1,1) = Unif(0,1)
n       = 10000 #tamanho da amostra
a       = 1
b       = 1
theta   = rbeta(n, a, b)
gtheta  = exp(-theta)/dbeta(theta, a, b)
Ichapeu  = mean(gtheta)
c(Ichapeu, -exp(-1) + exp(0))

## [1] 0.6321716 0.6321206

#avaliando a convergência

#definindo as variáveis das repetições do algoritmo
nite      = 2000
Ichapeu    = NULL

#executando o algoritmo
for(ite in 1:nite){
  theta   = rbeta(n, a, b)
  gtheta  = exp(-theta)/dbeta(theta, a, b)
  Ichapeu[ite]  = mean(gtheta)
}

#calculando a média do estimador e comparando com o valor verdadeiro
c(mean(Ichapeu), -exp(-1) + exp(0))

## [1] 0.6320965 0.6321206

#comparando a variancia do estimador com a estimativa tauchapeu
c(var(Ichapeu), sum((gtheta-mean(gtheta))^2)/(n*(n-1)) )

## [1] 3.188070e-06 3.302436e-06

#estimando o erro do estimador
alfa = 0.05
qnorm(alfa)*sqrt(var(Ichapeu))

## [1] -0.002936914

#Avaliando a convergência desta aproximação através de um histograma
hist(Ichapeu, freq=F, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     main="Histograma", ylab="densidade", xlab=expression(hat(I)))
curve(dnorm(x, mean(Ichapeu), sd(Ichapeu)), add=T, col="red", lwd=2)

#Testando normalidade
nivelsig = 0.05
T1 = shapiro.test(Ichapeu)
T2 = ks.test(Ichapeu, "pnorm", mean(Ichapeu), sd(Ichapeu))

c(T1$p.value, T2$p.value)

## [1] 0.1077854 0.6568932

if(T1$p.value<=nivelsig){print("Rejeito normalidade pelo teste de SW")}else{print("Não rejeito normalidade pelo teste de SW")}

## [1] "Não rejeito normalidade pelo teste de SW"

if(T2$p.value<=nivelsig){print("Rejeito normalidade pelo teste de KS")}else{print("Não rejeito normalidade pelo teste de KS")}

## [1] "Não rejeito normalidade pelo teste de KS"

1.3.2 Usando a distribuição Beta

#Gerando uma amostra da distribuição Beta(4,3)
n       = 10000 #tamanho da amostra
a       = 4
b       = 3
theta   = rbeta(n, a, b)
gtheta  = exp(-theta)/dbeta(theta, a, b)
Ichapeu  = mean(gtheta)
c(Ichapeu, -exp(-1) + exp(0))

## [1] 0.6092649 0.6321206

#avaliando a convergência

#definindo as variáveis das repetições do algoritmo
nite      = 2000
Ichapeu    = NULL

#executando o algoritmo
for(ite in 1:nite){
  theta   = rbeta(n, a, b)
  gtheta  = exp(-theta)/dbeta(theta, a, b)
  Ichapeu[ite]  = mean(gtheta)
}

#calculando a média do estimador e comparando com o valor verdadeiro
c(mean(Ichapeu), -exp(-1) + exp(0))

## [1] 0.6260131 0.6321206

#comparando a variancia do estimador com a estimativa tauchapeu
c(var(Ichapeu), sum((gtheta-mean(gtheta))^2)/(n*(n-1)) )

## [1] 0.0503875593 0.0002871065

#estimando o erro do estimador
alfa = 0.05
qnorm(alfa)*sqrt(var(Ichapeu))

## [1] -0.3692231

#Avaliando a convergência desta aproximação através de um histograma
hist(Ichapeu, freq=F, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     main="Histograma", ylab="densidade", xlab=expression(hat(I)))
curve(dnorm(x, mean(Ichapeu), sd(Ichapeu)), add=T, col="red", lwd=2)

#Testando normalidade
nivelsig = 0.05
T1 = shapiro.test(Ichapeu)
T2 = ks.test(Ichapeu, "pnorm", mean(Ichapeu), sd(Ichapeu))

c(T1$p.value, T2$p.value)

## [1] 6.620824e-70 0.000000e+00

if(T1$p.value<=nivelsig){print("Rejeito normalidade pelo teste de SW")}else{print("Não rejeito normalidade pelo teste de SW")}

## [1] "Rejeito normalidade pelo teste de SW"

if(T2$p.value<=nivelsig){print("Rejeito normalidade pelo teste de KS")}else{print("Não rejeito normalidade pelo teste de KS")}

## [1] "Rejeito normalidade pelo teste de KS"

1.4 Exemplo 3

Para entendermos o exemplo a ser mostrado, lembre que se \[ \left( \begin{array}{c} \theta_1 \\ \theta_2 \end{array} \right) \sim N_2\left( \left( \begin{array}{c} \mu_1 \\ \mu_2 \end{array} \right), \left( \begin{array}{c c} V_{1} & \rho \sqrt{V_1 V_2} \\ \rho \sqrt{V_1 V_2} & V_{2} \end{array} \right) \right),\] então sua f.d.p. é \[p(\theta_1, \theta_2) = \frac{1}{2\pi\sqrt{V_1 \; V_2 \; (1-\rho^2) }} \exp\left\{-\frac{1}{2} A\right\},\] sendo \(A = \frac{1}{1-\rho^2}\left[ \left(\frac{\theta_1-\mu_1}{\sqrt{V_1}}\right)^2 -2 \rho \left(\frac{\theta_1-\mu_1}{\sqrt{V_1}}\right) \left(\frac{\theta_2-\mu_2}{\sqrt{V_2}}\right) + \left(\frac{\theta_2-\mu_2}{\sqrt{V_2}}\right)^2\right]\).

Suponha que queremos estimar \[I = \int_{-\infty}^{+\infty}{\int_{-\infty}^{+\infty}{\theta_1 \theta_2 \frac{1}{\pi\sqrt{3}} \exp\left\{-\frac{2}{3}\left[\theta_1^2-\theta_1\theta_2+\theta_2^2\right]\right\} d\theta_1 d\theta_2}}.\]

Ou seja, \(I = E[\theta_1 \theta_2]\) sendo \(\boldsymbol{\theta} = (\theta_1, \theta_2)\) um vetor aleatório com distribuição normal bivariada com \(\mu_1=\mu_2=0\), \(V_1=V_2=1\) e correlação \(\rho=0,5\).

• Então podemos obter uma a.a.s. da distribuição \(N_2\left( \left( \begin{array}{c} 0 \\ 0 \end{array} \right), \left( \begin{array}{c c} 1 & 0,5 \\ 0,5 & 1 \end{array} \right) \right)\).
• Calculamos então \(g(\boldsymbol{\theta}_i) = \theta_{1i} \theta_{2i}\)
• E depois calculamos \(\hat{I} = \bar{g}\).

Note que \(cov(\theta_1, \theta_2) = 0,5 = E[\theta_1 \theta_2] - E[\theta_1]E[\theta_2] = E[\theta_1 \theta_2]\). Logo, \(\hat{I}\) tem que convergir para 0,5.

#Gerando uma amostra da distribuição normal bivariada
n       = 5000 #tamanho da amostra
K       = 2
mu      = rep(0,K)
Sigma   = matrix(c(1,0.5,0.5,1),2,2)
#Sigma
library(MASS)
theta   = mvrnorm(n, mu, Sigma)
#dim(theta)
gtheta  = theta[,1]*theta[,2]

Ichapeu  = mean(gtheta)
#comparando o valor estimado com o verdadeiro
c(Ichapeu, 0.5)

## [1] 0.4791058 0.5000000

#avaliando a convergência

#definindo as variáveis das repetições do algoritmo
nite      = 5000
Ichapeu    = NULL

#executando o algoritmo
for(ite in 1:nite){
  theta   = mvrnorm(n, mu, Sigma)
  gtheta  = theta[,1]*theta[,2]
  Ichapeu[ite]  = mean(gtheta)
}

#calculando a média do estimador e comparando com o valor verdadeiro
c(mean(Ichapeu), 0.5)

## [1] 0.500242 0.500000

#comparando a variancia do estimador com a estimativa tauchapeu
c(var(Ichapeu), sum((gtheta-mean(gtheta))^2)/(n*(n-1)) )

## [1] 0.0002446572 0.0002343796

#estimando o erro do estimador
alfa = 0.05
qnorm(alfa)*sqrt(var(Ichapeu))

## [1] -0.02572801

#Avaliando a convergência desta aproximação através de um histograma
hist(Ichapeu, freq=F, bty="n", lwd=2, cex.lab=1.4, cex.axis=1.4,
     main="Histograma", ylab="densidade", xlab=expression(hat(I)))
curve(dnorm(x, mean(Ichapeu), sd(Ichapeu)), add=T, col="red", lwd=2)

#Testando normalidade
nivelsig = 0.05
T1 = shapiro.test(Ichapeu)
T2 = ks.test(Ichapeu, "pnorm", mean(Ichapeu), sd(Ichapeu))

c(T1$p.value, T2$p.value)

## [1] 0.1611005 0.8032278

if(T1$p.value<=nivelsig){print("Rejeito normalidade pelo teste de SW")}else{print("Não rejeito normalidade pelo teste de SW")}

## [1] "Não rejeito normalidade pelo teste de SW"

if(T2$p.value<=nivelsig){print("Rejeito normalidade pelo teste de KS")}else{print("Não rejeito normalidade pelo teste de KS")}

## [1] "Não rejeito normalidade pelo teste de KS"

1.4.1 Observação

• Note que se quisermos calcular a probabilidade de \(\theta_1 < 1\) usando o exemplo acima, podemos gerar \(n\) valores da distribuição normal bivariada e calcular quantos destes valores tiveram \(\theta_1 < 1\).

#Gerando uma amostra da distribuição normal bivariada
n       = 10000 #tamanho da amostra
K       = 2
mu      = rep(0,K)
Sigma   = matrix(c(1,0.5,0.5,1),2,2)
#Sigma
library(MASS)
theta   = mvrnorm(n, mu, Sigma)
length(which(theta[,1]<1))/length(theta[,1])

## [1] 0.8405

pnorm(1,0,1)

## [1] 0.8413447

• Note que se quisermos calcular \(Pr(\theta_1 < 1|\theta_2<1)\) usando o exemplo acima, podemos gerar \(n\) valores da distribuição normal bivariada, selecionar somente os pontos nos quais \(\theta_2<1\) e calcular quantos pontos selecionados tiveram \(\theta_1 < 1\).

#Gerando uma amostra da distribuição normal bivariada
n       = 10000 #tamanho da amostra
K       = 2
mu      = rep(0,K)
Sigma   = matrix(c(1,0.5,0.5,1),2,2)
#Sigma
library(MASS)
theta   = mvrnorm(n, mu, Sigma)

theta1_rest = theta[which(theta[,2]<1), 1]
length(which(theta1_rest < 1)) / length(theta1_rest)

## [1] 0.8889152

library(mvtnorm)
pmvnorm(mean = mu, sigma = Sigma, lower = -Inf, upper = c(1,1)) /pnorm(1, mu[2], sqrt(Sigma[2,2]))

## [1] 0.8857292
## attr(,"error")
## [1] 1e-15
## attr(,"msg")
## [1] "Normal Completion"

#P(theta1<1, theta2<1)/P(theta2<1)