這裡介紹如何使用 C++11 標準中內建的亂數函式庫，產生各種機率分布的隨機亂數。

傳統上在 C/C++ 程式中若要產生亂數，大家最常用的就是標準的 rand 函數，它的用法簡單、快速又方便，雖然其功能比較陽春，生成的亂數品質也比較不好，但因為是標準的函數之一，所以還是非常多人喜歡用。

在 C++11 的標準中新增了 <random> 這個標準的亂數函式庫，它的功能比較完整，生成的隨機變數品質也會比 rand 函數來得好，以下是簡單的使用教學與範例程式碼。

<random> 概念

C++11 的 <random> 的主要功能可分為三大類：

  • 隨機亂數種子產生器：傳統上許多人都會習慣使用時間當作亂數種子，但其實時間的資訊並不是隨機的，<random> 中的隨機設備（std::random_device）可以讓我們產生接近隨機的亂數種子，讓生成的亂數無法被預測。
  • 亂數產生器：這個部分就是所謂的偽隨機亂數產生器（pseudorandom number generator），可依照給定的亂數種子產生一連串的亂數，<random> 中提供了好幾種不同的亂數產生器，例如 std::mt19937、std::mt19937_64 等，另外也有最陽春的 std::default_random_engine。
  • 機率分佈：機率分布的功能就是將亂數產生器所產生的均勻分布亂數，轉換為各種統計上常用的機率分布，例如常態分布（std::normal_distribution）、指數分布（std::exponential_distribution）、卜瓦松分布（std::poisson_distribution）等。

將以上三大類的工具組合之後，就可以很輕易的產生各種機率分布的隨機亂數。

Hello World

這是使用 <random> 亂數函式庫產生 [0, 1) 區間連續型均勻分布（uniform distribution）的範例程式：

#include <iostream>
#include <random>   /* 亂數函式庫 */
#include <ctime>
int main(){
  /* 亂數產生器 */
  std::default_random_engine generator( time(NULL) );

  /* 亂數的機率分布 */
  std::uniform_real_distribution<float> unif(0.0, 1.0);

  /* 產生亂數 */
  float x = unif(generator);

  std::cout << "x = " << x << std::endl;
  return 0;
}

這裡我們使用 std::default_random_engine 這個預設的亂數產生器，以時間的資訊作為亂數種子，然後透過 std::uniform_real_distribution 轉換成 [0, 1) 區間的連續型均勻分布。

這個 hello world 程式只是示範最簡單的亂數產生方式，這樣所產生的亂數品質就跟傳統上的 rand 函數差不多（一樣差），接下來我們會使用隨機設備產生更接近隨機的亂數種子，並且將 std::default_random_engine 換成更好一點的亂數產生器，最後套上各種機率分布的轉換，生成各種不同分布的亂數。

隨機設備

一般的亂數產生器都會需要指定亂數種子，如果亂數種子不是隨機的，那麼所產生出來的亂數也就不是隨機的，嚴格來說就算使用時間當作亂數種子，其所產生的亂數也不是隨機的，這種亂數在某些情況下甚至可以被預測。

隨機設備（std::random_device）是一個均勻分布（uniform distribution）的整數亂數產生器，可以產生比較接近隨機的亂數（例如從硬體上取得隨機的資料），讓人無法預測，以下是 std::random_device 的一些基本用法。

#include <iostream>
#include <random>
int main(){
  /* 隨機設備 */
  std::random_device rd;

  /* 隨機亂數的範圍 */
  std::cout << "Min = " << rd.min()
          << ", Max = " << rd.max() << std::endl;

  /* 產生隨機的亂數 */
  std::cout << "Random Number = " << rd() << std::endl;

  /* 隨機設備的熵值 */
  std::cout << "Entropy = " << rd.entropy() << std::endl;

  return 0;
}

由於 std::random_device 在產生亂數時會消耗掉系統的熵（entropy），所以當熵耗盡時就會影響到亂數的產生，不適合用來產生大量的亂數。

熵（entropy）的意思可以想像成系統上的「隨機性資訊總量」，電腦可以靠著一些有隨機性的設備來產生隨機的亂數，例如系統行程 ID、滑鼠移動軌跡、網路封包等，但是這些設備是有限的，而其所產生出來的隨機資訊也只能使用一次，若使用第二次就不是隨機的，因此系統上的隨機性資訊總量有一定的值，消耗太快而耗盡的話，就要等待這些設備再度產生。

實務上我們可以使用 std::random_device 所產生的亂數來做為亂數種子，這樣產生的亂數會比較有隨機性：

#include <iostream>
#include <random>
int main(){
  /* 隨機設備 */
  std::random_device rd;

  /* 亂數產生器 */
  std::default_random_engine generator( rd() );

  /* 亂數的機率分布 */
  std::uniform_real_distribution<float> unif(0.0, 1.0);

  /* 產生亂數 */
  float x = unif(generator);

  std::cout << "x = " << x << std::endl;
  return 0;
}

亂數產生器

偽隨機亂數產生器會依照給定的亂數種子，「計算」出一連串類似隨機的亂數，<random> 中提供好幾種不同演算法的亂數產生器：

