Dentre os métodos conhecidos para se estimar áreas, o método \(\textbf{hit or miss}\) estima a área de formas sobre superfícies gerando pontos aleatórios e verificando a proporção de pontos que recaem sobre a superfície. Usaremos o método em questão para a estimação do valor de \(\pi\). Seja \(S\) uma superfície bi-dimensional cuja área queremos estimar. Considere \(Q\) um retângulo de dimensões conhecidas circunscrito à superfície \(S\). Gerando pontos aleatórios e independentes sobre a superfície do retângulo \(Q\), teremos que cada ponto pode atingir ou não a superfície \(S\) nela circunscrita. Temos então a variável aletória: \[X = \left \{ 1\ se \ o\ ponto\ atinge\ S;\\\ \ 0\ caso\ contrario\right.\] Que claramente segue uma distribuição de Bernoulli com parâmetro \(p\), e também intuitivamente: \[p = P(X=1) = \frac{A_{S}}{A_{Q}}\] Em que \(A_{S}\) representa a área da superfície em questão e \(A_{Q}\) a área do retângulo. Repetindo-se a geração de pontos de forma independente \(n\) vezes, teremos então a nova variável aleatória: \[Y = \sum_{i=1}^{n}X_{i}\] no qual Y tem distribuição binomial com parâmetros \(n\) e \(p\), isto é: \[Y\sim bin(n,p)\] E finalmente, para um valor de \(n\) suficientemente grande, teremos \(Y\) convergindo em distribuição para uma variável \(Y{\underset{d}{\rightarrow}}Z\) tal que: \(Z\sim N(np,np(1-p))\). E portanto, a proporção de pontos que atingem a superfíce \(S\) terá distribuição assintótica: \[\frac{\sum_{i=1}^{n}X_{i}}{n} = \frac{Y}{n} \sim N(p,\frac{p(1-p)}{n})\] Logo, temos que a proporção \(\frac{\sum_{i=1}^{n}X_{i}}{n}\) é um estimador não-viesado e consistente para p, isto é: \[\widehat{p} = \frac{\sum_{i=1}^{n}X_{i}}{n}\\\] Assim: \(\widehat{\frac{A_{S}}{A_{Q}}} = \frac{\sum_{i=1}^{n}X_{i}}{n}\) e portanto: \(\widehat{A_{S}} = A_{Q}\frac{\sum_{i = 1}^{n}X_{i}}{n}\) (pois \(A_{Q}\) é um valor conhecido.) Usando os fatos acima, podemos estimar o valor \(\pi\) ao calcular a área de uma circuferência de raio \(1\). Como retangulo auxiliar, usaremos um quadrado circunscrito de lado igual a \(2\). De antemão sabemos que: \[\frac{A_{S}}{A_{Q}} = \frac{\pi}{4} \]

Implementações em R

Utilizando Loops

Com \(10000\) iterações

#número de iterações
N <- 10000

#variavel para acumular total de pontos dentro
totalDentro <- 0
set.seed(123)#condição inicial
for (i in 1:N){
  u.x <- runif(1, min =-1 , max = 1)#distribuição uniforme no eixo x, vai gerar um numero aleatório entre 0 e 1
  u.y <- runif(1, min = -1, max = 1)
  dentro <- u.x^2 + u.y^2 <= 1#teste se o ponto pertence ou não ao círculo limitado pela circuferência
  totalDentro <- totalDentro + dentro # somando 1 a cada condição verdadeira
}

piMat <- 4*totalDentro/N
cat("\n Pi =", round(piMat, digits = 5))
## 
##  Pi = 3.1416

Agora com \(100000\) iterações

N <- 100000

totalDentro <- 0
set.seed(123)
for (i in 1:N){
  u.x <- runif(1, min =-1 , max = 1)
  u.y <- runif(1, min = -1, max = 1)
  dentro <- u.x^2 + u.y^2 <= 1
  totalDentro <- totalDentro + dentro
}

piMat <- 4*totalDentro/N
cat("\n Pi =", round(piMat, digits = 5))
## 
##  Pi = 3.14044

O método por loops é bastante custoso computacionalmente e têm uma convergência lenta. Uma alternativa é usar uma abordagem vetorial: ###Abordagem Vetorial

library(tidyverse)
## -- Attaching packages --------------------- tidyverse 1.2.1 --
## v ggplot2 3.1.0       v purrr   0.3.0  
## v tibble  2.0.1       v dplyr   0.8.0.1
## v tidyr   0.8.2       v stringr 1.4.0  
## v readr   1.3.1       v forcats 0.4.0
## -- Conflicts ------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
N <- 10000

base <- tibble( #ou data.frame, ou data_frame
  Order = 1:N,
  U.x = runif(N, -1, 1),
  U.y = runif(N, -1, 1),
  Condicao = U.x^2 + U.y^2 <= 1,
  Fator = if_else(Condicao, "sim", "nao")
)


base2 <- base %>% mutate(Acumulado = cumsum(Condicao), 
                         PiHat = 4*Acumulado/Order)
cat("\n Pi =", round(piMat, digits = 5))
## 
##  Pi = 3.14044

Que com \(10.000\) interações apresenta resultado semelhante ao método por loops com \(100.000\) iterações.