在R中为用户提供的随机数生成器播种

Question

我在用R在R中播种用户定义的RNG时遇到了一些麻烦

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如下所示：

\n\n

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

\n\n

（链接到来源）

\n\n

尝试解读这将是找回原始种子的一种方法。

\n\n

更新

\n\n

可以通过69069 的乘法逆元来进行解密，如下所示：

\n\n

Int32 unscramble(Int32 scrambled)\n{\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82int j;\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82Int32 u = scrambled;\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82for (j=0; j<50; j++) {\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82u = ((u - 1) * 2783094533);\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82}\n\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82\xe2\x80\x82return u;\n}\n

\n\n

将其插入我的user_unif_init()函数中可以解决问题。

\n