亂數產生器說明
std::minstd_rand0Minimal standard，x = x * 16807 % 2147483647
std::minstd_randMinimal standard，x = x * 48271 % 2147483647
std::mt19937梅森旋轉演算法
std::mt19937_64梅森旋轉演算法
std::ranlux24_baseSubtract with carry
std::ranlux48_baseSubtract with carry
std::ranlux2424-bit RANLUX generator by Martin Lüscher and Fred James, 1994
std::ranlux4848-bit RANLUX generator by Martin Lüscher and Fred James, 1994
std::knuth_bKnuth-B generator
std::default_random_engine預設亂數產生器，通常是 LCG 類型的。

各種亂數產生器的用法都類似，只要將前面範例中的 std::default_random_engine 抽換掉即可：

#include <iostream>
#include <random>
int main(){
  std::random_device rd;

  /* 梅森旋轉演算法 */
  std::mt19937 generator( rd() );

  std::uniform_real_distribution<float> unif(0.0, 1.0);
  float x = unif(generator);
  std::cout << "x = " << x << std::endl;
  return 0;
}

其他的用法以此類推，若要查詢亂數的範圍，也是使用 min 與 max 兩個方法函數：

/* 產生亂數的範圍 */
std::cout << "Min = " << generator.min()
        << ", Max = " << generator.max() << std::endl;

講究亂數品質的人，請避免使用 std::minstd_rand0、std::minstd_rand 與 std::default_random_engine 這些太簡略的演算法，可以改用 std::mt19937 這類比較新的演算法。

機率分佈

亂數產生器只能產生固定範圍的整數亂數，若我們需要其他類型的亂數，就要靠機率分佈（probability distribution）的函數將基本的離散型均勻分布轉換成我們需要的分布。

以下是 <random> 中所提供的機率分布函數。

機率分布函數說明
std::uniform_int_distribution離散型均勻分布
std::uniform_real_distribution連續型均勻分布
std::bernoulli_distribution伯努利分布
std::binomial_distribution二項分布
std::geometric_distribution幾何分布
std::negative_binomial_distribution負二項分佈
std::poisson_distribution卜瓦松分布
std::exponential_distribution指數分布
std::gamma_distribution伽瑪分布
std::weibull_distribution韋伯分布
std::extreme_value_distributionGeneralized extreme value distribution
std::normal_distribution常態分布
std::lognormal_distribution對數常態分布
std::chi_squared_distribution卡方分布
std::cauchy_distribution柯西分布
std::fisher_f_distributionF-分布
std::student_t_distribution學生t分佈
std::discrete_distribution自訂離散分布
std::piecewise_constant_distribution自訂 piecewise constant 分布
std::piecewise_linear_distribution自訂 piecewise linear 分布

機率分布函數的使用也是非常單純，將前面範例中的 std::uniform_real_distribution 置換成自己需要的分布即可，另外再指定機率分布的參數即可：

#include <iostream>
#include <random>
int main(){
  std::random_device rd;
  std::mt19937 generator( rd() );

  /* 標準常態分布 */
  std::normal_distribution<double> norm(0.0, 1.0);

  /* 產生標準常態分布的亂數 */
  double x = norm(generator);
  std::cout << "x = " << x << std::endl;
  return 0;
}

常態分布需要的參數是平均值與標準差，這裡我們產生的是標準常態分布（平均值為 0、標準差為 1）的的亂數。

以下是一個產生常態分布亂數，並畫出簡易直方圖的範例：

#include <iostream>
#include <random>
int main(){
  std::random_device rd;
  std::mt19937 generator( rd() );

  /* 常態分布 */
  std::normal_distribution<double> norm(5.0, 2.0);

  /* 常態分布亂數直方圖 */
  const int n = 10000;
  const int nstars = 100;
  int p[10]={};
  for (int i=0; i < n; ++i) {
    double x = norm(generator);
    if ((x >= 0.0) && (x < 10.0)) ++p[int(x)];
  }
  for (int i=0; i<10; ++i) {
    std::cout << i << ": ";
    std::cout << std::string(p[i]*nstars/n, '*') << std::endl;
  }
  return 0;
}

執行後的輸出會類似這樣：

0: *
1: ****
2: *********
3: ***************
4: *******************
5: *******************
6: **************
7: ********
8: ****
9: *

只要將分布函數替換後，稍微修改一下分布的參數，就可以畫出各種分布的直方圖，以下是指數分布的範例：

#include <iostream>
#include <random>
int main(){
  std::random_device rd;
  std::mt19937 generator( rd() );

  /* 指數分布 */
  std::exponential_distribution<double> exp(0.3);

  /* 指數分布亂數直方圖 */
  const int n = 10000;
  const int nstars = 100;
  int p[10]={};
  for (int i=0; i < n; ++i) {
    double x = exp(generator);
    if ((x >= 0.0) && (x < 10.0)) ++p[int(x)];
  }
  for (int i=0; i<10; ++i) {
    std::cout << i << ": ";
    std::cout << std::string(p[i]*nstars/n, '*') << std::endl;
  }
  return 0;
}

執行後的輸出為：

0: **************************
1: *******************
2: *************
3: **********
4: ********
5: *****
6: ****
7: ***
8: **
9: *

參考資料