Estimação de Monte Carlo em C

Introdução

Estamos interessados em estimar \(\pi\) utilizando o Método de Monte Carlo assim como fizemos na atividade 1, porém dessa vez iremos implementar em C. 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}}\]

Implementação

#include <Rcpp.h>

// [[Rcpp::export]]
double mc_pi(const int n){
  int i, r = 0;
  double a, b;
  Rcpp::RNGScope scope;

  for (i=0; i<n; i++){
    a = R::runif(0, 1);
    b = R::runif(0, 1);

    if (a*a + b*b <= 1)
      r++;
  }

  return (double) 4.*r/n;
}

Resultados

Tentaremos estimar pi utilizando vários tamanhos de n: n = 100, n = 1000 e n = 5000. Para cada valor de n faremos 100 execuções, de forma a se ter uma certa consistência nos resultados.

result1 <- c()

n = 100
replicas = 100
for(i in 1:replicas) {
  a <- mc_pi(n)
  result1 <- c(result1, a)
}

result2 <- c()

n = 1000
replicas = 100
for(i in 1:replicas) {
  a <- mc_pi(n)
  result2 <- c(result2, a)
}


result3 <- c()

n = 5000
replicas = 100
for(i in 1:replicas) {
  a <- mc_pi(n)
  result3 <- c(result3, a)
}

result1 <- as.data.frame(result1)
result1 <- cbind(result1, rep(100, replicas))
colnames(result1) <- c("est", "n")

result2 <- as.data.frame(result2)
result2 <- cbind(result2, rep(1000, replicas))
colnames(result2) <- c("est", "n")

result3 <- as.data.frame(result3)
result3 <- cbind(result3, rep(5000, replicas))
colnames(result3) <- c("est", "n")

datapi <- rbind(result1, result2, result3)

Observando o gráfico de boxplot a seguir, a linha amarela representa o valor de pi, notamos que quanto maior o n mais próximo de \(\pi\) as estimativas estão e com menor variação também, assim como esperado.

datapi %>% ggplot() + geom_boxplot(aes(x = factor(n), y = est, fill = factor(n))) + 
  geom_hline(yintercept =pi, color = "yellow") + xlab("n") +ylab("Estimativa")