分布和内部状态

Question

在 Stackoverflow 上，有很多关于从先验未知范围生成均匀分布整数的问题。例如

• C++11 Generating random numbers from frequently changing range
• Vary range of uniform_int_distribution

典型的解决方案是这样的：

inline std::mt19937 &engine()
{
  thread_local std::mt19937 eng;
  return eng;
}

int get_int_from_range(int from, int to)
{
  std::uniform_int_distribution<int> dist(from, to);
  return dist(engine());
}

鉴于分发应该是一个轻量级对象并且多次重新创建它不存在性能问题，似乎即使是简单的分发也可能非常好并且通常会有some internal state。

所以我想知道是否通过不断重置它来干扰分布的工作方式（即在每次调用 get_int_from_range 时重新创建分布）我得到正确分布的结果。

Pete Becker 和 Steve Jessop 之间有一个long discussion，但没有最终决定。 在另一个问题（Should I keep the random distribution object instance or can I always recreate it?）中，内部状态的“问题​​”似乎不是很重要。

C++ 标准是否对此主题做出任何保证？

以下实现（来自N4316 - std::rand replacement）是否更可靠一些？

int get_int_from_range(int from, int to)
{
  using distribution_type = std::uniform_int_distribution<int>;
  using param_type = typename distribution_type::param_type;

  thread_local std::uniform_int_distribution<int> dist;
  return dist(engine(), param_type(from, to));    
}

编辑

这重用了分发的可能内部状态，但它很复杂，我不确定它是否值得麻烦：

int get_int_from_range(int from, int to)
{
  using range_t = std::pair<int, int>;
  using map_t = std::map<range_t, std::uniform_int_distribution<int>>;

  thread_local map_t range_map;

  auto i = range_map.find(range_t(from, to));
  if (i == std::end(range_map))
    i = range_map.emplace(
          std::make_pair(from, to),
          std::uniform_int_distribution<int>{from, to}).first;

  return i->second(engine());
}

