• Integração de Monte Carlo é um método estatístico baseado em amostragem aleatória;
• Em 1777, Comte de Buffon utilizou experimentos aleatórios para determinar as probabilidades envolvidas no experimento da Agulha de Buffon;
• Gosset utilizou amostragem aleatória para estudar a distribuição \(t\) de Student;

• Seja \(g(x)\) uma função;
• Deseja-se determinar \(\int_a^b g(x)dx < \infty\);
• Note que se \(X\) é uma v.a. com densidade \(f(x)\), então: \[E[g(X)] = \int_{-\infty}^\infty g(x) f(x) dx.\]
• Então, se \(X_1, \dots, X_n\) é uma amostra aleatória desta \(f(x)\), então um estimador não-viciado de \(E[g(X)]\) é a média amostral de \(g(X_1), \dots, g(X_n)\).\[\hat{\theta} = \overline{g(X)} = \frac{1}{n} \sum_{i=1}^n g(X_i)\]

• Deseja-se estimar \(\int_0^1 g(x) dx\);
• Se \(X_1, \dots, X_n\) é uma a.a. de Uniforme(0,1), então\[\hat{\theta} = \overline{g(X)} = \frac{1}{n} \sum_{i=1}^n g(X_i)\]converge em probabilidade para \(E[g(X)] = \theta\).
• Assim, o estimador simples de Monte Carlo para \(\int_0^1 g(x) dx\) é \[\hat{\theta}=\overline{g(X)}\].

Determine uma estimativa de Monte Carlo para \(\int_0^1 e^{-x} dx\).

Determine uma estimativa de Monte Carlo para \(\int_0^1 e^{-x} dx\).

n = 1e5
unifs = runif(n)
gx = exp(-unifs)
mean(gx)

## [1] 0.6325936

Compare com o valor exato:

mean(gx)

## [1] 0.6325936

(gxExact = 1-exp(-1))

## [1] 0.6321206

Compare também com uma aproximação:

integrate(function(x) exp(-x), 0, 1)

## 0.6321206 with absolute error < 7e-15

estimaGx = function(n){
  unifs = runif(n)
  gx = exp(-unifs)
  return(mean(gx))
}
erroEstimativa = function(n){
  gxExact = 1-exp(-1)
  abs(gxExact-estimaGx(n))
}

n = 10^(0:7)
erro = sapply(n, erroEstimativa)
plot(n, erro, xlab='Tamanho da Amostra', log='xy',
     ylab='Erro Absoluto', main=expression(e^{-x}))

• Para estimar \(\int_a^b g(x)dx\), faça uma mudança de variáveis de modo que o novo intervalo seja \((0,1)\).
• Qual é essa transformação linear?

• Para estimar \(\int_a^b g(x)dx\), faça uma mudança de variáveis de modo que o novo intervalo seja \((0,1)\).
• Qual é essa transformação linear?\[y = \frac{x-a}{b-a}\]\[dy = \frac{1}{b-a} dx\]

\[y = \frac{x-a}{b-a}\]\[dy = \frac{1}{b-a} dx\]

\[ \int_a^b g(x)dx = \int_0^1 g(y(b-a)+b)(b-a)dy \\ \]Alternativamente, utilize uma Uniforme(a,b):\[\int_a^b g(x)dx = (b-a) \int_a^b g(x) \frac{1}{b-a}dx\]

Determine uma estimativa de Monte Carlo para \(\int_1^2 e^{-x} dx\).

n = 1e5
unifs = runif(n, 1, 2)
gx = exp(-unifs)
mean(gx)

## [1] 0.2323687

exp(-1) - exp(-2)

## [1] 0.2325442

Determine uma estimativa de Monte Carlo para\[\int_{-\infty}^x \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt.\]

Note que o intervalo não-limitado impede o uso direto do algoritmo dado. Entretanto, use a simetria da distribuição normal e divida o problema em duas partes:

• \(x \geq 0\)
• \(x < 0\)

Assim o problema se resume a \(\int_0^x \frac{1}{\sqrt{2\pi}} e^{-\frac{t^2}{2}} dt\).

Assim o problema se resume a \(\int_0^x e^{-\frac{t^2}{2}} dt\).

Efetuando a transformação de variáveis, temos:\[ \begin{aligned} y & = \frac{t}{x} \\ dt & = x dy \\ \theta & = \int_0^1 x e^{-\frac{(xy)^2}{2}} dy \\ \hat{\theta} & = E_Y\left[ x e^{-\frac{(xY)^2}{2}} \right] \end{aligned} \]

