[R] Code of sample() in C
Berwin A Turlach
berwin at maths.uwa.edu.au
Mon Apr 6 03:15:42 CEST 2009
G'day Ranjan,
On Sun, 5 Apr 2009 17:00:56 -0500
Ranjan Maitra <maitra at iastate.edu> wrote:
> [...[
> I haven't tried this function out much myself, so YMMV.
There is an obvious problem in this code since len is not decreased but
n is when an index is sampled. The last line in the last for loop
should be
x[j] = x[--len];
Otherwise this algorithm will produce samples in which an index can
be repeated. And the probabilities with which an index is sampled would
be non-trivial to determine; they would definitely not be uniform.
HTH.
Cheers,
Berwin
> #include <stdlib.h>
> #ifndef USING_RLIB
> #define MATHLIB_STANDALONE 1 /*It is essential to have this before
> the call to the Rmath's header file because this decides
> the definitions to be set. */
> #endif
> #include <Rmath.h>
>
> /* Note that use of this function involves a prior call to the Rmath
> library to get the seeds in place. It is assumed that this is done in
> the calling function. */
>
> /* Equal probability sampling; without-replacement case */
> /* Adapted from the R function called SampleNoReplace */
>
> /**
> * Stores k out of n indices sampled at random without replacement
> * in y.
> */
> int srswor(int n, int k, int *y)
> {
> if (k > n) {
> return 1;
> }
> else {
> const double len = (double) n;
> int i;
> int* x = malloc(n * sizeof(int));
> if (!x) {
> return 2;
> }
>
> for (i = 0; i < n; ++i) x[i] = i;
>
> for (i = 0; i < k; ++i) {
> const int j = (int)(len * runif(0.0, 1.0));
> y[i] = x[j];
> x[j] = x[--n];
> }
> free(x);
> }
> return 0;
> }
=========================== Full address =============================
Berwin A Turlach Tel.: +65 6516 4416 (secr)
Dept of Statistics and Applied Probability +65 6516 6650 (self)
Faculty of Science FAX : +65 6872 3919
National University of Singapore
6 Science Drive 2, Blk S16, Level 7 e-mail: statba at nus.edu.sg
Singapore 117546 http://www.stat.nus.edu.sg/~statba
More information about the R-help
mailing list