[Rd] Bias in R's random integers?

Dirk Eddelbuettel edd @ending from debi@n@org
Fri Sep 21 16:20:00 CEST 2018


On 21 September 2018 at 09:50, Steve Grubb wrote:
| Hello,
| 
| Top posting. Several people have asked about the code to replicate my 
| results. I have cleaned up the code to remove an x/y coordinate bias for 
| displaying the results directly on a 640 x 480 VGA adapter. You can find the 
| code here:
| 
| http://people.redhat.com/sgrubb/files/vseq.c
| 
| To collect R samples:
| X <- runif(10000, min = 0, max = 65535)
| write.table(X, file = "~/r-rand.txt", sep = "\n", row.names = FALSE)
| 
| Then:
| cat ~/r-rand.txt | ./vseq > ~/r-rand.csv
| 
| And then to create the chart:
| 
| library(ggplot2);
| num.csv <- read.csv("~/random.csv", header=T)
| qplot(X, Y, data=num.csv);
| 
| Hope this helps sort this out.

Here is a simpler version in one file, combining the operations in a bit of
C++ and also an auto-executed piece of R to create the data, call the
transformation you created, and plot.  Store as eg /tmp/stevegrubb.cpp
and then call   Rcpp::sourceCpp("/tmp/stevegrubb.cpp")

Dirk

--- snip --------------------------------------------------------------------
#include <Rcpp.h>

#define SCALING 0.5

// [[Rcpp::export]]
Rcpp::NumericMatrix stevePlot(Rcpp::NumericVector itable) {
  size_t cnt = itable.size();
  Rcpp::NumericMatrix M(cnt,2);
  size_t num = 0;
  double nth, n1th = 0, n2th = 0, n3th = 0;
  double x, y;
  double usex, usey, pha = 0;

  
  usex = sin(pha);
  usey = cos(pha);
  while (num < cnt) {
    double x1, y1, z1;
    nth = itable[num];
    x1 = (nth - n1th)  * SCALING;
    y1 = (n1th - n2th) * SCALING;
    z1 = (n2th - n3th) * SCALING;
    x = (x1*usey+z1*usex);
    y = y1 + (z1*usey-x1*usex) / 2;
    if (num > 4) {
      M(num,0) = x;
      M(num,1) = y;
    } else {
      M(num,0) = M(num,1) = NA_REAL;
    }
    num++;
    n3th=n2th; n2th=n1th; n1th=nth;
  }
  Rcpp::colnames(M) = Rcpp::CharacterVector::create("X", "Y");
  return(M);
}

/*** R
library(ggplot2);
Xin <- runif(10000, min = 0, max = 65535)
Xout <- stevePlot(Xin);
qplot(X, Y, data=as.data.frame(Xout));
*/
--- snip --------------------------------------------------------------------


-- 
http://dirk.eddelbuettel.com | @eddelbuettel | edd using debian.org



More information about the R-devel mailing list