用C++ TR1 產生隨機數

用C++ TR1 產生隨機數

原作 ：[英文原文]
翻譯 ：Orbit(阿貓阿狗)

介紹

本文將介紹如何使用C++ TR1(C++ 標準委員會 Technical Report 1) 提供的隨機數產生功能
除了介紹基本的一致隨機數產生之外，還會介紹隨機樣本的可能分布，包括：柏努利分布，二項分布，指數分布，伽馬分布，幾何分布，常態分佈和泊松分布。我會指出一些針對特定分布需要注意的地方，比如參數約定，最後還會給出一些提示，比如如何使用TR1不直接支援的分布產生隨機數，比如柯西分布，chi-squared分布和Student t。

Visual Studio 2008 現在通過

feature pack支援TR1擴充（orbit註：微軟已經發布了Visual Studio 2008的Services Pack 1，它包含了此前發布的feature pack，以及完整的TR1支援），其它支援TR1的軟體或庫有
Boost和Dinkumware。
GCC 從4.3版本開始支援 C++ TR1。

為了清楚起見，本文的例子代碼都使用了完整的名字空間(namespace)限定，你可以使用using語句來消除這些限定。想要瞭解更多的例子代碼可以看看這篇文章
測試 C++ TR1 隨機數產生機制 ，此文還介紹了一個不太成熟的測試架構。

 

提綱

• 開始
• 標頭檔和名字空間(namespace)
• 核心引擎
• 設定種子
• 從均勻分布產生隨機數
• 離散整數
• 實數
• 非均勻分布產生隨機數
• 直接支援的分布
• 柏努利分布
• 二項分布
• 指數分布
• 伽馬分布
• 幾何分布
• 常態分佈(Gaussian)
• 泊松分布
• 其它分布
• 柯西分布(Cauchy)
• 卡方分布(Chi squared)
• T分布(Student t)

開始標頭檔和名字空間(namespace)

C++ 的隨機數產生類和函數都在 <random> 標頭檔中定義，並包含在名字空間(namespace)std::tr1中，注意在C++中
tr是小寫字母，在英語文本中“TR”是大寫的。

核心引擎

任何偽隨機數產生軟體的核心就是一個產生均勻分布的隨機整數的常式，它被用在引導過程中產生均勻分布的實數，這些均勻分布的實數再通過各種變換（演算法）產生其它分布，例如“接受-拒絕”演算法，不過這些都需要一個原始的隨機數做為起始。

在 C++ TR1 中，有幾種不同的（隨機數）產生核心，或稱之為“引擎”供你選擇，下面就是Visual
Studio 2008 feature pack所支援的四個引擎類。

• linear_congruential 迴圈使用公式 x(i) = (A * x(i-1) + C) mod M 線性疊加
• mersenne_twister 使用了
  
  馬特賽旋轉演算 演算法
• subtract_with_carry （帶進位線性同餘演算法）對整數演算法迴圈使用公式
  x(i) = (x(i - R) - x(i - S) - cy(i - 1)) mod M
• subract_with_carry_01 （改進的帶進位線性同餘演算法）對浮點數演算法迴圈使用公式
  x(i) = (x(i - R) - x(i - S) - cy(i - 1)) mod 1

設定種子

每種引擎都有一個 seed() 函數，它接受一個
unsigned long 型參數作為隨機數產生的種子。 也可以通過引擎之間獨特的模版參數設定種子。

 

從均勻分布產生隨機數離散整數

uniform_int 類以相同的機率在一個範圍內抽取整數，它的建構函式有兩個參數，分別表示抽取範圍的最大值和最小值，需要注意的是抽取範圍是個閉區間，也就是這兩個值也可能被抽取到。

例子:

下面面的代碼將列印5個從集合：1,2,3,...,52隨即抽取的數字，這可能就是一個從一副牌中抽取五張牌的演算法模型，每次抽取一張牌，之後歸還、洗牌，再抽取下一張牌。

std::tr1::uniform_int<int> unif(1, 52);
for (int i = 0; i > 5; ++i)
    std::cout << unif(eng) << std::endl;

實數

