Método

Montaremos um quadrado de lado \(l\), com área: \(A_{quad} = l^2\). Depois, circunscreveremos um círculo com \(r = l/2\) e área: \(A_{circ} = \pi r^2 = \frac{\pi l^2}{4}\).

Proporcionalmente, teremos que:

\[\frac{A_{circ}}{A_{quad}} = \frac{\pi}{4} \Rightarrow \pi = 4*\frac{A_{circ}}{A_{quad}}\] Podemos ver que \(p = \frac{A_{circ}}{A_{quad}}\) pode ser visto como uma proporção da quantidade de pontos dentro de cada área.

Sendo assim, construiremos duas variáveis aleatórias i.i.d X e Y, que seguem a distribuição uniforme com parâmetros \((-l, l)\).

\[ X,Y \sim U_{(-l, l)}\] Então, o vetor aleatório está definido no quadrado. Cada realização de tal vetor é um ponto jogado uniformemente aleatório dentro do quadrado.

Construamos, agora, duas variáveis aleatórias, tais que:

\(C_{n}\) = se o ponto n, ou n-ésima observação do vetor (X,Y), está dentro do círrculo.

\[C_{n} = \begin{cases} 1, & \text{se \(X^2 + Y^2 \leq (l/2)^2 \)} \\ 0, & \text{c.c.} \end{cases}\]

\(Q_{n}\) = se o ponto n, ou n-ésima observação do vetor \((X,Y)\), está dentro do quadrado.

\(Q_{n}\) é análogo a \(C_{n}\) em sua definição matemática.

Temos, então, que \(S(C_{n})\) = quantidade de pontos dentro do círculo depois de n realizações do vetor \((X,Y)\).

Porém, assim sendo, tempos que \(p\) pode ser reescrito como:

\[p = \frac{A_{circ}}{A_{quad}} = \frac{S(C_{n})}{S(Q_{n})}\] Contudo, nosso vetor \((X,Y)\) sempre irá jogar pontos dentro do quadrado, ou seja, \(S(Q_{n}) = n\), pois \(P(Q_{n} = 1) = 1\) \(\forall{n}\).

Enfim,

\[p = \frac{S(C_{n})}{n}\]

Como a sequência \((C_{n})_{n > 0}\) é i.i.d., teremos que pela Lei dos grandes números de Kogomorov,

\[p = \frac{S(C_{n})}{n} \rightarrow E(C_{n}) \\ \text{quase certamente ao} \space n \rightarrow \propto\] Ou seja,

\[ \pi \rightarrow 4*E(C_{n})\\ \text{quase certamente}\].

Ao aumentarmos o número de pontos aleatórios jogados dentro do quadrado, \(n\), mais perto \(\pi\) chegará do verdadeiro valor. Ou seja,

\[ \pi \approx 4*\frac{S(C_{n})}{n}\] Especialmente, quando \(n \rightarrow \propto\).

Código por Loop

#Pacote
#library(tidyverse)

#N = numero de interacoes

pi_loop <- function(N = 1000)
{
  
  #S(C_n)
  S_cn <- 0
  
  l = 2 #lado do nosso quadrado
  
  #Rodando N vezes
  for (i in 1:N) 
    {
    #Se quisermos testar:
    #cat("\nIteracao: ", i)
    
    u.x <- runif(n = 1, min = -l/2, max = l/2) #gerando nossos pares aleatorios,
    u.y <- runif(n = 1, min = -l/2, max = l/2) #realizacoes dos vetor aleatorio,
                                               #(X,Y)
    
    C_n <-  u.x^2 + u.y^2 <= (l/2)^2 #Valor de C_n ; l/2 = 1
    
    S_cn <- S_cn + C_n #Adicione 1 caso verdade, montando o S_cn.
    }
  
  pihat <- 4*(S_cn/N)
  
  return(pihat)
}

Código Vetorial

#library(tidyverse)


pi_vec <- function(N = 1000)
{
  
  l = 2 #lado do quadrado, poderia ser qualquer l > 0.
  
  dt <- tibble(
    Index = 1:N,
    U.x = runif(N, -l/2, l/2), #gerando nossos pares aleatorios,
    U.y = runif(N, -l/2, l/2), #do vetor aleatorio (X,Y)
    S_qn = 1:N,
    Circ = U.x^2 + U.y^2 <= (l/2)^2 , #C_n, porem 1 = TRUE, 0 = FALSE.
    S_cn = cumsum(Circ), #somando os C_n.
    Pihat = 4*S_cn/S_qn #podemos ver Pihat convergindo para pi.
  )
  
  a <- dt$Pihat[N]
  
  res <- list(Pi = a, Dt = dt)
  return(res)
}

Gráficos

#Plotting os graficos
#library(tidyverse)

N = 1000 

l = 2

dt <- tibble( 
  Index = 1:N,
  U.x = runif(N, -l/2, l/2),
  U.y = runif(N, -l/2, l/2),
  S_qn = 1:N,
  Circ = U.x^2 + U.y^2 <= (l/2)^2 ,
  S_cn = cumsum(Circ),
  Pihat = 4*S_cn/S_qn
)

#Plotaremos o gráfico para vermos essa convergencia de Pihat.

#Utilizaremos o ggplot que funciona por adicionar camadas daquilo que se deseja.

dt %>% ggplot(aes(x = Index, y = Pihat)) + #camada basica
  geom_line(color = "blue") + #vermos onde esta o nosso Pihat
  xlim(0,N) + ylim(0, 4) + # a nossa formula pihat = 4*S_cn/S_qn, mostra os limites de y.
  geom_hline(yintercept = pi, color = "red") + #verdadeiro valor de pi
  ggtitle("Estimação de Pi", subtitle = "Método de Monte Carlo") 

get_circle_tibble <- function(raio = l/2, start = 0, ...){ #queremos os valores tabelados que formam um circulo
  tibble(
    Angle = seq(from = start, to = 2*pi, ...), #montaremos a sequencia que discretizara o circulo inteiro
    x = cos(Angle)*raio, #escrevemos x e y em radianos
    y = sin(Angle)*raio
  )
}

ggplot() + 
  
  #Criando o nosso circulo.
  geom_path(aes(x, y), 
            color = "black",
            data = get_circle_tibble(raio = l/2, 
                                     start = 0, .01)) + #discretizando com 0.01 de diferenca no angulo
  #Criando nosso quadrado
  geom_path(aes(x, y),
            color = "green", 
            data = get_circle_tibble(raio = sqrt(l), # o nosso quadrado tem l = 2
                                     start = -pi/4, pi/2)) + #so queremos os quantros pontos do quadrado
  geom_point(aes(x = U.x, y = U.y, color = Circ), data = dt) #nossas realizacoes de (X,Y)