为< random>创建PRNG引擎.在C ++ 11中与PRNG匹配的结果在R中 [英] Creating a PRNG engine for <random> in C++11 that matches PRNG results in R
问题描述
这个问题是双重的.我正在将R脚本转换为C ++，使用L'Ecuyer组合多重递归生成器(CMRG)作为其引擎(尤其是MRG32k3a)，然后从间隔(0，1)的均匀分布中返回一个随机数. R中的一个最小示例如下所示:
This question is twofold. I am translating an R script into C++ that uses the L'Ecuyer combined multiple recursive generator (CMRG) as it's engine (in particular, MRG32k3a), which then returns a random number from the uniform distribution over the interval (0, 1). A minimal example in R is shown below:
seednum<-100 # set seed
set.seed(seednum, kind="L'Ecuyer-CMRG") # set RNG engine
runif(1) # set distribution
我希望能够在R脚本和C ++代码之间验证我的结果(因为生成的随机数只是开始).我发现在不同语言中具有相同种子的PRNG不一定会产生相同的结果(因为它们可能具有编译器可以自由指定的参数)，如SO帖子此处.也就是说，根据PRNG的特定实现，使用相同的种子，相同的引擎和相同的分布可能会导致不同的随机数.下面是R和C ++ 11之间的一个相关示例.在R中使用无所不在的Mersenne-Twister PRNG:
I want to be able to validate my results between the R script and the C++ code (as the random numbers generated are only the beginning). I have found that PRNG's with the same seeds across different languages do not necessarily produce the same result (since they may have parameters that the compiler is free to specify) as seen in the SO posts here and here. That is to say, using the same seed, the same engine, and the same distribution may result in different random numbers depending on the particular implementation of the PRNG. A pertinent example between R and C++11 is below. Using the ubiquitous Mersenne-Twister PRNG in R:
seednum<-100
set.seed(seednum, kind="Mersenne-Twister")
runif(1)
得出随机数0.3077661.在C ++ 11中做同样的事情:
Results in a random number of 0.3077661. Doing the same thing in C++11:
#include <iostream>
#include <random>
int main()
{
unsigned seed = 100;
std::mt19937 generator (seed);
std::uniform_real_distribution<double> distribution (0.0, 1.0);
std::cout << distribution(generator) << std::endl;
return 0;
}
得出随机数0.671156.我最初对这个结果感到困惑，但是以前的SO问题为我澄清了这一点(如上所述).似乎有一些参数传递给R中的MRG32k3a，我需要在C ++中复制这些参数才能生成相同的随机数.因此，第一个问题是，在哪里可以找到有关R中指定这些参数的MRG32k3a实现的文档?
Results in a random number of 0.671156. I was originally confused over this result, but previous SO questions clarified this for me (as linked above). It would appear that there are parameters being passed to MRG32k3a in R that I need to replicate in C++ in order to generate the same random numbers. The first question is thus, where can I find the documentation on the MRG32k3a implementation in R that specifies these parameters?
第二个问题涉及在C ++ 11中实现此生成器.此生成器未出现在此处找到，并在C中实现了MRG32k3a的示例，如下所示:
The second question regards implementing this generator in C++11. This generator does not appear in the list of pre-configured engine types within the <random> library of C++11 listed here. An example of MRG32k3a being implemented in C can be found here and is shown below:
/*
32-bits Random number generator U(0,1): MRG32k3a
Author: Pierre L'Ecuyer,
Source: Good Parameter Sets for Combined Multiple Recursive Random
Number Generators,
Shorter version in Operations Research,
47, 1 (1999), 159--164.
---------------------------------------------------------
*/
#include <stdio.h>
#define norm 2.328306549295728e-10
#define m1 4294967087.0
#define m2 4294944443.0
#define a12 1403580.0
#define a13n 810728.0
#define a21 527612.0
#define a23n 1370589.0
/***
The seeds for s10, s11, s12 must be integers in [0, m1 - 1] and not all 0.
The seeds for s20, s21, s22 must be integers in [0, m2 - 1] and not all 0.
***/
#define SEED 100
static double s10 = SEED, s11 = SEED, s12 = SEED,
s20 = SEED, s21 = SEED, s22 = SEED;
double MRG32k3a (void)
{
long k;
double p1, p2;
/* Component 1 */
p1 = a12 * s11 - a13n * s10;
k = p1 / m1;
p1 -= k * m1;
if (p1 < 0.0)
p1 += m1;
s10 = s11;
s11 = s12;
s12 = p1;
/* Component 2 */
p2 = a21 * s22 - a23n * s20;
k = p2 / m2;
p2 -= k * m2;
if (p2 < 0.0)
p2 += m2;
s20 = s21;
s21 = s22;
s22 = p2;
/* Combination */
if (p1 <= p2)
return ((p1 - p2 + m1) * norm);
else
return ((p1 - p2) * norm);
}
int main()
{
double result = MRG32k3a();
printf("Result with seed 100 is: %f\n", result);
return (0);
}
如前所述，我需要使用此生成器来创建一个可以馈入统一实分布的引擎.问题是我不知道该如何完成，而且我似乎在任何地方都找不到任何信息(除了知道引擎是类之外).是否有可用的C ++ 11资源可以帮助我完成此类任务?我不是在寻求解决问题的方法，而是在寻求一些可以帮助我自己实现这一目标的指针.
As previously noted, I need to use this generator to create an engine that can be fed into the uniform real distribution. The problem is I have no idea how this is done and I can't seem to find any information anywhere (aside from knowing that engines are classes). Are there any C++11 resources available that might help me in such a task? I am not asking for a solution to the problem, but rather pointers that would help me in implementing this myself.
推荐答案
因此，第一个问题是，在哪里可以找到有关R中指定这些参数的MRG32k3a实现的文档?
The first question is thus, where can I find the documentation on the MRG32k3a implementation in R that specifies these parameters?
我将使用源代码: https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/RNG.c#L143
问题是我不知道该如何完成，而且似乎无法在任何地方找到任何信息(除了知道引擎是类).
The problem is I have no idea how this is done and I can't seem to find any information anywhere (aside from knowing that engines are classes).
RandomNumberEngine的要求可以在这里找到:
https://en.cppreference.com/w/cpp/named_req/RandomNumberEngine
如果要使用uniform_real_distribution，虽然足以满足 UniformRandomBitGenerator :
The requirements for a RandomNumberEngine can be found here:
https://en.cppreference.com/w/cpp/named_req/RandomNumberEngine
Although it is sufficient to fulfill UniformRandomBitGenerator if you want to use uniform_real_distribution:
Expression Return type Requirements
G::result_type T T is an unsigned integer type
G::min() T Returns the smallest value that G's operator()
may return. The value is strictly less than
G::max().
G::max() T Returns the largest value that G's operator() may
return. The value is strictly greater than
G::min()
g() T Returns a value in the closed interval [G::min(),
G::max()]. Has amortized constant complexity.
主要问题是MRG32k3a打算返回(0,1)中的浮点数，而C ++ UniformRandomBitGenerator返回一个整数类型.为什么要与<random>标头集成?
Main problem is that MRG32k3a is meant to return a floating point number in (0,1), while a C++ UniformRandomBitGenerator returns an integer type. Why do you want to integrate with the <random> header?
您必须考虑的其他困难:
Additional difficulties you would have to take into account:
- R，c.f.使用的播种策略 https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/RNG.c#L293
- R扰乱用户提供的种子，请参见c.f. https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/RNG.c#L272
- Seeding strategy used by R, c.f. https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/RNG.c#L293
- R scrambles the user supplied seed, c.f. https://github.com/wch/r-source/blob/5a156a0865362bb8381dcd69ac335f5174a4f60c/src/main/RNG.c#L272
替代方法包括直接使用R源代码，而不与<random>标头集成或链接到libR.
Alternatives would include using R source code directly without integration with the <random> header or link to libR.
这篇关于为< random>创建PRNG引擎.在C ++ 11中与PRNG匹配的结果在R中的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
