几个随机数 C++

several random numbers c++

本文关键字:随机数 C++ 几个      更新时间:2023-10-16

我是一名物理学家,编写一个程序,涉及从高斯分布中提取几个(数十亿量级(随机数。我正在尝试使用 C++11。这些随机数的生成由一个操作分隔,该操作应该花费很少的时间。我最担心的是,我生成如此多的随机数,并且时间间隔如此之小,是否可能导致次优性能。我正在测试某些统计属性,这些属性在很大程度上依赖于数字随机性的独立性,因此,我的结果对这些问题特别敏感。我的问题是,对于我在代码中提到的数字类型(我实际代码的简化版本(,我是否做了明显(甚至微妙(错误的事情?

#include <random>
// Several other includes, etc.
int main () {
  int dim_vec(400), nStats(1e8);
  vector<double> vec1(dim_vec), vec2(dim_vec);
  // Initialize the above vectors, which are order 1 numbers.
  random_device rd;
  mt19937 generator(rd());
  double y(0.0);
  double l(0.0);
  for (int i(0);i<nStats;i++)
    {
      for (int j(0);j<dim_vec;j++)
        {
          normal_distribution<double> distribution(0.0,1/sqrt(vec1[j]));
          l=distribution(generator);
          y+=l*vec2[j];
        }
      cout << y << endl;
      y=0.0;
    }
}

允许normal_distribution具有状态。 对于这种特定的分布,通常每隔一次调用成对生成数字,并且在奇数调用中,返回第二个缓存的数字。 通过在每个调用上构造一个新的发行版,您将丢弃该缓存。

幸运的是,您可以通过调用不同的 normal_distribution::p aram_type 来"塑造"单个发行版:

 normal_distribution<double> distribution;
 using P = normal_distribution<double>::param_type;
 for (int i(0);i<nStats;i++)
    {
      for (int j(0);j<dim_vec;j++)
        {
          l=distribution(generator, P(0.0,1/sqrt(vec1[j])));
          y+=l*vec2[j];
        }
      cout << y << endl;
      y=0.0;
    }

我不熟悉std::normal_distribution的所有实现。 但是我为libc ++编写了一个。 因此,我可以肯定地告诉您,我对代码的轻微重写将产生积极的性能影响。 我不确定它会对质量产生什么影响,只能说我知道它不会降低它。

更新

关于 Severin Pappadeux 下面关于在发行版中一次生成数字对的合法性的评论: 请参阅 N1452,其中讨论并允许这种技术:

分布有时会存储来自其关联源的值 调用其运算符((的随机数。例如,一个常见的 生成正态分布随机数的方法是 检索两个均匀分布的随机数并计算两个 正态分布的随机数。为了重置 分布的随机数缓存到定义的状态,每个 分发具有重置成员功能。它应该在 每当交换或恢复其关联的引擎时分发。

关于优秀 HH 答案的一些想法

  1. 正态分布 (mu,sigma( 由正态分布 (0,1( 通过偏移和尺度生成:

N(mu, sigma( = mu + N(0,1(*sigma

如果你的平均值(MU(始终为零,你可以通过做一些事情来简化和加速(通过不添加0.0(你的代码,

比如
normal_distribution<double> distribution;
for (int i(0);i<nStats;i++)
{
  for (int j(0);j<dim_vec;j++)
    {
      l  = distribution(generator);
      y += l*vec2[j]/sqrt(vec1[j]);
    }
  cout << y << endl;
  y=0.0;
}
  1. 如果速度是最重要的,我会尝试在主 10^8 循环之外预先计算我所能计算的一切。是否可以预先计算 sqrt(vec1[j]( 以便节省 sqrt(( 调用?是否有可能将 vec2[j]/sqrt(vec1[j]( 作为单个向量?

  2. 如果无法预先计算这些向量,我会尝试节省内存访问。将 vec2[j] 和 vec1[j] 的片段放在一起可能有助于获取一个缓存行而不是两个。因此声明vector<pair<double,double>> vec12(dim_vec);并在采样中使用y+=l*vec12[j].first/sqrt(vec12[j].second)