uniform_real 類用於從一個區間產生一個浮點數，它的建構函式有兩個參數，就是區間的兩個端點，預設值是從0到1。這個類的（隨機）值是從半開區間產生的，對於參數 a 和 b，隨機值產生的範圍是區間：[a, b)，換句話說，就是產生的值x滿足條件： a <= x < b。

下面的代碼從區間 [3, 7) 產生隨機數。

std::tr1::uniform_real<double> reif(3,7);

 非均勻分布產生隨機數

C++ TR1 庫還支援根據各種分布類（演算法）產生非均勻分布的隨機數，這些類通過
operator() 方法返回隨機樣本值。這個方法使用一個引擎類作為參數，例如下面的代碼就示範了使用常態分佈產生10個隨機數。

std::tr1::mt19937 eng;  // a core engine class
std::tr1::normal_distribution<double> dist; 
for (int i = 0; i < 10; ++i)
std::cout << dist(eng) << std::endl;

儘管 TR1 的規範和Visual Studio的文檔都明確表示分布類內建一個預設的模版參數，但事實是必須顯式指定這個參數，因此上面的 normal_distribution 類被聲明為帶一個 double 類型的傳回值，即使double 就是預設類型。

下面的例子都假設引擎類 eng 正在使用。

直接支援的分布柏努利分布

柏努利隨機值以機率 p 返回1， 以機率 1-p 返回0值，這個類從柏努利分布類bernoulli_distribution產生樣本，它的建構函式只有一個參數，p，這個參數的預設值是0.5。有趣的是bernoulli_distribution類傳回值類型不是 int，而是bool 類型， 所以它返回true表示機率返回false表示機率1-p，這是唯一的一個返回布爾值的分布類，也正式這個原因，它成為唯一的一個沒有傳回值類型模版參數的分布類。

例子:

下面的代碼假設已經定義了引擎 eng。
這個例子將在四分之一的時間列印 “six of one” ，剩下的時間列印“half dozen of another”。

std::tr1::bernoulli_distribution bern(0.25);
std::cout << (bern(eng) ? "six of one" : "half dozen of another") << std::endl;

二項分布

二項分布將n次獨立柏努利測試中成功的次數作為傳回值，成功的機率是p，產生二項分布隨機數的類是
binomial_distribution，它代兩個模版參數，就是n的類型和p的類型，預設值分別是int 和 double，建構函式使用參數 n = 1 和 p = 0.5做為預設值。

例子:

下面的代碼將類比擲5次骰子，統計骰子上面是1的次數。（所有的結果都很相似，你可以用來計算其它值出現的次數）

std::tr1::binomial_distribution<int, double> roll(1.0/6.0);
std::cout << roll(eng) << std::endl;

指數分布

指數分布有兩個公用的參數表示: 斜率和平均值，注意C++ TR1 標準已經指定了斜率的參數表示，使用指定的參數表示，密度函數就可以表示為 f(x) = λ e-λ x ，並且平局值就是 1/λ。

用於產生指數分布的類是
exponential_distribution，如果沒有指定斜率，建構函式就使用預設值 λ = 1。

例子:

下面的代碼就是從平均值是13的指數分布中找到10個值，並求和。

const double mean = 13.0;
double sum = 0.0;
std::tr1::exponential_distribution<double< exponential(1.0/mean);
for (int i = 0; i < 10; ++i)
    sum += exponential(eng);

伽馬分布

伽馬分布在 TR1 已經被參數化，它的密度函數是f(x) = xα-1e-x/Γ(α)。
α 參數就是形態參數，伽馬分布的更普遍的參數化表示應該包含比例參數 β，這樣伽馬分布的密度函數實際上是 f(x) = (x/β)α-1e-x/β/(βΓ(α))，TR1的參數化實際上就是令β = 1而已。

還有一個不太普遍的參數化表示是使用參數1/β，我猜想標準委員會接受只有一個參數的伽馬分布版本是為了避免與第二個參數混淆，除了這個動機之外，沒有理由要採用只有一個參數的版本。要使用比例參數 β 產生伽馬分布的隨機樣本，只要使用比例參數1產生伽馬分布的隨機樣本，然後乘以β就行了。

用於產生伽馬分布的類是 gamma_distribution，建構函式有一個參數，就是形態參數 α，這個參數預設是1。