\[\hat{\theta} = E_Y\left[ x e^{-\frac{(xY)^2}{2}} \right]\]

n = 1e6
unifs = runif(n)
xs = seq(0, 3, by=.5)
Fxs = function(x, unifs)
  1/2+mean(x*exp(-(x*unifs)^2/2))/sqrt(2*pi)
gx = sapply(xs, Fxs, unifs)
Phix = pnorm(xs)
out = rbind(gx, Phix)
colnames(out) = xs
round(out, 4)

##        0    0.5      1    1.5      2    2.5      3
## gx   0.5 0.6915 0.8414 0.9333 0.9774 0.9940 0.9989
## Phix 0.5 0.6915 0.8413 0.9332 0.9772 0.9938 0.9987

Encontre um estimador de Monte Carlo para \(F_X(x)\) para \(X\sim N(0,1)\).

\[ F_X(x) = P(X \leq x) = E \left[ I(X \leq x) \right] \]

Então, é suficiente gerar \(X_1, \dots, X_n\) segundo uma Normal(0,1) e determinar a proporção de observações pertencentes ao evento \(\{X \leq x\}\), para qualquer \(x\) pré-especificado.

\(F_X(x)\) para \(X\sim N(0,1)\):

n = 1e6
xs = seq(0, 3, by=.5)
z = rnorm(n)
gx = sapply(xs, function(x, z) mean(z <= x), z)
Phix = pnorm(xs)
out = rbind(gx, Phix)
colnames(out) = xs
round(out, 4)

##           0    0.5      1    1.5      2    2.5      3
## gx   0.5001 0.6906 0.8412 0.9328 0.9772 0.9939 0.9986
## Phix 0.5000 0.6915 0.8413 0.9332 0.9772 0.9938 0.9987

\[ \mbox{V}(\hat{\theta}) = \frac{1}{n^2} \sum \mbox{V}_f\left[ g(x_i) \right] = \frac{\sigma^2}{n} \]

Ao utilizarmos a distribuição empírica de \(X\), temos:\[ \mbox{V}_f\left[ g(X) \right] = \frac{1}{n} \sum_{i=1}^n \left[ g(x_i) - \overline{g(x)}\right]^2 \]

Assim:

\[ \begin{aligned} \mbox{se}(\hat{\theta}) & = \frac{\sigma}{\sqrt{n}} = \frac{1}{n} \sqrt{\sum_{i=1}^n \left[ g(x_i) - \overline{g(x)}\right]^2} \end{aligned} \]

Segundo o Teorema do Limite Central, temos: \[\frac{\hat{\theta_n} - E(\hat{\theta_n})}{\sqrt{V(\hat{\theta_n})}} \sim N(0, 1), \mbox{quando } n \rightarrow \infty\]

Determine o IC(95%) para o estimador de Monte Carlo de \(F_X(x)\), quando \(X \sim Normal(0,1)\) e \(x=1.75\), utilizando a esperança segundo uma distribuição Uniforme. Note que \(F_X(1.75)=\) 0.9599.

n = 1e6
x = 1.75
unifs = runif(n)
gx = 1/2+x*exp(-(x*unifs)^2/2)/sqrt(2*pi)
v = mean((gx-mean(gx))^2)/n
ic = round(mean(gx) + c(-1, 1) * 1.96*sqrt(v), 4)
ic

## [1] 0.9596 0.9603

Determine o IC(95%) para o estimador de Monte Carlo de \(F_X(x)\), quando \(X \sim Normal(0,1)\) e \(x=1.75\), utilizando a esperança de uma função indicadora. Note que \(F_X(1.75)=\) 0.9599.

n = 1e6
x = 1.75
z = rnorm(n)
gx = z <= x
v = mean((gx-mean(gx))^2)/n
ic = round(mean(gx) + c(-1, 1) * 1.96*sqrt(v), 4)
ic

## [1] 0.9595 0.9602

Vimos que \[\theta = \int_a^b g(x) dx\] pode ser vista como \[(b-a) E[g(X)], \mbox{ quando } X\sim U(a,b).\]

Assim, o algoritmo para estimar \(\theta\) consiste em:

• \(X_1, \dots, X_n \sim U(a,b);\)
• \(\overline{g(X)} = \frac{\sum_i g(x_i)}{n};\)
• \(\hat{\theta} = (b-a) \overline{g(X)}.\)

