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

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

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


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

<random> 概念

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

隨機亂數種子產生器
傳統上許多人都會習慣使用時間當作亂數種子,但其實時間的資訊並不是隨機的,<random> 中的隨機設備(std::random_device)可以讓我們產生接近隨機的亂數種子,讓生成的亂數無法被預測。
亂數產生器
這個部分就是所謂的偽隨機亂數產生器(pseudorandom number generator),可依照給定的亂數種子產生一連串的亂數,<random> 中提供了好幾種不同的亂數產生器,例如 std::mt19937std::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_rand0 Minimal standardx = x * 16807 % 2147483647
std::minstd_rand Minimal standardx = x * 48271 % 2147483647
std::mt19937 梅森旋轉演算法
std::mt19937_64 梅森旋轉演算法
std::ranlux24_base Subtract with carry
std::ranlux48_base Subtract with carry
std::ranlux24 24-bit RANLUX generator by Martin Lüscher and Fred James, 1994
std::ranlux48 48-bit RANLUX generator by Martin Lüscher and Fred James, 1994
std::knuth_b Knuth-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;
}

其他的用法以此類推,若要查詢亂數的範圍,也是使用 minmax 兩個方法函數:

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

講究亂數品質的人,請避免使用 std::minstd_rand0std::minstd_randstd::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_distribution Generalized extreme value distribution
std::normal_distribution 常態分布
std::lognormal_distribution 對數常態分布
std::chi_squared_distribution 卡方分布
std::cauchy_distribution 柯西分布
std::fisher_f_distribution F-分布
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: *

參考資料:GeeksforGeeksGeeksforGeeksStackOverflowdreamincodePCG

程式設計

1 Comment

  1. -

    原來C++也有這麼方便的函式庫!
    謝謝你詳細的教學~

Leave a Reply