Atividade: Refazer o algoritmo da seção 22.3 (nova edição) em Rcpp. A função deve retornar um dataframe com as variâncias para os casos de N 10, 100, 200, 1000. Tanto para o tradicional e com usando a técnica de redução de variância.
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat Ginv(arma::mat u) {
arma::mat X = u.transform( [](double val) { return (sqrt(2)*tan(val*atan(1/sqrt(2)))); } );
return(X);
}
// [[Rcpp::export]]
arma::mat psi1(arma::mat x){
arma::mat P = x.transform([](double val) {return (exp(-(pow(val,2))/2)*sqrt(2)*atan(1/sqrt(2))*(1+((pow(val,2))/2)));});
return(P);
}
// [[Rcpp::export]]
arma::mat psi2(arma::mat x){
arma::mat P2 = x.transform([](double val) {return ((1- (pow(val,2))/2)*sqrt(2)*atan(1/sqrt(2))*(1+((pow(val,2))/2)));});
return(P2);
}
// [[Rcpp::export]]
arma::vec var(int N=10000,int n =50){
arma::mat commonG = Ginv(arma::randu<arma::mat>(n,N));
arma::mat p1g = psi1(commonG);
arma::mat theta_hat = arma::mean(p1g,0);
arma::mat p2g = psi2(commonG);
arma::mat mu_hat = arma::mean(p2g,0);
arma::mat samplecov = arma::sum((p1g - arma::repelem(theta_hat,50,1)) % (p2g - ((double)5/(double)6)),0)/n;
arma::mat samplevar = arma::sum(arma::square(p2g - ((double)5/(double)6)),0)/n;
arma::mat alphastar = samplecov/samplevar;
arma::mat theta_hat_c = theta_hat - (alphastar % (mu_hat - ((double)5/(double)6)));
arma::vec aux = arma::var(theta_hat,0,1);
double var1 = aux(0);
arma::vec aux2 = arma::var(theta_hat_c,0,1);
double varc = aux2(0);
double reduction = (100*(var1 - varc)/var1);
arma::vec variances = {var1,varc,reduction};
return(variances);
}
// [[Rcpp::export]]
Rcpp::DataFrame variancias(){
arma::Col<int> Nsize = {10,100,200,1000};
arma::vec size1 = var(Nsize(0),50);
arma::vec size2 = var(Nsize(1),50);
arma::vec size3 = var(Nsize(2),50);
arma::vec size4 = var(Nsize(3),50);
arma::vec v1 = {size1(0),size2(0),size3(0),size4(0)};
arma::vec vc = {size1(1),size2(1),size3(1),size4(1)};
arma::Col<int> redu = {(int)round(size1(2)),(int)round(size2(2)),(int)round(size3(2)),(int)round(size4(2))};
Rcpp::DataFrame result = Rcpp::DataFrame::create(Rcpp::Named("N") = Nsize,
Rcpp::Named("Var.1") = v1,
Rcpp::Named("Var.C") = vc,
Rcpp::Named("Reduction") = redu);
return (result);
}
O código escrito acima é a versão em Rcpp da mesma função descrita no Jones da técnica de redução de variancia usando variável de controle. O código original em R se encontra no livro do Jones. Como saída, geramos um data frame como foi solicitado. As funções ‘Ginv’, ‘psi1’ e ‘psi2’ transformam os dados. A função ‘var’ estima as variancias e a redução utilizando a tecnica. Já a função ‘variancias’ gera o dataframe com as informações calculadas. Ao chamar a função variancias obtemos o dataframe que pode ser visto a seguir.
variancias()
N Var.1 Var.C Reduction
1 10 9.641807e-06 7.554533e-08 99
2 100 7.262149e-06 5.104257e-08 99
3 200 8.146245e-06 5.300702e-08 99
4 1000 8.176335e-06 1.063355e-07 99