Metropolis–Hastings algorithm

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

Usa-se o algoritmo Metropolis-Hastings para gerar uma amostra da distribuição Rayleigh.
A seguinte função (em R) analisa a densidade da Rayleigh:

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

Agora podemos observar a seguir a função escrita 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)
  );
}

Realizando a comparação solicitada dos 2 algoritmos através do 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:

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

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

N = 1.000.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:

Analisando os microbenchmarks e os plots gerados, vemos que a função em sua versão em Rcpp é mais rápida, sendo assim, menos custosa.