例子:

下面的代碼用形態參數3和比例1產生10個伽馬分布的隨機數。

std::tr1::gamma_distribution<double> gamma(3.0);
for (int i = 0; i < 10; ++i)
    std::cout << gamma(eng) << std::endl;

例子:

下面的代碼用形態參數3和比例5產生10個伽馬分布的隨機數。

std::tr1::gamma_distribution<double> gamma(3.0);
for (int i = 0; i < 10; ++i)
    std::cout << 5.0*gamma(eng) << std::endl;

幾何分布

幾何分布有一個參數p，使用這個p作為柏努利測試的成功機率，連續進行柏努利測試直到第一次出現成功時統計之前所進行的柏努利測試的次數。

產生幾何分布的類是 geometric_distribution，它的建構函式有一個參數，就是成功的機率p。

例子:

下面的代碼示範如何使用十分之一的成功機率擷取一個幾何分布的隨機樣本。

std::tr1::geometric_distribution<int, double> geometric(0.1);
std::cout << geometric(eng) << std::endl;

常態分佈 (高斯分布)

常態分佈 (高斯分布)有兩個參數，平局值 μ 和標準差 σ。

產生常態分佈的類是
normal_distribution，它的建構函式有兩個參數， μ 和 σ，預設值分別是 0 和 1。
注意常態分佈類使用的是標準差，不是方差

例子:

下面的代碼用平均值3和標準差4（如果是方差就是16）。

std::tr1::normal_distribution<double> normal(3.0, 4.0);
std::cout << normal(eng) << std::endl;

泊松分布

泊松分布用於表示在一段給定的時間間隔內，各種可能的比率或獨立事件可能出現的機率。

用於產生泊松分布的類是
poisson_distribution，它的建構函式有一個參數，就是分布的平均值 λ 。

例子:

下面的代碼示範了用平均值7產生一個泊松分布的樣本

std::tr1::poisson_distribution<double> poisson(7.0);
std::cout << poisson(eng) << std::endl;

其它分布

C++ TR1 規範只是包含了一些常見的分布，其它分布可以通過直接支援的分布相互組合實現，下面的幾種分布就是這樣組合的例子，統計學書籍中可能會有更多的分布。

柯西分布(Cauchy)

要產生柯西分布的樣本，首先要產生(-π/2, π/2)區間上的平均分布樣本，然後取正切值（tangent）。

const double pi = 3.14159265358979;
std::tr1::uniform_real<double> uniform(-0.5*pi, 0.5*pi);
double u = uniform(eng);
std::cout << tan(u) << std::endl;

上面的代碼有些詭秘的地方需要說明一下，由於uniform_real 產生的值是一個半開區間，它就有可能返回區間最左端的值，正切函數(tangent)對於精確的值-π/2會溢出。然而，我們的值 pi 要比真的 π 值小那麼一點，正式這被階段的一點使得tan 在uniform 返回最小值時也沒有問題。

卡方分布(Chi-squared)

卡方分布 (χ2) 是伽馬分布的一種特殊形式，有 ν 自由度的卡方分布和形態參數是ν/2 ，比例參數是2 的伽馬分布是一樣的，由於C++ TR1的實現只是直接支援比例參數是1的伽馬分布，這就需要先產生形態參數是 ν/2，比例參數是1的伽馬分布，然後乘以2。

下面的代碼示範可使用5自由度的卡方分布。

double nu = 5.0;
std::tr1::gamma_distribution<double> gamma(0.5*nu);
std::cout << 2.0*gamma(eng) << std::endl;

T分布(Student t)

要產生 ν 自由度的T分布，首先產生一個區間(0,1)的常態分佈樣本 X 和一個ν 自由度的卡方分布樣本 Y，然後用 X/sqrt(Y/ν)計算出T分布樣本。

下面的代碼示範了使用7自由度的T分布產生樣本

double nu = 7.0;
std::tr1::normal_distribution<double> normal;
std::tr1::gamma_distribution<double> gamma(0.5*nu);
double x = normal(eng);
double y = 2.0*gamma(eng);
std::cout << x/sqrt(y/nu) << std::endl;