Inconsistent results between Rcpp and R code

Question

UPDATE Previous example is complicated, hence please allow me to use a simpler example as shown below:

Here is the Rcpp code:

#include <RcppArmadillo.h>
#include <RcppArmadilloExtensions/sample.h>
#include <Rmath.h>
#include <Rcpp.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp ;
using namespace arma;
using namespace std;

// [[Rcpp::export]]  
double chooseC(double n, double k) {
  return Rf_choose(n, k);
}
// [[Rcpp::export]]
double function3(double n, double m, double beta) {
  double prob;
  NumericVector k(m);
  NumericVector k_vec(m);
  if(n<m){prob=0;}
  else{
    if(chooseC(n,m)==R_PosInf){
      k=seq_len(m)-1;
      k_vec= (n-k)/(m-k)*std::pow((1-beta),(n-m)/m)*beta;
          prob=std::accumulate(k_vec.begin(),k_vec.end(), 1, std::multiplies<double>())*beta;
    }
    else{ 
      prob = beta * chooseC(n,m) * std::pow(beta,m) * std::pow((1-beta),(n-m));
    }

  }
  return(prob);
}

Here is the R code:

function4 <- function ( n , m , beta )
{
  if ( n < m )
  {
    prob <- 0.0
  }
  else
  {
    if (is.infinite(choose(n,m))){
      k<-0:(m-1)
      prob <- beta *prod((n-k)/(m-k)*(1-beta)^((n-m)/m)*beta)
    }
    else{
      prob <- beta * choose(n,m) * beta^m * (1-beta)^(n-m)
    }
  }
  prob
}

Comparison:

input<-619
beta<-0.09187495

x<-seq(0, (input+1)/beta*3)
yy<-sapply(x,function(n)function3(n,input, beta=beta))
yy2<-sapply(x,function(n)function4(n,input, beta=beta))
sum(yy)=0
sum(yy2)=1

However， with other input:

input<-1
beta<-0.08214248

Both results are the same, sum(yy)=sum(yy2)=0.9865887.

I used double in Rcpp code, I don't know what else could cause the inconsistent precision between Rcpp and R code.

Thanks a lot!

Hi! @coatless Thank you for your comment. I run that code in Rstudio, it returned the following error: > options(digits = 22) > sum(yy) [1] 0.019373457517486758 > sum(yy2) [1] 0.019373457517486748 Actually I expected the Rcpp returns nonzero value as R code does. Do I make any mistake in Rcpp code? Thank you very much! — Song Tang
– Song Tang, Commented May 21, 2017 at 22:36
There are only 16 digits of precision, so don't listen to people telling you to print at 22. The last six are random. — Dirk is no longer here
– Dirk is no longer here, Commented May 21, 2017 at 22:55
The question is far from minimal. Reduce it further, please. — Dirk is no longer here
– Dirk is no longer here, Commented May 21, 2017 at 23:06
Hi @DirkEddelbuettel, thank you for your comment. I updated my question, I used a simpler example. These are the cases where I found the Rcpp result is not consistent with the R code result. I appreciate for your help! — Song Tang
– Song Tang, Commented May 22, 2017 at 7:35

Song Tang · Accepted Answer · 2017-05-22 11:59:30Z

I think I fix the Rcpp code, so right now both Rcpp and R code produce the same results when the results are very small values. The solution is shown as below:

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

using namespace Rcpp ;
using namespace arma;
using namespace std;

// [[Rcpp::export]]  
double chooseC(double n, double k) {
  return Rf_choose(n, k);
}

// [[Rcpp::export]]
double function3(double n, double m, double beta) {
  double prob;
  arma::vec k = arma::linspace<vec>(0, m-1, m);
  arma::vec k_vec;

  if(n<m){prob=0;}
  else{
    if(chooseC(n,m)==R_PosInf){ 
      k_vec= (n-k)/(m-k)*pow((1-beta),(n-m)/m)*beta;  
      prob=arma::prod(k_vec)*beta;
    }
    else{ 
      prob = beta * chooseC(n,m) * pow(beta,m) * pow((1-beta),(n-m));
    }

  }
  return(prob);
}

However, I still do not understand why by writing code in this way will fix the precision inconsistent. Rcpp and RcppArmadillo still look like black boxes to me.

Collectives™ on Stack Overflow

Inconsistent results between Rcpp and R code

1 Answer 1

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Comments

Your Answer

Sign up or log in

Post as a guest

Related