\[E(\hat{\theta}) = (b-a) \frac{\theta}{b-a} = \theta\]

\[V(\hat{\theta}) = \frac{(b-a)^2}{m} Var[g(X)]\]

\[TLC: \hat{\theta} \rightarrow N\left(\theta, \frac{(b-a)^2}{m} Var[g(X)] \right)\]

Se \(X\) tem densidade \(f(x)\), uma estimativa de \(F(x) = \int_{-\infty}^x f(t)dt\) é:

• \(X_1, \dots, X_n \sim F_X\);
• \(g(X_i) = I(X_i \leq x)\)
• \(\hat{F(x)} = \overline{g(X)} = \frac{1}{n} \sum_{i=1}^n I(X_i \leq x)\)

Observe que:

• \(g(X_i) \sim B(F(x))\);
• \(E(\hat{F(x)}) = F(x) = p\);
• \(V(\hat{F(x)}) = \frac{F(x)[1-F(x)]}{n} = \frac{p(1-p)}{n} \leq \frac{1}{4n}\)

Se temos 2 estimadores, \(\theta_1\) e \(\theta_2\), para \(\theta\), então \(\theta_1\) é mais eficiente que \(\theta_2\) se:

\[\frac{V(\hat{\theta}_1)}{V(\hat{\theta}_2)} \leq 1.\]

Note que:

• Se estas variâncias forem desconhecidas, utilizamos suas estimativas;
• A variância é reduzida com o aumento do tamanho amostral, \(n\), então eficiência computacional é também importante.

Se temos 2 estimadores, \(\theta_1\) e \(\theta_2\), e \(V(\theta_2) < V(\theta_1)\), então o percentual de redução de variância é:

\[ 100 \left( \frac{V(\hat{\theta}_1)-V(\hat{\theta}_2)}{V(\hat{\theta}_1)} \right) \%\]

Se \(U_1\) e \(U_2\) são independentes e identicamente distribuídas:

\[V\left( \frac{U_1+U_2}{2} \right) = \frac{V(U_1)+V(U_2)}{4}.\]

Se \(U_1\) e \(U_2\) forem negativamente correlacionadas, então:

\[V\left( \frac{U_1+U_2}{2} \right) = \frac{V(U_1)+V(U_2)+2Cov(U_1, U_2)}{4},\]

note que \(Cov(U_1,U_2) < 0\).

Se \(X_1, \dots, X_n\) são independentes e \(f\) e \(g\) são funções crescentes, então:

\[E\left[f(X)g(X)\right] \geq E\left[f(X)\right]E\left[g(X)\right].\]

Se \(g=g(X_1, \dots, X_n)\) é monótona, então:

\[Y=g(F_X^{-1}(U_1), \dots, F_X^{-1}(U_n))\] e \[W=g(F_X^{-1}(1-U_1), \dots, F_X^{-1}(1-U_n))\]

são negativamente correlacionadas.

Se \(m\) replicatas de Monte Carlo são necessárias, então:

1. Gere \(m/2\) replicatas utilizando \(Y_j=g(F_X^{-1}(U_1^{(j)}), \dots, F_X^{-1}(U_n)^{(j)})\) e
2. Gere \(m/2\) replicatas utilizando \(W_j=g(F_X^{-1}(1-U_1^{(j)}), \dots, F_X^{-1}(1-U_n)^{(j)})\),

com \(U_i^{(j)} \sim U(0,1)\). Assim, o estimador antitético é:

\[ \begin{aligned} \hat{\theta} & = \frac{Y_1 + W_1 + Y_2 + W_2 + \dots + Y_{m/2} + W_{m/2}}{m} \\ & = \frac{1}{m}\sum_{i=1}^{m/2} \left(Y_i+W_i\right) = \frac{2}{m}\sum_{i=1}^{m/2} \frac{Y_i+W_i}{2} \end{aligned} \]

Anteriormente, estimamos:

\[\Phi(x) = \int_{-\infty}^x \frac{1}{\sqrt{2\pi}}e^{-t^2/2}dt\]

utilizando:

\[\theta = E_U\left[xe^{-(xU)^2/2}\right],\]

com \(U\sim U(0,1)\).

Restringindo a simulação para a cauda superior, \(g(.)\) é monótona. Então, podemos utilizar o método de variáveis antitéticas.

