为什么新的随机库比std :: rand()更好？

Question

所以,我看到一个谈话叫兰特()是有害的,它提倡使用随机数生成的发动机分布模式与简单的std::rand()加模模式.

但是,我想看到std::rand()第一手的失败,所以我做了一个快速的实验:

1. 基本上,我编写了2个函数getRandNum_Old(),getRandNum_New()它们分别使用std::rand()和std::mt19937+ 生成0到5之间的随机数std::uniform_int_distribution.
2. 然后我使用"旧"方式生成了960,000(可被6整除)的随机数,并记录了数字0-5的频率.然后我计算了这些频率的标准偏差.我正在寻找的是尽可能低的标准差,因为如果分布是真正均匀的那样会发生什么.
3. 我运行了1000次模拟并记录了每次模拟的标准偏差.我还记录了它花费的时间,以毫秒为单位.
4. 之后,我再次完全相同,但这次生成随机数字的"新"方式.
5. 最后,我计算出的平均值和标准偏差为新老两种方式,平均和采取两种新老方法次名单的标准偏差的列表的标准偏差.

结果如下:

[OLD WAY]
Spread
       mean:  346.554406
    std dev:  110.318361
Time Taken (ms)
       mean:  6.662910
    std dev:  0.366301

[NEW WAY]
Spread
       mean:  350.346792
    std dev:  110.449190
Time Taken (ms)
       mean:  28.053907
    std dev:  0.654964

令人惊讶的是,两种方法的卷的总扩散是相同的.也就是说,std::mt19937+ std::uniform_int_distribution并不比简单std::rand()+ 更"统一" %.我做的另一个观察是新的比旧的方式慢了大约4倍.总的来说,似乎我付出了巨大的速度成本,几乎没有质量上的提升.

我的实验在某种程度上有缺陷吗？或者std::rand()真的不是那么糟糕,甚至可能更好？

作为参考,这是我完整使用的代码:

#include <cstdio>
#include <random>
#include <algorithm>
#include <chrono>

int getRandNum_Old() {
    static bool init = false;
    if (!init) {
        std::srand(time(nullptr)); // Seed std::rand
        init = true;
    }

    return std::rand() % 6;
}

int getRandNum_New() {
    static bool init = false;
    static std::random_device rd;
    static std::mt19937 eng;
    static std::uniform_int_distribution<int> dist(0,5);
    if (!init) {
        eng.seed(rd()); // Seed random engine
        init = true;
    }

    return dist(eng);
}

template <typename T>
double mean(T* data, int n) {
    double m = 0;
    std::for_each(data, data+n, [&](T x){ m += x; });
    m /= n;
    return m;
}

template <typename T>
double stdDev(T* data, int n) {
    double m = mean(data, n);
    double sd = 0.0;
    std::for_each(data, data+n, [&](T x){ sd += ((x-m) * (x-m)); });
    sd /= n;
    sd = sqrt(sd);
    return sd;
}

int main() {
    const int N = 960000; // Number of trials
    const int M = 1000;   // Number of simulations
    const int D = 6;      // Num sides on die

    /* Do the things the "old" way (blech) */

    int freqList_Old[D];
    double stdDevList_Old[M];
    double timeTakenList_Old[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_Old, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_Old();
            freqList_Old[roll] += 1;
        }
        stdDevList_Old[j] = stdDev(freqList_Old, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_Old[j] = timeTaken;
    }

    /* Do the things the cool new way! */

    int freqList_New[D];
    double stdDevList_New[M];
    double timeTakenList_New[M];

    for (int j = 0; j < M; j++) {
        auto start = std::chrono::high_resolution_clock::now();
        std::fill_n(freqList_New, D, 0);
        for (int i = 0; i < N; i++) {
            int roll = getRandNum_New();
            freqList_New[roll] += 1;
        }
        stdDevList_New[j] = stdDev(freqList_New, D);
        auto end = std::chrono::high_resolution_clock::now();
        auto dur = std::chrono::duration_cast<std::chrono::microseconds>(end-start);
        double timeTaken = dur.count() / 1000.0;
        timeTakenList_New[j] = timeTaken;
    }

    /* Display Results */

    printf("[OLD WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_Old, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_Old, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_Old, M));
    printf("\n");
    printf("[NEW WAY]\n");
    printf("Spread\n");
    printf("       mean:  %.6f\n", mean(stdDevList_New, M));
    printf("    std dev:  %.6f\n", stdDev(stdDevList_New, M));
    printf("Time Taken (ms)\n");
    printf("       mean:  %.6f\n", mean(timeTakenList_New, M));
    printf("    std dev:  %.6f\n", stdDev(timeTakenList_New, M));
}

Answer 1

几乎任何"旧"的实现都rand()使用LCG ; 虽然它们通常不是最好的发电机,但通常你不会在这样的基本测试中看到它们失败 - 即使是最差的PRNG,平均值和标准偏差通常也是正确的.

"坏"的常见缺陷 - 但很常见 - rand()实现是:

• 低阶位的低随机性;
• 短期内;
• 低RAND_MAX;
• 连续提取之间的一些相关性(通常,LCG产生有限数量的超平面上的数字,尽管这可以以某种方式减轻).

但是,这些都不是特定于API的rand().一个特定的实现可以将xorshift-family生成器置于srand/ rand并且在算法上说,获得最先进的PRNG而不改变接口,因此没有像你所做的那样的测试会显示输出中的任何弱点.

编辑: @R.正确地注意到rand/ srand接口是由如下事实的限制srand需要一个unsigned int,所以任何发生器一种实施方式可以把它们后面固有地限于UINT_MAX可能起始种子(并且因此产生的序列).这的确也是如此,虽然API可以平凡扩展,使srand采取unsigned long long,或添加一个单独的srand(unsigned char *, size_t)过载.

实际上,实际问题原则上rand()并不是很多,而是:

• 向后兼容; 许多当前的实现使用次优的生成器,通常具有错误选择的参数; 一个臭名昭着的例子是Visual C++,其中RAND_MAX只有32767.但是,这不容易改变,因为它会破坏与过去的兼容性 - 使用srand固定种子进行可重复模拟的人不会太高兴(事实上,IIRC)上述实现可以追溯到Microsoft C早期版本 - 甚至是八十年代中期的Lattice C);

