C++11 內建亂數函式庫使用教學：隨機亂數產生器與機率分布

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

傳統上在 C/C++ 程式中若要產生亂數，大家最常用的就是標準的 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），所以當熵耗盡時就會影響到亂數的產生，不適合用來產生大量的亂數。

實務上我們可以使用 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::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_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: *

參考資料：GeeksforGeeks、GeeksforGeeks、StackOverflow、dreamincode、PCG