• Gerar \(u_1, \dots, u_{m/2} \sim U(0,1)\);
• Gerar \(Y_i = x e^{-(u_i x)^2/2}\);
• Gerar \(W_i = x e^{-((1-u_i) x)^2/2}\);
• Retornar \(\hat{\theta} = \frac{2}{m} \sum_i^{m/2} \frac{Y_i + W_i}{2}\).

n = 1e6
unifs = runif(n)
xs = seq(0, 3, by=.5)
Fxs = function(x, unifs)
  1/2+mean(x*exp(-(x*unifs)^2/2))/sqrt(2*pi)
gx = sapply(xs, Fxs, unifs)
Phix = pnorm(xs)
out = rbind(gx, Phix)
colnames(out) = xs
round(out, 4)

##        0    0.5      1    1.5      2    2.5      3
## gx   0.5 0.6915 0.8414 0.9332 0.9773 0.9938 0.9987
## Phix 0.5 0.6915 0.8413 0.9332 0.9772 0.9938 0.9987

n = 1e6
unifs = runif(n/2)
xs = seq(0, 3, by=.5)
Gxs = function(x, unifs)
  1/2+ mean((x*exp(-(x*unifs)^2/2)+x*exp(-(x*(1-unifs))^2/2))/2)/sqrt(2*pi)
gx2 = sapply(xs, Gxs, unifs)
Phix = pnorm(xs)
out = rbind(gx2, Phix)
colnames(out) = xs
round(out, 4)

##        0    0.5      1    1.5      2    2.5      3
## gx2  0.5 0.6915 0.8414 0.9332 0.9773 0.9938 0.9986
## Phix 0.5 0.6915 0.8413 0.9332 0.9772 0.9938 0.9987

mc1 = replicate(1000, Fxs(1.96, runif(1000)))
mc2 = replicate(1000, Gxs(1.96, runif(1000)))
round(sd(mc1), 5)

## [1] 0.00703

round(sd(mc2), 5)

## [1] 0.00033

round(100*(var(mc1)-var(mc2))/var(mc1), 2)

## [1] 99.79

Nosso objetivo é estimar \(\theta = E[g(X)]\). Suponha que exista \(f\) para a qual você conheça \(\mu = E\left[f(X)\right]\) e que ela seja correlacionado com \(g(X)\). Então, um estimador não-viciado de \(\theta\) é:

\[\hat{\theta}_c = g(X) + c(f(X)-\mu).\]

\[V(\hat{\theta}_c) = V(g(X)) + c^2 V(f(X)) + 2c Cov(g(X), f(X))\]

\(V(\hat{\theta}_c)\) é quadrática em \(c\). Portanto, seu valor de mínimo é atingido quando:

\[c^{*} = -\frac{Cov(g(X), f(X))}{V(f(X))}.\]

Utilizando o resultado anterior, então a variância mínima é:

\[ V(\hat{\theta}_c) = V(g(X)) - \frac{Cov^2(g(X), f(X))}{V(f(X))} \]

Assim, a redução de variância é:

\[ 100 \times \left\{ \frac{Cov^2(g(X), f(X))}{V(g(X))V(f(X))} \right\} = 100 \times Cor^2(g(X), f(X)). \]

• O método é melhor para casos de alta correlação entre \(f(X)\) e \(g(X)\);
• O ganho é nulo se estas funções forem não-correlacionadas;
• Para determinar \(c^{*}\) são necessárias \(Cov(f(X), g(X))\) e \(V(f(X))\);
• Se estas características forem desconhecidas, então podemos utilizar um experimento de Monte Carlo para estimá-las;

\[\theta = E\left[ e^U \right] = \int_0^1 e^u du,\]

expressão na qual \(U\sim U(0,1)\).

Se escolhemos \(f(U) = U\), então \(E(U)=1/2\) e \(V(U)=1/12\).

• \(Cov(f(U), g(U)) = ?\)
• \(c^{*} = -\frac{Cov(e^U, U)}{V(U)} = ?\)
• \(\hat{\theta}_c = g(U) + c^{*} (f(U)-\mu)\)

n = 1e6
U = runif(n)
T1 = exp(U)
ce = -12 + 6*(exp(1)-1)
T2 = exp(U) + ce*(U-0.5)
c(mean(T1), mean(T2), (var(T1)-var(T2))/var(T1), cor(exp(U), U)^2)

## [1] 1.718663 1.718287 0.983720 0.983720

\[\theta = \int_0^1 \frac{e^{-x}}{1+x^2} dx\]