^{（来自https://stackoverflow.com/a/30097323/3235496）}

Answer 1

有趣的问题。

所以我想知道是否会干扰分发的工作方式 不断重置它（即在每个 get_int_from_range 的调用）我得到了正确分布的结果。

我已经编写了代码来使用uniform_int_distribution 和poisson_distribution 进行测试。如果您愿意，可以很容易地扩展它来测试另一个发行版。答案似乎是是的。

样板代码：

#include <random>
#include <memory>
#include <chrono>
#include <utility>

typedef std::mt19937_64 engine_type;

inline size_t get_seed()
    { return std::chrono::system_clock::now().time_since_epoch().count(); }

engine_type& engine_singleton()
{  
    static std::unique_ptr<engine_type> ptr;

    if ( !ptr ) 
        ptr.reset( new engine_type(get_seed()) );
    return *ptr;
}

// ------------------------------------------------------------------------

#include <cmath>
#include <cstdio>
#include <vector>
#include <string>
#include <algorithm>

void plot_distribution( const std::vector<double>& D, size_t mass = 200 )
{
    const size_t n = D.size();
    for ( size_t i = 0; i < n; ++i ) 
    {
        printf("%02ld: %s\n", i, 
            std::string(static_cast<size_t>(D[i]*mass),'*').c_str() );
    }
}

double maximum_difference( const std::vector<double>& x, const std::vector<double>& y )
{
    const size_t n = x.size(); 

    double m = 0.0;
    for ( size_t i = 0; i < n; ++i )
        m = std::max( m, std::abs(x[i]-y[i]) );

    return m;
}

实际测试代码：

#include <iostream>
#include <vector>
#include <cstdio>
#include <random>
#include <string>
#include <cmath>

void compare_uniform_distributions( int lo, int hi )
{
    const size_t sample_size = 1e5;

    // Initialize histograms
    std::vector<double> H1( hi-lo+1, 0.0 ), H2( hi-lo+1, 0.0 );

    // Initialize distribution
    auto U = std::uniform_int_distribution<int>(lo,hi);

    // Count!
    for ( size_t i = 0; i < sample_size; ++i )
    {
        engine_type E(get_seed());

        H1[ U(engine_singleton())-lo ] += 1.0;
        H2[ U(E)-lo ] += 1.0;
    }

    // Normalize histograms to obtain "densities"
    for ( size_t i = 0; i < H1.size(); ++i )
    {
        H1[i] /= sample_size; 
        H2[i] /= sample_size; 
    }

    printf("Engine singleton:\n"); plot_distribution(H1);
    printf("Engine creation :\n"); plot_distribution(H2);
    printf("Maximum difference: %.3f\n", maximum_difference(H1,H2) );
    std::cout<< std::string(50,'-') << std::endl << std::endl;
}

void compare_poisson_distributions( double mean )
{
    const size_t sample_size = 1e5;
    const size_t nbins = static_cast<size_t>(std::ceil(2*mean));

    // Initialize histograms
    std::vector<double> H1( nbins, 0.0 ), H2( nbins, 0.0 );

    // Initialize distribution
    auto U = std::poisson_distribution<int>(mean);

    // Count!
    for ( size_t i = 0; i < sample_size; ++i )
    {
        engine_type E(get_seed());
        int u1 = U(engine_singleton());
        int u2 = U(E);

        if (u1 < nbins) H1[u1] += 1.0;
        if (u2 < nbins) H2[u2] += 1.0;
    }

    // Normalize histograms to obtain "densities"
    for ( size_t i = 0; i < H1.size(); ++i )
    {
        H1[i] /= sample_size; 
        H2[i] /= sample_size; 
    }

    printf("Engine singleton:\n"); plot_distribution(H1);
    printf("Engine creation :\n"); plot_distribution(H2);
    printf("Maximum difference: %.3f\n", maximum_difference(H1,H2) );
    std::cout<< std::string(50,'-') << std::endl << std::endl;

}

// ------------------------------------------------------------------------

int main()
{
    compare_uniform_distributions( 0, 25 );
    compare_poisson_distributions( 12 );
}

运行它here。

---

C++ 标准是否对此主题做出任何保证？

我不知道。但是，我会说该标准隐含地建议不要每次都重新创建引擎。对于任何发行版Distrib，Distrib::operator() 的原型采用引用 URNG& 而不是 const 引用。这是可以理解的，因为引擎可能需要更新其内部状态，但这也意味着代码看起来像这样

auto U = std::uniform_int_distribution(0,10);
for ( <something here> ) U(engine_type());

不编译，这对我来说是一个明确的动机，不要编写这样的代码。

---

我确信有很多关于如何正确使用随机库的建议。如果您必须处理使用random_devices 并允许确定性播种用于测试目的的可能性，它确实会变得复杂，但我认为将我自己的建议也提出来可能很有用：

#include <random>
#include <chrono>
#include <utility>
#include <functional>

inline size_t get_seed()
    { return std::chrono::system_clock::now().time_since_epoch().count(); }

template <class Distrib>
using generator_type = std::function< typename Distrib::result_type () >;

template <class Distrib, class Engine = std::mt19937_64, class... Args>
inline generator_type<Distrib> get_generator( Args&&... args )
{ 
    return std::bind( Distrib( std::forward<Args>(args)... ), Engine(get_seed()) ); 
}

// ------------------------------------------------------------------------

#include <iostream>

int main()
{
    auto U = get_generator<std::uniform_int_distribution<int>>(0,10);
    std::cout<< U() << std::endl;
}

运行它here。希望这会有所帮助！

编辑我的第一个建议是一个错误，对此我深表歉意；我们不能像上面的测试那样使用单例引擎，因为这意味着两个统一的 int 分布会产生相同的随机序列。相反，我依赖于std::bind 使用自己的种子在本地复制std::function 中新创建的引擎这一事实，这会产生预期的行为；具有相同分布的不同生成器产生不同的随机序列。