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:
urand.c
R CMD SHLIB urand.c
R
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!
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.
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With