为什么C++和R中的这些RNG不会产生类似的结果？

Question

请原谅这篇文章令人作呕的noobish性质,但我对那些在个人电脑上使用C++和R编程的人有一个疑问.

问题:为什么下面两个程序产生的这些随机数不相等,我该如何解决这个问题？

1. 首先,我怀疑我在R程序中误用了local函数和<<-运算符.
2. 其次,我怀疑这可能是一个浮动准确性问题.这对我来说并不是很明显两个程序是如何不同的,所以我不知道如何解决这个问题.

我已经尝试将我在C++中的所有计算转换为double/ float(偶数long double),并使用fmod而不是模数运算符%:再次输出不同的输出,但仍然与R中的输出不相似.我不知道它是否具有任何重要意义,但我我想补充说我正在使用G ++编译器编译C++代码.

算法:以下算法可用于任何标准个人计算机.建议并行使用三个字生成器,

• m _k = 171 m _k-1(mod 30269)
• M ' _ķ 172例M' _K-1(MOD 30307)
• 米'' _ķ 172例M '' _k-1个(MOD 30323)

并用作伪随机数的小数部分

• 克_ķ = { 米_ķ/30269 + M" _ķ/30307 + M '' _ķ/30323}

我已经使用的初始值中号₀ = 5,M" ₀ = 11,和M '' ₀ = 17.

程序:我在C++中有以下程序:

//: MC:Uniform.cpp
// Generate pseudo random numbers uniformly between 0 and 1
#include <iostream>
#include <math.h> // For using "fmod()"
using namespace std;

float uniform(){
    // A sequence of initial values
    static int x = 5;
    static int y = 11;
    static int z = 17;

    // Some integer arithmetic required
    x = 171 * (x % 177) - 2 * (x / 177);
    y = 172 * (x % 176) - 35 * (y / 176);
    z = 170 * (x % 178) - 63 * (z / 178);

    /* If both operands are nonnegative then the
    remainder is nonnegative; if not, the sign of
    the remainder is implementation-defined. */
    if(x < 0)
        x = x + 30269;
    if(y < 0)
        y = y + 30307;
    if(z < 0)
        z = z + 30323;

    return fmod(x / 30269. + y / 30307. + z / 30323., 1.);
}

int main(){
    // Print 5 random numbers
    for(int i = 0; i < 5; i++){
        cout << uniform() << ", ";
    }               
}///:~

该程序使用代码进行扩展并输出以下内容:

0.686912, 0.329174, 0.689649, 0.753722, 0.209394,

我在R中也有一个程序,如下所示:

## Generate pseudo random numbers uniformly between 0 and 1
uniform <- local({
    # A sequence of initial values
    x = 5
    y = 11
    z = 17

    # Use the <<- operator to make x, y and z local static 
    # variables in R.
    f <- function(){
        x <<- 171 * (x %% 177) - 2 * (x / 177)
        y <<- 172 * (y %% 176) - 35 * (y / 176)
        z <<- 170 * (z %% 178) - 63 * (z / 178)

        return((x / 30269. + y / 30307. + z / 30323.)%%1.)
    }
})

# Print 5 random numbers
for(i in 1:5){
    print(uniform())
}

该程序也使用代码进行扩展并生成输出

[1] 0.1857093
[1] 0.7222047
[1] 0.05103441
[1] 0.7375034
[1] 0.2065817

任何建议都表示赞赏,提前谢谢.

Answer 1

你%/%的R代码需要更多(整数除法).请记住,默认情况下,R中的数值变量是浮点数,而不是整数.所以/将使用非整数商进行普通除法.您还省略了处理否定x/ y/的部分z.

f <- function(){
    x <<- 171 * (x %% 177) - 2 * (x %/% 177)
    y <<- 172 * (y %% 176) - 35 * (y %/% 176)
    z <<- 170 * (z %% 178) - 63 * (z %/% 178)

    if(x < 0)
        x <<- x + 30269;
    if(y < 0)
        y <<- y + 30307;
    if(z < 0)
        z <<- z + 30323;

    return((x / 30269. + y / 30307. + z / 30323.)%%1)
}

进行这些更改后,结果似乎没有任何严重错误.一个100000随机抽取的快速直方图看起来非常均匀,我找不到自相关.仍然不符合您的C++结果....

Answer 2

您的C++代码中存在简单的复制/粘贴错误.这个

x = 171 * (x % 177) -  2 * (x / 177);
y = 172 * (x % 176) - 35 * (y / 176);
z = 170 * (x % 178) - 63 * (z / 178);

应该是这个.

x = 171 * (x % 177) -  2 * (x / 177);
y = 172 * (y % 176) - 35 * (y / 176);
z = 170 * (z % 178) - 63 * (z / 178);