• 简单的界面; rand()为整个程序提供具有全局状态的单个生成器.虽然这对于许多简单的用例来说非常好(并且实际上非常方便),但它会带来问题:

  • 使用多线程代码:为了修复它,你需要一个全局互斥体 - 它会无缘无故地减慢一切并消除任何重复性的可能性,因为调用序列本身就是随机的 - 或线程本地状态; 最后一个已经被几个实现采用(特别是Visual C++);
  • 如果您想要一个"私有",可重现的序列到您的程序的特定模块中,不会影响全局状态.

最后,rand事态:

• 没有指定实际的实现(C标准只提供了一个示例实现),因此任何旨在跨不同编译器生成可重现输出(或期望某些已知质量的PRNG)的程序必须使用自己的生成器;
• 没有提供任何跨平台的方法来获得合适的种子(time(NULL)不是,因为它不够精细,而且经常 - 认为没有RTC的嵌入式设备 - 甚至不够随机).

因此新的<random>标头,试图修复这个混乱提供以下算法:

• 完全指定(因此您可以使用交叉编译器可重现的输出和保证的特性 - 例如,发电机的范围);
• 通常具有最先进的质量(从设计图书馆时开始 ;见下文);
• 封装在类中(因此不会强制使用全局状态,这可以避免完全线程化和非局部性问题);

...以及默认random_device以及播种它们.

现在,如果你问我,我也喜欢也建立在此之上的"易","猜一个数字"的情况下(类似Python不如何提供"复杂"的API,也是微不足道的一个简单的API random.randint&CO.使用全球预播种PRNG给我们那些不喜欢随机设备/引擎/适配器/无论何时我们想要为宾果卡提取数字的简单人群,但是你可以很容易地通过当前的设施自己构建它(虽然构建"完整"API而不是简单的API)是不可能的.

最后,为了回到你的性能比较:正如其他人所指出的那样,你正在将快速LCG与较慢(但通常被认为质量较好)的Mersenne Twister进行比较; 如果您对LCG的质量感到满意,可以使用std::minstd_rand而不是std::mt19937.

确实,经过调整你的函数使用std::minstd_rand并避免无用的静态变量进行初始化

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    static std::uniform_int_distribution<int> dist{0, 5};
    return dist(eng);
}

我得到9毫秒(旧)对比21毫秒(新); 最后,如果我摆脱dist(与经典的模运算符相比,处理输出范围的分布偏差而不是输入范围的多倍)并回到你正在做的事情getRandNum_Old()

int getRandNum_New() {
    static std::minstd_rand eng{std::random_device{}()};
    return eng() % 6;
}

我把它降到6毫秒(所以,快了30%),可能是因为,与调用不同rand(),std::minstd_rand更容易内联.

顺便说一下,我使用手动滚动(但几乎符合标准库接口)进行相同的测试XorShift64*,并且比rand()(3.68 ms vs 8.61 ms)快2.3倍; 因为,不像梅森倍捻机和各种提供的LCG,它通过目前的随机性测试套件出色 和它的速度极快,它使你想知道为什么它不是在标准库中包含的呢.

Answer 2

如果您使用大于5的范围重复实验,那么您可能会看到不同的结果.当您的范围明显小于RAND_MAX大多数应用程序时没有问题.

例如,如果我们有RAND_MAX25,那么rand() % 5将生成具有以下频率的数字:

0: 6
1: 5
2: 5
3: 5
4: 5

由于RAND_MAX保证大于32767并且最不可能和最可能之间的频率差异仅为1,对于小数字,对于大多数用例,分布几乎足够随机.