Siew użytkownikowi dostarczana generatora liczb losowych w R

Question

Mam pewne problemy z zaszczepienie zdefiniowane przez użytkownika RNG w R. Wydaje się, że

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

faktycznie nie przechodzą 123 do inicjalizacji RNG zdefiniowane przez użytkownika.

Powróciłem do dokumentacji dostępnej pod numerem ?Random.user i wypróbowałem podany tam przykładowy kod, z drobną modyfikacją, którą wydrukuję, przekazując materiał siewny do funkcji user_unif_init (pełny kod poniżej).

Kroki do odtworzenia:

1. Wklej poniższy kod w urand.c
2. Run R CMD SHLIB urand.c
3. Otwórz R

4. uruchom następujące polecenia:

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

Oto pełny kod użyłem w 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; 
} 

Każda pomoc będzie bardzo mile widziane!

Answer 1

Wydaje się, że R szyfruje dostarczonego materiału siewnego użytkownika w RNG.c następująco:

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

(link to source)

Próbując rozszyfrować byłoby to sposób na uzyskanie oryginalnego materiału siewnego powrotem.

UPDATE

unscrambling można zrobić poprzez multiplicative inverse z 69069 w następujący sposób:

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

Podłączenie to w moim user_unif_init() funkcji rozwiązuje problem.

Answer 2

Nasiono, które jest przekazywane do RNG różni się od dostarczonego materiału siewnego, jednak jest odtwarzalne, gdy używany jest "normalny" przepływ pracy. To z kolei daje powtarzalne liczb losowych:

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 

(Zauważ, że zmieniłeś urand.c nieznacznie używać Rprintf z R_ext/Print.h.)

---

Edit: Jeśli zajdzie potrzeba kontroli nad nasion (dlaczego?), niż możesz zrobić to sam: zastąp user_unif_init, user_unif_nseed i user_unif_seedloc z

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

i nazywają to wyraźnie:

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: Czy alook do źródła w 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); 
    [...] 

te 50 rund LCG są odpowiedzialne za różnice. Domyślam się, że autorzy R zakładają, że typowe dostarczone przez użytkownika nasiona są małe i dlatego nie są przypadkowe dla nasion.