[Rcpp-devel] Difference between runif and unif_rand?

[Dirk Eddelbuettel]

On 12 July 2013 at 12:28, Neal Fultz wrote:
| 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();

There was a reason we never ever used this one; I think it is meant for
standalone use via Rmathlib.
 
|   do { u -= p[i++]; } while(u  > 0);
| 
|   return i;
| }
| ')
| 
| 
| cppFunction('
| int f_broke(NumericVector p) {
| 
|   int i = 0;
|   
|   double u = Rcpp::runif(1)[0];

If you really want just oone draw, use R::runif(...)

If you set / reset the RNG seed you very clearly get the exact same draws
from R and Rcpp.

So maybe it is just your assumption that unif_rand() and runif() should give
identical results?

Also, never hurts to do add an explicit RNGScope object even though
cppFunction() will add it for you too...

|   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?

I am not entirely sure.  My head is a little foggy froim traveling and I am
about to get onto another plane.  Maybe I'll take a look...

Hope this helps so far,  Dirk

| 
| Thanks much,
| 
| Neal
| 
| 
| 
| 
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com