Skip to content
Prev 2290 / 10988 Next

[Rcpp-devel] Rcpp equivalent of sample()?

For now I am using the following, which seems to work.  I didn't really need
all the functionality from sample(); I just needed to get a random index
(but with non-uniform probabilities).  It might be cleaner to use
std::accumulate in rcategorical, but I got this working first.  Also,
rcategorical requires a double[K]; how do I cast my probabilities argument P
as a doubles array?  I naively thought P would already be an array of
doubles, but rcategorical(P,K) doesn't compile.

Any thoughts/suggestions on how to improve on this first attempt?
Thanks in advance,
Chris

###

inc <-
'int rcategorical(double probs[],int K) {
  double total = 0;
  double cumpr = 0;
  int k = 0;
  for (k = 0; k < K; k++) {
    total += probs[k];
  }
  double u = unif_rand();
  for (k = 0; k < K; k++) {
    cumpr += probs[k] / total;
    if (u < cumpr) break;
  }
  return(k);
}'

src <- '
Rcpp:NumericVector probs(P);
int K = probs.size();
double pr[K];
for (int k = 0; k < K; k++) {  // copy elements over one-by-one?  Noob
status, but works. :-)
  pr[k] = probs[k];
}
int R = rcategorical(pr,K);
return wrap(R);
'
fun <-  cxxfunction(signature(P="numeric"),includes=inc,body = src,
plugin="Rcpp")
p <- c(.2,.3,.4,.1)
fun(p)
On Sat, May 14, 2011 at 6:19 PM, Dirk Eddelbuettel <edd at debian.org> wrote:

            
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20110514/e7cd3a4f/attachment-0001.htm>