[Rcpp-devel] Difference between runif and unif_rand?

[Neal Fultz]

I've been updating the C in an r package to use Rcpp and started seeing 
so odd results from a function that samples from a discrete
distribution.

I think I've narrowed it down to the sugar runif. The two programs below
are identical except f_works uses unif_rand and f_broke goes through runif. 

---

library(Rcpp)

cppFunction('
int f_works(NumericVector p) {

  int i = 0;

  double u = unif_rand();

  do { u -= p[i++]; } while(u  > 0);

  return i;
}
')


cppFunction('
int f_broke(NumericVector p) {

  int i = 0;
  
  double u = Rcpp::runif(1)[0];
  
  do { u -= p[i++]; } while(u  > 0);
  
  return i;
}
')

---

When I run f_broke, sometimes there is a number that doesn't belong:

R>p <- 1:4 / 10;
R>table(replicate(10000, f_works(p)))

   1    2    3    4 
 961 2045 2984 4010 

R>table(replicate(10000, f_broke(p)))

0.427225098479539 0.458629200235009 0.687304416438565 0.735608365153894 0.978602966060862 
                1                 1                 1                 1                 1 

                1                 2                 3                 4  
             1053              2010              3021              3911  


I'm pretty stumped about why I consistently see odd numbers every couple
thousand draws, becuase runif just wraps unif_rand anyway. My
functions are declared as int, so it should be impossible to
return a double. 

What is happening here that I don't get?

Thanks much,

Neal