Neste experimento, deseja-se calcular o valor da integral \(\int_{0}^{1} e^{-x^2} dx\) utilizando uma chuva de pontos aleatórios sobre o gráfico da função. Trata-se de um método probabilístico para calcular integrais. Para fins de comparação, pode-se calcular a integral acima analiticamente, com ajuda da tabela da distribuição normal padrão (a tabela não computa os valores exatos, então, a rigor, este não é o valor exato da integral), segue:

Resolvendo \(\int_{0}^{1} e^{-x^2} dx\) com ajuda na Normal padrão:

Faça \(y = \sqrt{2}x\)
Temos que \(-x^2 = \frac{-y^2}{2}\) e \(dx = \frac{dy}{\sqrt{2}}\)
Substituindo na integral e modificando os limites de integração, temos:

\[\int_{0}^{\sqrt{2}} e^{-\frac{y^2}{2}}\frac{1}{\sqrt{2}} dy\]

Agora, multiplicamos e dividimos a integral por um termo que não depende de \(y\),
a saber, \(\sqrt{2\pi}\):

\[\sqrt{2\pi} \int_{0}^{\sqrt{2}} \frac{1}{\sqrt{2\pi}}e^{-\frac{y^2}{2}}\frac{1}{\sqrt{2}} dy\]

Com o que obtemos
\[\sqrt{\pi} \int_{0}^{\sqrt{2}} \frac{1}{\sqrt{2\pi}} e^{-\frac{y^2}{2}} dy\]

Note que a integral corresponde à \(P(0 < Z < \sqrt{2})\), em que \(Z\) é uma variável aleatória com distribuição Normal(0,1). Daí, a integral se resume a

\[\sqrt{\pi} (\Phi(\sqrt{2}) - \Phi(0))\]

Em que \(\Phi(z)\) é a densidade acumulada em \(z\).

Com auxílio do R, encontramos o valor da integral:

sqrt(pi)*(pnorm(sqrt(2)) - pnorm(0))  
## [1] 0.7468241

A técnica para calcular uma área utilizando uma “chuva de pontos” consiste em

  1. Gerar uma grande quantidade de pontos uniformente distribuídos no intervalo de integração desejado, em \(x\) e \(y\) e
  2. Computar a proporção de pontos que caíram dentro da área que se deseja calcular.

Considerando a função \(f(x) = e^{-x^2}\) os seguintes experimentos serão realizados:

  1. Variar a quantidade de pontos utilizados na aproximação para os valores \(N = 100, 1000, 10000, 100000\) com \(1000\) replicações e representar num boxplot;
  2. Computar a diferença entre os valores encontrados pelo método e o valor da integral encontrado analiticamente, assim como computar o erro percentual (\(\frac{calculado - real}{real}\)) e representar graficamente.

Inicialmente, para fazer as replicações para cada \(N\), escrevemos a seguinte função, que retorna os resultados para cada replicação:

integral <- function(N,rep)

replicate(rep, 
{
a <- runif(N)
b <- runif(N)

funcao <- function(x){
  exp(-x^2)
}

vetor = b <= funcao(a)
area = sum(vetor)/N
}
)

Agora, para criar o boxplot, criamos um data frame com os resultados para os valores de N desejados:

dados <- data.frame(N100 = integral(100,1000), N1000 = integral(1000,1000),
                    N10000 = integral(10000,1000), N100000 = integral(100000,1000))
colnames(dados) = c("N100", "N1000", "N10000", "N100000")

Para criar o boxplot que permitirá a comparação, utilizamos os seguintes comandos:

suppressPackageStartupMessages(library(ggplot2))
library(tidyr)
library(dplyr, warn.conflicts =  FALSE)

p <- dados %>% gather("replication","results") %>% ggplot(aes(x = replication, y = results, fill = replication))

p <- p + geom_boxplot(alpha=0.5)
p

Agora, faremos as comparações envolvendo o valor da integral calculado analiticamente \((0.7468241)\) e os valores obtidos pelo método probabilístico. Inicialmente, Criamos um data frame para guardar a diferença absoluta:

attach(dados)
diff <- abs(dados - 0.7468241) 
colnames(diff) <- c("diff100", "diff1000", "diff10000", "diff100000")

d <- diff %>% gather("replication","difference") %>% ggplot(aes(x = replication, y = difference, fill = replication))

d <- d + geom_boxplot(alpha=0.5)
d

Repetimos o procedimento, agora considerando o erro percentual:

percent <- diff/0.7468241 
colnames(percent) <- c("perc100", "perc1000", "perc10000", "perc100000")

r <- percent %>% gather("replication","percentual") %>% ggplot(aes(x = replication, y = percentual, fill = replication))

r <- r + geom_boxplot(alpha=.5)
r

#  Variar N e fazer, dentro de cada N, 100 replicações. ok
#  Fazer um boxplot com ggplot2 e escrever no rmarkdown ok
#  valor | diferença(calculado - valor real) | error percentual(calc - Valor real)/valor real 
#  Fazer boxplots