Então, \(g(X) = \frac{e^{-x}}{1+x^2}\) e \(X \sim U(0, 1)\).

Escolhendo-se \(f(X) = \frac{e^{-0.5}}{1+x^2}\), temos \(E[f(X)] = \frac{\pi}{4\sqrt{e}}\).

• \(c^{*} = -\frac{Cov(f(X), g(X))}{V(f(X))} = ?\)
• \(\hat{\theta}_c = g(X) + c^{*} (f(X)-\mu)\)

• Em Cálculo, o valor médio de \(g(x)\) no intervalo \((a,b)\) é:

\[\frac{1}{b-a}\int_a^b g(x)dx;\]

• Isso é análogo à \(E[g(X)]\) quando \(X \sim U(a,b)\);
• Trata-se apenas da média de \(g(x)\) no referido intervalo utilizando a uniforme como função de ponderação;

• Assim, se \(X_1, \dots, X_n \sim U(a,b)\), então o estimador de Monte Carlo de \(\int_a^b g(x) dx\) é:

\[ \frac{b-a}{n}\sum_{i=1}^n g(X_i)\]

• Este método habitualmente falha em intervalos não-limitados;
• O método também pode ser ineficiênte se \(g(X)\) não for muito uniforme;

Ao considerarmos o problema de calcular uma integral como sendo uma determinação de esperanças, torna-se razoável considerar outras funções de ponderação (densidade). Esta estratégia nos leva ao Método de Amostragem por Importância.

• \(X \sim f(x)\);
• \(f(x) > 0, \forall \{x : g(x) > 0\}\);
• Seja \(Y\) a variável aleatória \(Y = g(X)/f(X)\), então:

\[E(Y) = \int \frac{g(x)}{f(x)} f(x) dx = \int g(x) dx\]

Assim, estimamos \(E(Y)\) utilizando:

\[\frac{1}{n} \sum_{i=1}^n \frac{g(x_i)}{f(x_i)}.\]

• A variância deste estimador é \(V(Y)/n\), portanto deve ser pequena;
• Quanto mais constante for \(Y\), menor será a variância;
• Desta maneira, devemos escolher uma \(f(x)\) bastante próxima de \(g(x)\).

Considere a integral

\[ \int_0^1 \frac{e^{-x}}{1+x^2} dx\]

e as seguintes funções candidatas:

\[ \begin{aligned} f_0(x) & = 1, & \forall & x \in (0,1); \\ f_1(x) & = e^{-x}, & \forall & x \in (0, \infty); \\ f_2(x) & = (1+x^2)^{-1}/\pi, & \forall & x \in (-\infty, \infty); \\ f_3(x) & = e^{-x}/(1-e^{-1}), & \forall & x \in (0, 1); \\ f_4(x) & = 4(1+x^2)^{-1}/\pi, & \forall & x \in (0, 1). \end{aligned} \]

• Todas as funções são positivas no intervalo \((0,1)\), onde \(g(x)>0\);
• As funções \(f_1\) e \(f_2\) tem intervalos maiores e vários dos valores simulados contarão como zero. Assim, elas são ineficientes;
• São densidades fáceis de serem simuladas;
• Exercício 1: fazer gráficos de \(g(x)\) e \(f_i(x)\);
• Exercício 2: estimar a integral usando amostragem por importância;
• Exercício 3: estimar o erro-padrão das estimativas acima;

\[ \theta = \int_A g(x) dx = \int_A \frac{g(x)}{\phi(x)} dx = E\left[ \frac{g(x)}{\phi(x)}\right]\]

Se \(X_1, \dots, X_n \sim F_X\), o estimador de MC é:

\[\frac{1}{n} \sum_{i=1}^n \frac{g(x_i)}{\phi(x_i)}\].

Assim, a variância:

\[ V(\hat{\theta}) = E[\hat{\theta}^2] - E^2[\theta] = \int_A \frac{g^2(x)}{\phi(x)} dx - \theta^2.\]

\[ V(\hat{\theta}) = E[\hat{\theta}^2] - E^2[\theta] = \int_A \frac{g^2(x)}{\phi(x)} dx - \theta^2\]

Podemos escolher a distribuição de \(X\) para minimizar a variância, que será obtida quando:

\[\phi(x) = \frac{|g(x)|}{\int_A |g(x)| dx}.\]

Regra geral: \(\phi(x) \approx |g(x)| f(x), \forall x \in A\).

Capítulo 5 - Maria Rizzo - Statistical Computing with R