7

I'm having some trouble with seeding a user defined RNG in R. It seems that

set.seed(123, kind='user', normal.kind='user')

Does not actually pass 123 to the user defined RNG initialization.

I went back to the documentation available at ?Random.user and tried the example code given there, with the minor modification that I print the seed passed to the user_unif_init function (full code below).

Steps to reproduce:

  1. Paste code below in urand.c
  2. Run R CMD SHLIB urand.c
  3. Open R
  4. Run the following commands:

    > dyn.load('urand.so')
    > set.seed(123, kind='user', normal.kind='user')
    Received seed: 720453763
    Received seed: 303482705 // any other numbers than 123
    

Here's the full code I used in urand.c:

// ##  Marsaglia's congruential PRNG

#include <stdio.h>
#include <R_ext/Random.h>

static Int32 seed;
static double res;
static int nseed = 1;

double * user_unif_rand()
{
    seed = 69069 * seed + 1;
    res = seed * 2.32830643653869e-10;
    return &res;
}

void  user_unif_init(Int32 seed_in) {
    printf("Received seed: %u\n", seed_in);
    seed = seed_in;
}
int * user_unif_nseed() { return &nseed; }
int * user_unif_seedloc() { return (int *) &seed; }

/*  ratio-of-uniforms for normal  */
#include <math.h>
static double x;

double * user_norm_rand()
{
    double u, v, z;
    do {
        u = unif_rand();
        v = 0.857764 * (2. * unif_rand() - 1);
        x = v/u; z = 0.25 * x * x;
        if (z < 1. - u) break;
        if (z > 0.259/u + 0.35) continue;
    } while (z > -log(u));
    return &x;
}

Any help would be greatly appreciated!

2 Answers 2

2

It appears that R scrambles the user supplied seed in RNG.c as follows:

for(j = 0; j < 50; j++)
    seed = (69069 * seed + 1)

(link to source)

Trying to unscramble this would be a way to get the original seed back.

UPDATE

Unscrambling can be done through the multiplicative inverse of 69069 as follows:

Int32 unscramble(Int32 scrambled)
{
        int j;
        Int32 u = scrambled;
        for (j=0; j<50; j++) {
                u = ((u - 1) * 2783094533);
        }
        return u;
}

Plugging this in my user_unif_init() function solves the problem.

2
  • 1
    This would seem to invite undefined behaviour.
    – Hugh
    Commented May 16 at 23:46
  • @Hugh Indeed. On should replace 2783094533 with 2783094533UL. Commented May 29 at 7:36
2

The seed that is forwarded to the RNG is different from the provided seed, however, it is reproducible when a "normal" workflow is used. This then gives reproducible random numbers:

dyn.load('urand.so')
RNGkind("user", "user")
#> Received seed: 1844983443
set.seed(123)
#> Received seed: 303482705
runif(10)
#>  [1] 0.42061954 0.77097033 0.14981063 0.27065365 0.77665767 0.96882090
#>  [7] 0.49077135 0.08621131 0.52903479 0.90398294
set.seed(123)
#> Received seed: 303482705
runif(10)
#>  [1] 0.42061954 0.77097033 0.14981063 0.27065365 0.77665767 0.96882090
#>  [7] 0.49077135 0.08621131 0.52903479 0.90398294

(Note that I have changed your urand.c slightly to use Rprintf from R_ext/Print.h.)


Edit: If you need control over the seed (why?), than you can do it yourself: replace user_unif_init, user_unif_nseed and user_unif_seedloc with

void set_seed(int * seed_in) {
    Rprintf("Received seed: %u\n", *seed_in);
    seed = *seed_in;
}

And call it explicitly:

dyn.load('urand.so')
RNGkind("user", "user")
set_seed <- function(seed) {
  invisible(.C("set_seed", seed_in = as.integer(seed)))
}
set_seed(123)
#> Received seed: 123
runif(10)
#>  [1] 0.00197801 0.61916849 0.34846373 0.04152509 0.09669026 0.29923760
#>  [7] 0.04184693 0.32557942 0.44473242 0.22339845
set_seed(123)
#> Received seed: 123
runif(10)
#>  [1] 0.00197801 0.61916849 0.34846373 0.04152509 0.09669026 0.29923760
#>  [7] 0.04184693 0.32557942 0.44473242 0.22339845

Edit 2: Have alook into the source at https://svn.r-project.org/R/trunk/src/main/RNG.c:

static void RNG_Init(RNGtype kind, Int32 seed)
{
    int j;

    BM_norm_keep = 0.0; /* zap Box-Muller history */

    /* Initial scrambling */
    for(j = 0; j < 50; j++)
    seed = (69069 * seed + 1);
    [...]

These 50 LCG rounds are responsible for the difference. My guess is that the authors of R assume that typical user supplied seeds are small and are therefore not random enough for a seed.

2
  • Thanks for taking a look! It's nice to know that it's consistent, but this doesn't really solve my problem. I need to be able to get 123 in the C function. Is there any way to bypass or undo the black magic that R inserts into this seemingly straighforward process?
    – GjjvdBurg
    Commented Dec 15, 2017 at 13:44
  • Thanks for the update! I need control over the seed for a specific application. Your solution is a nice work-around, but it doesn't quite answer my question. I wonder what R is doing to the seed and for what purpose it would be useful.
    – GjjvdBurg
    Commented Dec 15, 2017 at 19:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.