Посев пользовательского генератора случайных чисел в R

Question

У меня возникли проблемы с заполнением пользовательского RNG в R. Кажется, что

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

На самом деле не проходит 123 пользователь определил инициализацию RNG.

Я вернулся к документации, доступной на ?Random.user и попробовал приведенный там пример кода с незначительной модификацией, которую я распечатаю user_unif_init функция (полный код ниже).

Действия по воспроизведению:

1. Вставьте код ниже в urand.c
2. Бежать R CMD SHLIB urand.c
3. открыто R

4. Запустите следующие команды:

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

Вот полный код, который я использовал в 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;
}

Любая помощь будет принята с благодарностью!

Answer 1

Похоже, что R скремблирует предоставленное пользователем семя в RNG.c следующее:

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

( ссылка на источник)

Попытка расшифровать это было бы способом вернуть оригинальное семя.

ОБНОВИТЬ

Расшифровка может быть выполнена с помощью обратного умножения 69069 следующим образом:

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

Подключая это в моем user_unif_init() Функция решает проблему.

Answer 2

Начальное число, которое передается в ГСЧ, отличается от предоставленного начального числа, однако оно воспроизводимо при использовании "нормального" рабочего процесса. Это тогда дает воспроизводимые случайные числа:

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

(Обратите внимание, что я изменил ваш urand.c немного использовать Rprintf от R_ext/Print.h.)

---

Изменить: Если вам нужен контроль над семенем (почему?), Чем вы можете сделать это самостоятельно: заменить user_unif_init, user_unif_nseed а также user_unif_seedloc с

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

И назовите это явно:

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

---

Изменить 2: Посмотрите на источник по адресу 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);
    [...]

Эти 50 раундов LCG несут ответственность за разницу. Я предполагаю, что авторы R предполагают, что типичные семена, поставляемые пользователем, малы и, следовательно, не достаточно случайны для семян.