Objetivo:

Comparar algoritmo implementado na última aula com sua versão em R (use n=1.000,100.000,1.000.000).

Função em R

f_rayleighR = function(x, sigma) {
  if (any(x < 0)) return (0)
  stopifnot(sigma > 0)
  return((x / sigma^2) * exp(-x^2 / (2*sigma^2)))
}

Função em RCPP

#include <RcppArmadillo.h>
#// [[Rcpp::depends(RcppArmadillo)]]

#// [[Rcpp::export]]
#double f_rayleigh(double x, double sigma) {
#  double y = 0; // y = f(x)
#  if (sigma <= 0.0) {// div by zero or not defined here
#    throw std::range_error("Inadmissible value");
#  }
#  // pow(x, a) // power // x^a
#  y = (x/pow(sigma, 2))*exp(-pow(x,2)/(2*pow(sigma,2)));
#  return(y);
#}

#// [[Rcpp::export]]
#double d_chi2(double x, double k) {
#  double y = 0; // y = f(x)
#  
#  // pow(x, a) // power // x^a
#  y = (1/(std::tgamma(k/2)*pow(2,k/2)))*pow(x,(k/2) - 1)*exp(-x/2);
#  return(y);
#}

#// Rcpp:Rcout << "Testo" << std::endl;
#// Entrada: m n?mero de intera??es
#// Entrada: sigma dispers?o
#// Sa?da: n?meros seguindo a distrui??o de Rayleigh (via MCMC)
#// [[Rcpp::export]]
#Rcpp::DataFrame esqueleto(int m, double sigma) {
#  // arma::vec zeros = arma::zeros<arma::vec>(m);
#  // arma::vec uns = arma::ones<arma::vec>(m);
#  // arma::vec inteiros = arma::randi<arma::vec>(m, arma::distr_param(0, 9));
#  // arma::vec uniformes = arma::randu<arma::vec>(m);
#  // arma::vec normais = arma::randn<arma::vec>(m);
  
#  arma::vec x = arma::zeros<arma::vec>(m);
#  x(0) = arma::chi2rnd(1.0);
#  int k = 0;
#  arma::vec u = arma::randu<arma::vec>(m);
  
  
#  for(int i = 1; i < m; i++){
#    x(i) = arma::chi2rnd(1.0);
#    double xt = x(i-1);
#    double y = arma::chi2rnd(xt);
#    double num = f_rayleigh(y, sigma)*d_chi2(xt,y);
#    double den = f_rayleigh(xt, sigma)*d_chi2(y,xt);
    
#    if (u(i) <= num/den)
#      x(i) = y;
#    else {
#      x(i) = xt;
#      k++; // y is rejected
#    }
#  }
  
#  return(Rcpp::DataFrame::create(Rcpp::Named("Rayleigh")=x,
#                                 Rcpp::Named("Uniforme")=u)
#  );
#}

Comparando os dois algoritmos utilizando o microbenchmarking.

library(microbenchmark)
library(ggplot2)


# n = 1.000
#m1 = microbenchmark(
#     f_rayleighR(runif(n = 1, min = 1, max = 10),1),
#     f_rayleigh(runif(n = 1, min = 1, max = 10),1),
#     times = 1000,
#     setup=set.seed(21) )

#m1
#ggplot2::autoplot(m1)


# n = 100.000
#m2 = microbenchmark(
#     f_rayleighR(runif(n = 1, min = 1, max = 10),1),
#     f_rayleigh(runif(n = 1, min = 1, max = 10),1),
#     times = 100000,
#     setup=set.seed(21) )

#m2
#ggplot2::autoplot(m2)


# n = 1.000.000
#m3 = microbenchmark(
#     f_rayleighR(runif(n = 1, min = 1, max = 10),1),
#     f_rayleigh(runif(n = 1, min = 1, max = 10),1),
#     times = 1000000,
#     setup=set.seed(21) )

#m3
#ggplot2::autoplot(m3)

Resultados

Para N = 1000

#m1:

#Unit: microseconds
#                                        expr     min     lq      mean    median #     uq     max   neval
#f_rayleighR(runif(n = 1, min = 1, max = 10), 1) 72.664 78.822 144.03858 103.864 138.348 24584.504 1000
#f_rayleigh(runif(n = 1, min = 1, max = 10), 1)  11.905 15.190  27.24689  19.706  27.095  3767.832 1000

Para N = 100.000

#m2:

#Unit: microseconds
#                                          expr   min    lq     mean    median   uq     max    neval
#f_rayleighR(runif(n = 1, min = 1, max = 10), 1) 30.380 39.411 96.04647 60.759 91.548 61236.50 1e+05
#f_rayleigh(runif(n = 1, min = 1, max = 10), 1)  4.516  7.390  17.05373 11.084 16.011 19696.75 1e+05

Para N = 1000.000

#m3:

#Unit: microseconds
#                                        expr     min     lq     mean    median    uq      max    neval
#f_rayleighR(runif(n = 1, min = 1, max = 10), 1) 30.379 58.296 147.98050 73.074 112.075 240572.21 1e+06
#f_rayleigh(runif(n = 1, min = 1, max = 10), 1)  4.516  10.674 26.63541  13.137  19.705  57955.14 1e+06

Conclusão:

Ao implementar o algoritmo acima, iremos Observar que a função na versão Rcpp é mais veloz.