[Rcpp-devel] problem with sample() implemented in Rcpp/RcppArmadillo

[Christian Gunning]

Update -- I've implemented all of R's sample (with identical results)
except for the walker alias method.  I used Armadillo for ProbSampleReplace
and ProbSampleNoReplace, since there's no good STL replacement for R's
revsort.  I'll submit a patch to RcppArmadillo when I get a chance.  If
anyone wants this in the meantime, I can email the source as a package that
contains some example glue and unit tests.

best,
Christian

On Sat, Nov 3, 2012 at 4:52 AM, Christian Gunning <xian at unm.edu> wrote:

> Dear all,
>
> Code for uniform sampling is below, following the prototype from earlier
> post.  Results/unit-tests (not included) are identical to R-core sample().
>
> -Christian
>
> template <class T>
> T sample(const T &x, const unsigned int size, const bool replace ) {
>     int ii, jj;
>     int nOrig = x.size();
>     T ret(size);
>     if ( size > nOrig && !replace) throw std::range_error( "Tried to
> sample more elements than in x without replacement" ) ;
>     IntegerVector index(size);
>     if (replace) {
>         SampleReplace(index, nOrig);
>     } else {
>         SampleNoReplace(index, nOrig);
>     }
>
>     for (ii=0; ii<size; ii++) {
>         jj = index[ii];
>         ret[ii] = x[jj];
>     }
>     return(ret);
> }
>
> void SampleReplace( IntegerVector index, unsigned int nOrig) {
>     int ii;
>     int  nSample = index.size();
>     for (ii = 0; ii < nSample; ii++) {
>         index[ii] = nOrig * unif_rand();
>     }
> }
>
> void SampleNoReplace( IntegerVector index, unsigned int nOrig) {
>     int ii, jj;
>     int nSample = index.size();
>     IntegerVector sub(nOrig);
>     for (ii = 0; ii < nOrig; ii++) {
>         sub[ii] = ii;
>     }
>
>     for (ii = 0; ii < nSample; ii++) {
>         jj = nOrig * unif_rand();
>         index[ii] = sub[jj];
>         // replace sampled element with last, decrement
>         sub[jj] = sub[--nOrig];
>     }
> }
>



-- 
A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama!
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121103/10add7d7/attachment.html>