Introdução

Desejamos estimar o valor de \(\pi\), para isso vamos considerar um quadrado de lado \(L\) com uma circunferência interna de raio \(r = \frac{L}{2}\). A Área do quadrado é dada por \(A_{q} = L^2\) enquanto a Área da circunferência é dada por: \[A_{c} = \pi.r^2 = \frac{\pi.L^2}{4}\] Portanto, a razão entre a área da circunferência e a do quadrado é dada por: \[\frac{A_{c}}{A_{q}} = \frac{\frac{\pi.L^2}{4}}{L^2} = \frac{\pi}{4}\] Logo, temos que \[\pi = \frac{4.A_{c}}{A_{q}}\]

Estimação usando for

#pacotes
library(tidyverse)

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


#variavel para acumular total de pontos dentro
totalDentro <- 0

#inicialização de semente
set.seed("123")

for (i in 1:N)
{
  u.x <- runif(n = 1, min = -0.5, max = 0.5)
  
  u.y <- runif(n = 1, min = -0.5, max = 0.5)
  
  dentro <- u.x^2 + u.y^2<= 0.25
  
  totalDentro <- totalDentro + dentro 
  
}

#apaga i índice pois não será mais necessário
rm(i)

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

Estimação usando vetor

N <- 1000

set.seed("123")

#Dataframe
base <- tibble(Order = 1:N,
               U.x = runif(N, -1, 1),
               U.y = runif(N, -1, 1),
               Condicao = U.x^2 + U.y^2 <= 1, 
               Fator = ifelse(Condicao, "sim", "não")
               )

base2 <- base %>% 
  mutate(Acumulado = cumsum(Condicao),
         PiHat = 4*Acumulado/Order)
                
#Gráfico de convergência
base2 %>% ggplot(aes(x=Order, y = PiHat)) + 
          geom_line(color= "red") +
          xlim(0,N) + ylim(0,4) +
          geom_hline(yintercept =pi, color = "blue") + #reta em pi
          ggtitle("Estimação de Pi", 
                  subtitle = "Método de Monte Carlo")

Diagrama representativo da simulação

#Gráfico da frequencia
get_circle_tibble <- function(raio= 1, start=0,...){tibble(Angle = seq(from = start, to = 2*pi, ...), 
                                                    x = cos(Angle)*raio,
                                                    y = sin(Angle)*raio)}

ggplot() + 
  geom_path(aes(x,y),  
            color = "black", 
            data = get_circle_tibble(raio = 1, start = 0, .01)) +
  geom_path(aes(x,y),  
            color = "green",
            data = get_circle_tibble(raio = sqrt(2), start = -pi/4, pi/2)) + 
  geom_point(aes(x = U.x, y = U.y, color = Fator), data = base2)

Conclusão

Mostramos dois códigos para estimar \(\pi\). Observando os resultados podemos verificar que a estimativa se aproxima cada vez mais do valor real a medida que o número de iterações aumenta. Portanto, o valor estimado de \(\pi\) converge para o valor real de \(\pi\), como esperado.